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          69 : System<Derived>::System(std::shared_ptr<const basis_t> basis)
      29          69 :     : hamiltonian(std::make_unique<typename System<Derived>::operator_t>(std::move(basis))) {}
      30             : 
      31             : template <typename Derived>
      32          28 : System<Derived>::System(const System &other)
      33          28 :     : hamiltonian(std::make_unique<typename System<Derived>::operator_t>(*other.hamiltonian)),
      34          28 :       hamiltonian_requires_construction(other.hamiltonian_requires_construction),
      35          28 :       hamiltonian_is_diagonal(other.hamiltonian_is_diagonal),
      36          56 :       blockdiagonalizing_labels(other.blockdiagonalizing_labels) {}
      37             : 
      38             : template <typename Derived>
      39          16 : System<Derived>::System(System &&other) noexcept
      40          16 :     : hamiltonian(std::move(other.hamiltonian)),
      41          16 :       hamiltonian_requires_construction(other.hamiltonian_requires_construction),
      42          16 :       hamiltonian_is_diagonal(other.hamiltonian_is_diagonal),
      43          32 :       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         113 : 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          26 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_basis() const {
      77          26 :     if (hamiltonian_requires_construction) {
      78           0 :         construct_hamiltonian();
      79           0 :         hamiltonian_requires_construction = false;
      80             :     }
      81          26 :     return hamiltonian->get_basis();
      82             : }
      83             : 
      84             : template <typename Derived>
      85          55 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_eigenbasis() const {
      86          55 :     if (hamiltonian_requires_construction) {
      87           0 :         construct_hamiltonian();
      88           0 :         hamiltonian_requires_construction = false;
      89             :     }
      90          55 :     if (!hamiltonian_is_diagonal) {
      91           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
      92             :     }
      93          55 :     return hamiltonian->get_basis();
      94             : }
      95             : 
      96             : template <typename Derived>
      97          71 : Eigen::VectorX<typename System<Derived>::real_t> System<Derived>::get_eigenenergies() const {
      98          71 :     if (hamiltonian_requires_construction) {
      99           0 :         construct_hamiltonian();
     100           0 :         hamiltonian_requires_construction = false;
     101             :     }
     102          71 :     if (!hamiltonian_is_diagonal) {
     103           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
     104             :     }
     105         122 :     return hamiltonian->get_matrix().diagonal().real();
     106             : }
     107             : 
     108             : template <typename Derived>
     109             : const Eigen::SparseMatrix<typename System<Derived>::scalar_t, Eigen::RowMajor> &
     110           9 : System<Derived>::get_matrix() const {
     111           9 :     if (hamiltonian_requires_construction) {
     112           4 :         construct_hamiltonian();
     113           4 :         hamiltonian_requires_construction = false;
     114             :     }
     115           9 :     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          42 : Sorting System<Derived>::get_sorter(const std::vector<TransformationType> &labels) const {
     140          42 :     if (hamiltonian_requires_construction) {
     141           0 :         construct_hamiltonian();
     142           0 :         hamiltonian_requires_construction = false;
     143             :     }
     144          42 :     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          42 : System<Derived> &System<Derived>::transform(const Sorting &transformation) {
     174          42 :     if (hamiltonian_requires_construction) {
     175           0 :         construct_hamiltonian();
     176           0 :         hamiltonian_requires_construction = false;
     177             :     }
     178          42 :     hamiltonian = std::make_unique<operator_t>(hamiltonian->transformed(transformation));
     179             : 
     180          42 :     return *this;
     181             : }
     182             : 
     183             : template <typename Derived>
     184          57 : 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          57 :     if (hamiltonian_requires_construction) {
     188          56 :         construct_hamiltonian();
     189          56 :         hamiltonian_requires_construction = false;
     190             :     }
     191             : 
     192          57 :     if (hamiltonian_is_diagonal) {
     193           1 :         return *this;
     194             :     }
     195             : 
     196          56 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenvectors;
     197          56 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenenergies;
     198             : 
     199             :     // Sort the Hamiltonian according to the block structure
     200          56 :     if (!blockdiagonalizing_labels.empty()) {
     201          30 :         auto sorter = hamiltonian->get_sorter(blockdiagonalizing_labels);
     202          30 :         hamiltonian = std::make_unique<operator_t>(hamiltonian->transformed(sorter));
     203          30 :     }
     204             : 
     205             :     // Get the indices of the blocks
     206          56 :     auto blocks = hamiltonian->get_indices_of_blocks(blockdiagonalizing_labels);
     207             : 
     208          56 :     assert((blockdiagonalizing_labels.empty() && blocks.size() == 1) ||
     209             :            !blockdiagonalizing_labels.empty());
     210             : 
     211         112 :     SPDLOG_DEBUG("Diagonalizing the Hamiltonian with {} blocks.", blocks.size());
     212             : 
     213             :     // Diagonalize the blocks in parallel
     214          56 :     std::vector<Eigen::VectorX<real_t>> eigenenergies_blocks(blocks.size());
     215          56 :     std::vector<Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>> eigenvectors_blocks(blocks.size());
     216          56 :     oneapi::tbb::parallel_for(
     217         202 :         oneapi::tbb::blocked_range<size_t>(0, blocks.size()), [&](const auto &range) {
     218         292 :             for (size_t idx = range.begin(); idx != range.end(); ++idx) {
     219         470 :                 auto eigensys = min_eigenenergy.has_value() || max_eigenenergy.has_value()
     220         145 :                     ? diagonalizer.eigh(
     221          98 :                           hamiltonian->get_matrix().block(blocks[idx].start, blocks[idx].start,
     222          66 :                                                           blocks[idx].size(), blocks[idx].size()),
     223             :                           min_eigenenergy, max_eigenenergy, rtol)
     224         113 :                     : diagonalizer.eigh(
     225         485 :                           hamiltonian->get_matrix().block(blocks[idx].start, blocks[idx].start,
     226         224 :                                                           blocks[idx].size(), blocks[idx].size()),
     227             :                           rtol);
     228         146 :                 eigenvectors_blocks[idx] = eigensys.eigenvectors;
     229         147 :                 eigenenergies_blocks[idx] = eigensys.eigenvalues;
     230             :             }
     231             :         });
     232             : 
     233             :     // Get the number of non-zeros per row of the combined eigenvector matrix
     234          56 :     std::vector<Eigen::Index> non_zeros_per_inner_index;
     235          56 :     non_zeros_per_inner_index.reserve(hamiltonian->get_matrix().rows());
     236          56 :     Eigen::Index num_rows = 0;
     237          56 :     Eigen::Index num_cols = 0;
     238         203 :     for (const auto &matrix : eigenvectors_blocks) {
     239        2485 :         for (int i = 0; i < matrix.outerSize(); ++i) {
     240        4683 :             non_zeros_per_inner_index.push_back(matrix.outerIndexPtr()[i + 1] -
     241        2343 :                                                 matrix.outerIndexPtr()[i]);
     242             :         }
     243         147 :         num_rows += matrix.rows();
     244         147 :         num_cols += matrix.cols();
     245             :     }
     246             : 
     247          56 :     assert(static_cast<size_t>(num_rows) == hamiltonian->get_basis()->get_number_of_kets());
     248          56 :     assert(static_cast<size_t>(num_cols) <= hamiltonian->get_basis()->get_number_of_states());
     249             : 
     250          56 :     eigenvectors.resize(num_rows, num_cols);
     251          55 :     eigenenergies.resize(num_cols, num_cols);
     252             : 
     253          56 :     if (num_cols > 0) {
     254             :         // Get the combined eigenvector matrix (in case of an restricted energy range, it is not
     255             :         // square)
     256          54 :         eigenvectors.reserve(non_zeros_per_inner_index);
     257          53 :         Eigen::Index offset_rows = 0;
     258          53 :         Eigen::Index offset_cols = 0;
     259         198 :         for (const auto &matrix : eigenvectors_blocks) {
     260        2423 :             for (Eigen::Index i = 0; i < matrix.outerSize(); ++i) {
     261        2277 :                 for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     262             :                          matrix, i);
     263       43599 :                      it; ++it) {
     264       41336 :                     eigenvectors.insert(it.row() + offset_rows, it.col() + offset_cols) =
     265       41300 :                         it.value();
     266             :                 }
     267             :             }
     268         144 :             offset_rows += matrix.rows();
     269         145 :             offset_cols += matrix.cols();
     270             :         }
     271          54 :         eigenvectors.makeCompressed();
     272             : 
     273          54 :         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          54 :         eigenenergies.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     279          54 :         Eigen::Index offset = 0;
     280         199 :         for (const auto &matrix : eigenenergies_blocks) {
     281        2225 :             for (int i = 0; i < matrix.size(); ++i) {
     282        2079 :                 eigenenergies.insert(i + offset, i + offset) = matrix(i);
     283             :             }
     284         145 :             offset += matrix.size();
     285             :         }
     286          54 :         eigenenergies.makeCompressed();
     287             : 
     288             :         // Fix phase ambiguity
     289          54 :         std::vector<scalar_t> map_col_to_max(num_cols, 0);
     290        2334 :         for (int row = 0; row < eigenvectors.outerSize(); ++row) {
     291        2280 :             for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     292             :                      eigenvectors, row);
     293       43585 :                  it; ++it) {
     294       41258 :                 if (std::abs(it.value()) > std::abs(map_col_to_max[it.col()])) {
     295       11621 :                     map_col_to_max[it.col()] = it.value();
     296             :                 }
     297             :             }
     298             :         }
     299             : 
     300          54 :         Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> phase_matrix;
     301          54 :         phase_matrix.resize(num_cols, num_cols);
     302          54 :         phase_matrix.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     303        2134 :         for (int i = 0; i < num_cols; ++i) {
     304        2080 :             phase_matrix.insert(i, i) = std::abs(map_col_to_max[i]) / map_col_to_max[i];
     305             :         }
     306          54 :         phase_matrix.makeCompressed();
     307             : 
     308          54 :         eigenvectors = eigenvectors * phase_matrix;
     309          54 :     }
     310             : 
     311             :     // Store the diagonalized hamiltonian
     312          56 :     hamiltonian->get_matrix() = eigenenergies;
     313          56 :     hamiltonian->get_basis() = hamiltonian->get_basis()->transformed(eigenvectors);
     314             : 
     315          56 :     hamiltonian_is_diagonal = true;
     316             : 
     317          56 :     return *this;
     318          56 : }
     319             : 
     320             : template <typename Derived>
     321          26 : bool System<Derived>::is_diagonal() const {
     322          26 :     if (hamiltonian_requires_construction) {
     323           6 :         construct_hamiltonian();
     324           6 :         hamiltonian_requires_construction = false;
     325             :     }
     326          26 :     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