Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/system/System.hpp"
       5             : 
       6             : #include "pairinteraction/basis/BasisAtom.hpp"
       7             : #include "pairinteraction/basis/BasisPair.hpp"
       8             : #include "pairinteraction/enums/TransformationType.hpp"
       9             : #include "pairinteraction/interfaces/DiagonalizerInterface.hpp"
      10             : #include "pairinteraction/operator/OperatorAtom.hpp"
      11             : #include "pairinteraction/operator/OperatorPair.hpp"
      12             : #include "pairinteraction/system/SystemAtom.hpp"
      13             : #include "pairinteraction/system/SystemPair.hpp"
      14             : #include "pairinteraction/utils/eigen_assertion.hpp"
      15             : #include "pairinteraction/utils/eigen_compat.hpp"
      16             : 
      17             : #include <Eigen/SparseCore>
      18             : #include <complex>
      19             : #include <limits>
      20             : #include <memory>
      21             : #include <numeric>
      22             : #include <oneapi/tbb.h>
      23             : #include <optional>
      24             : #include <spdlog/spdlog.h>
      25             : 
      26             : namespace pairinteraction {
      27             : template <typename Derived>
      28         185 : System<Derived>::System(std::shared_ptr<const basis_t> basis)
      29         185 :     : hamiltonian(std::make_unique<typename System<Derived>::operator_t>(std::move(basis))) {}
      30             : 
      31             : template <typename Derived>
      32         414 : System<Derived>::System(const System &other)
      33         414 :     : hamiltonian(std::make_unique<typename System<Derived>::operator_t>(*other.hamiltonian)),
      34         414 :       hamiltonian_requires_construction(other.hamiltonian_requires_construction),
      35         414 :       hamiltonian_is_diagonal(other.hamiltonian_is_diagonal),
      36         828 :       blockdiagonalizing_labels(other.blockdiagonalizing_labels) {}
      37             : 
      38             : template <typename Derived>
      39          60 : System<Derived>::System(System &&other) noexcept
      40          60 :     : hamiltonian(std::move(other.hamiltonian)),
      41          60 :       hamiltonian_requires_construction(other.hamiltonian_requires_construction),
      42          60 :       hamiltonian_is_diagonal(other.hamiltonian_is_diagonal),
      43         120 :       blockdiagonalizing_labels(std::move(other.blockdiagonalizing_labels)) {}
      44             : 
      45             : template <typename Derived>
      46           0 : System<Derived> &System<Derived>::operator=(const System &other) {
      47           0 :     if (this != &other) {
      48           0 :         hamiltonian = std::make_unique<operator_t>(*other.hamiltonian);
      49           0 :         hamiltonian_requires_construction = other.hamiltonian_requires_construction;
      50           0 :         hamiltonian_is_diagonal = other.hamiltonian_is_diagonal,
      51           0 :         blockdiagonalizing_labels = other.blockdiagonalizing_labels;
      52             :     }
      53           0 :     return *this;
      54             : }
      55             : 
      56             : template <typename Derived>
      57           0 : System<Derived> &System<Derived>::operator=(System &&other) noexcept {
      58           0 :     if (this != &other) {
      59           0 :         hamiltonian = std::move(other.hamiltonian);
      60           0 :         hamiltonian_requires_construction = other.hamiltonian_requires_construction;
      61           0 :         hamiltonian_is_diagonal = other.hamiltonian_is_diagonal,
      62           0 :         blockdiagonalizing_labels = std::move(other.blockdiagonalizing_labels);
      63             :     }
      64           0 :     return *this;
      65             : }
      66             : 
      67             : template <typename Derived>
      68         659 : System<Derived>::~System() = default;
      69             : 
      70             : template <typename Derived>
      71           0 : const Derived &System<Derived>::derived() const {
      72           0 :     return static_cast<const Derived &>(*this);
      73             : }
      74             : 
      75             : template <typename Derived>
      76         249 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_basis() const {
      77         249 :     if (hamiltonian_requires_construction) {
      78         134 :         construct_hamiltonian();
      79         134 :         hamiltonian_requires_construction = false;
      80             :     }
      81         249 :     return hamiltonian->get_basis();
      82             : }
      83             : 
      84             : template <typename Derived>
      85          48 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_eigenbasis() const {
      86          48 :     if (hamiltonian_requires_construction) {
      87           0 :         construct_hamiltonian();
      88           0 :         hamiltonian_requires_construction = false;
      89             :     }
      90          48 :     if (!hamiltonian_is_diagonal) {
      91           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
      92             :     }
      93          48 :     return hamiltonian->get_basis();
      94             : }
      95             : 
      96             : template <typename Derived>
      97         162 : Eigen::VectorX<typename System<Derived>::real_t> System<Derived>::get_eigenenergies() const {
      98         162 :     if (hamiltonian_requires_construction) {
      99           0 :         construct_hamiltonian();
     100           0 :         hamiltonian_requires_construction = false;
     101             :     }
     102         162 :     if (!hamiltonian_is_diagonal) {
     103           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
     104             :     }
     105         277 :     return hamiltonian->get_matrix().diagonal().real();
     106             : }
     107             : 
     108             : template <typename Derived>
     109             : const Eigen::SparseMatrix<typename System<Derived>::scalar_t, Eigen::RowMajor> &
     110          58 : System<Derived>::get_matrix() const {
     111          58 :     if (hamiltonian_requires_construction) {
     112          31 :         construct_hamiltonian();
     113          31 :         hamiltonian_requires_construction = false;
     114             :     }
     115          58 :     return hamiltonian->get_matrix();
     116             : }
     117             : 
     118             : template <typename Derived>
     119             : const Transformation<typename System<Derived>::scalar_t> &
     120           0 : System<Derived>::get_transformation() const {
     121           0 :     if (hamiltonian_requires_construction) {
     122           0 :         construct_hamiltonian();
     123           0 :         hamiltonian_requires_construction = false;
     124             :     }
     125           0 :     return hamiltonian->get_transformation();
     126             : }
     127             : 
     128             : template <typename Derived>
     129             : Transformation<typename System<Derived>::scalar_t>
     130           0 : System<Derived>::get_rotator(real_t alpha, real_t beta, real_t gamma) const {
     131           0 :     if (hamiltonian_requires_construction) {
     132           0 :         construct_hamiltonian();
     133           0 :         hamiltonian_requires_construction = false;
     134             :     }
     135           0 :     return hamiltonian->get_rotator(alpha, beta, gamma);
     136             : }
     137             : 
     138             : template <typename Derived>
     139         112 : Sorting System<Derived>::get_sorter(const std::vector<TransformationType> &labels) const {
     140         112 :     if (hamiltonian_requires_construction) {
     141           0 :         construct_hamiltonian();
     142           0 :         hamiltonian_requires_construction = false;
     143             :     }
     144         112 :     return hamiltonian->get_sorter(labels);
     145             : }
     146             : 
     147             : template <typename Derived>
     148             : std::vector<IndicesOfBlock>
     149           0 : System<Derived>::get_indices_of_blocks(const std::vector<TransformationType> &labels) const {
     150           0 :     if (hamiltonian_requires_construction) {
     151           0 :         construct_hamiltonian();
     152           0 :         hamiltonian_requires_construction = false;
     153             :     }
     154           0 :     return hamiltonian->get_indices_of_blocks(labels);
     155             : }
     156             : 
     157             : template <typename Derived>
     158           0 : System<Derived> &System<Derived>::transform(const Transformation<scalar_t> &transformation) {
     159           0 :     if (hamiltonian_requires_construction) {
     160           0 :         construct_hamiltonian();
     161           0 :         hamiltonian_requires_construction = false;
     162             :     }
     163           0 :     hamiltonian = std::make_unique<operator_t>(hamiltonian->transformed(transformation));
     164             : 
     165             :     // A transformed system might have lost its block-diagonalizability if the
     166             :     // transformation was not a sorting
     167           0 :     blockdiagonalizing_labels.clear();
     168             : 
     169           0 :     return *this;
     170             : }
     171             : 
     172             : template <typename Derived>
     173         112 : System<Derived> &System<Derived>::transform(const Sorting &transformation) {
     174         112 :     if (hamiltonian_requires_construction) {
     175           0 :         construct_hamiltonian();
     176           0 :         hamiltonian_requires_construction = false;
     177             :     }
     178         112 :     hamiltonian = std::make_unique<operator_t>(hamiltonian->transformed(transformation));
     179             : 
     180         112 :     return *this;
     181             : }
     182             : 
     183             : template <typename Derived>
     184         112 : System<Derived> &System<Derived>::diagonalize(const DiagonalizerInterface<scalar_t> &diagonalizer,
     185             :                                               std::optional<real_t> min_eigenenergy,
     186             :                                               std::optional<real_t> max_eigenenergy, double rtol) {
     187         112 :     if (hamiltonian_requires_construction) {
     188         105 :         construct_hamiltonian();
     189         105 :         hamiltonian_requires_construction = false;
     190             :     }
     191             : 
     192         112 :     if (hamiltonian_is_diagonal) {
     193           1 :         return *this;
     194             :     }
     195             : 
     196         111 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenvectors;
     197         111 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenenergies;
     198             : 
     199             :     // Sort the Hamiltonian according to the block structure
     200         111 :     if (!blockdiagonalizing_labels.empty()) {
     201          71 :         auto sorter = hamiltonian->get_sorter(blockdiagonalizing_labels);
     202          71 :         hamiltonian = std::make_unique<operator_t>(hamiltonian->transformed(sorter));
     203          71 :     }
     204             : 
     205             :     // Get the indices of the blocks
     206         111 :     auto blocks = hamiltonian->get_indices_of_blocks(blockdiagonalizing_labels);
     207             : 
     208         111 :     assert((blockdiagonalizing_labels.empty() && blocks.size() == 1) ||
     209             :            !blockdiagonalizing_labels.empty());
     210             : 
     211         222 :     SPDLOG_DEBUG("Diagonalizing the Hamiltonian with {} blocks.", blocks.size());
     212             : 
     213             :     // Diagonalize the blocks in parallel
     214         111 :     std::vector<Eigen::VectorX<real_t>> eigenenergies_blocks(blocks.size());
     215         111 :     std::vector<Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>> eigenvectors_blocks(blocks.size());
     216         111 :     oneapi::tbb::parallel_for(
     217         451 :         oneapi::tbb::blocked_range<size_t>(0, blocks.size()), [&](const auto &range) {
     218         680 :             for (size_t idx = range.begin(); idx != range.end(); ++idx) {
     219        1063 :                 auto eigensys = min_eigenenergy.has_value() || max_eigenenergy.has_value()
     220         340 :                     ? diagonalizer.eigh(
     221         129 :                           hamiltonian->get_matrix().block(blocks[idx].start, blocks[idx].start,
     222          86 :                                                           blocks[idx].size(), blocks[idx].size()),
     223             :                           min_eigenenergy, max_eigenenergy, rtol)
     224         297 :                     : diagonalizer.eigh(
     225        1231 :                           hamiltonian->get_matrix().block(blocks[idx].start, blocks[idx].start,
     226         594 :                                                           blocks[idx].size(), blocks[idx].size()),
     227             :                           rtol);
     228         340 :                 eigenvectors_blocks[idx] = eigensys.eigenvectors;
     229         340 :                 eigenenergies_blocks[idx] = eigensys.eigenvalues;
     230             :             }
     231             :         });
     232             : 
     233             :     // Get the number of non-zeros per row of the combined eigenvector matrix
     234         111 :     std::vector<Eigen::Index> non_zeros_per_inner_index;
     235         111 :     non_zeros_per_inner_index.reserve(hamiltonian->get_matrix().rows());
     236         111 :     Eigen::Index num_rows = 0;
     237         111 :     Eigen::Index num_cols = 0;
     238         451 :     for (const auto &matrix : eigenvectors_blocks) {
     239        7194 :         for (int i = 0; i < matrix.outerSize(); ++i) {
     240       13708 :             non_zeros_per_inner_index.push_back(matrix.outerIndexPtr()[i + 1] -
     241        6854 :                                                 matrix.outerIndexPtr()[i]);
     242             :         }
     243         340 :         num_rows += matrix.rows();
     244         340 :         num_cols += matrix.cols();
     245             :     }
     246             : 
     247         111 :     assert(static_cast<size_t>(num_rows) == hamiltonian->get_basis()->get_number_of_kets());
     248         111 :     assert(static_cast<size_t>(num_cols) <= hamiltonian->get_basis()->get_number_of_states());
     249             : 
     250         111 :     eigenvectors.resize(num_rows, num_cols);
     251         111 :     eigenenergies.resize(num_cols, num_cols);
     252             : 
     253         111 :     if (num_cols > 0) {
     254             :         // Get the combined eigenvector matrix (in case of an restricted energy range, it is not
     255             :         // square)
     256         109 :         eigenvectors.reserve(non_zeros_per_inner_index);
     257         109 :         Eigen::Index offset_rows = 0;
     258         109 :         Eigen::Index offset_cols = 0;
     259         447 :         for (const auto &matrix : eigenvectors_blocks) {
     260        7112 :             for (Eigen::Index i = 0; i < matrix.outerSize(); ++i) {
     261        6774 :                 for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     262             :                          matrix, i);
     263      200952 :                      it; ++it) {
     264      194178 :                     eigenvectors.insert(it.row() + offset_rows, it.col() + offset_cols) =
     265      194178 :                         it.value();
     266             :                 }
     267             :             }
     268         338 :             offset_rows += matrix.rows();
     269         338 :             offset_cols += matrix.cols();
     270             :         }
     271         109 :         eigenvectors.makeCompressed();
     272             : 
     273         109 :         assert(
     274             :             eigenvectors.nonZeros() ==
     275             :             std::accumulate(non_zeros_per_inner_index.begin(), non_zeros_per_inner_index.end(), 0));
     276             : 
     277             :         // Get the combined eigenenergy matrix
     278         109 :         eigenenergies.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     279         109 :         Eigen::Index offset = 0;
     280         447 :         for (const auto &matrix : eigenenergies_blocks) {
     281        6671 :             for (int i = 0; i < matrix.size(); ++i) {
     282        6333 :                 eigenenergies.insert(i + offset, i + offset) = matrix(i);
     283             :             }
     284         338 :             offset += matrix.size();
     285             :         }
     286         109 :         eigenenergies.makeCompressed();
     287             : 
     288             :         // Fix phase ambiguity
     289         109 :         std::vector<scalar_t> map_col_to_max(num_cols, 0);
     290        6883 :         for (int row = 0; row < eigenvectors.outerSize(); ++row) {
     291        6774 :             for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     292             :                      eigenvectors, row);
     293      200952 :                  it; ++it) {
     294      194178 :                 if (std::abs(it.value()) > std::abs(map_col_to_max[it.col()])) {
     295       37158 :                     map_col_to_max[it.col()] = it.value();
     296             :                 }
     297             :             }
     298             :         }
     299             : 
     300         109 :         Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> phase_matrix;
     301         109 :         phase_matrix.resize(num_cols, num_cols);
     302         109 :         phase_matrix.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     303        6442 :         for (int i = 0; i < num_cols; ++i) {
     304        6333 :             phase_matrix.insert(i, i) = std::abs(map_col_to_max[i]) / map_col_to_max[i];
     305             :         }
     306         109 :         phase_matrix.makeCompressed();
     307             : 
     308         109 :         eigenvectors = eigenvectors * phase_matrix;
     309         109 :     }
     310             : 
     311             :     // Store the diagonalized hamiltonian
     312         111 :     hamiltonian->get_matrix() = eigenenergies;
     313         111 :     hamiltonian->get_basis() = hamiltonian->get_basis()->transformed(eigenvectors);
     314             : 
     315         111 :     hamiltonian_is_diagonal = true;
     316             : 
     317         111 :     return *this;
     318         111 : }
     319             : 
     320             : template <typename Derived>
     321          51 : bool System<Derived>::is_diagonal() const {
     322          51 :     if (hamiltonian_requires_construction) {
     323          11 :         construct_hamiltonian();
     324          11 :         hamiltonian_requires_construction = false;
     325             :     }
     326          51 :     return hamiltonian_is_diagonal;
     327             : }
     328             : 
     329             : // Explicit instantiation
     330             : template class System<SystemAtom<double>>;
     331             : template class System<SystemAtom<std::complex<double>>>;
     332             : template class System<SystemPair<double>>;
     333             : template class System<SystemPair<std::complex<double>>>;
     334             : } // namespace pairinteraction