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/system/SystemAtom.hpp"
      11             : #include "pairinteraction/system/SystemPair.hpp"
      12             : #include "pairinteraction/utils/TaskControl.hpp"
      13             : #include "pairinteraction/utils/eigen_assertion.hpp"
      14             : #include "pairinteraction/utils/eigen_compat.hpp"
      15             : 
      16             : #include <Eigen/SparseCore>
      17             : #include <algorithm>
      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        1228 : System<Derived>::System(std::shared_ptr<const basis_t> basis)
      29        1228 :     : basis(std::move(basis)),
      30        1228 :       matrix(static_cast<Eigen::Index>(this->basis->get_number_of_states()),
      31        3684 :              static_cast<Eigen::Index>(this->basis->get_number_of_states())) {}
      32             : 
      33             : template <typename Derived>
      34        1531 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_basis() const {
      35        1531 :     if (hamiltonian_requires_construction) {
      36          83 :         construct_hamiltonian();
      37          83 :         hamiltonian_requires_construction = false;
      38             :     }
      39        1531 :     return basis;
      40             : }
      41             : 
      42             : template <typename Derived>
      43         437 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_eigenbasis() const {
      44         437 :     if (hamiltonian_requires_construction) {
      45           2 :         construct_hamiltonian();
      46           2 :         hamiltonian_requires_construction = false;
      47             :     }
      48         437 :     if (!this->is_diagonal()) {
      49           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
      50             :     }
      51         437 :     return basis;
      52             : }
      53             : 
      54             : template <typename Derived>
      55        1898 : Eigen::VectorX<typename System<Derived>::real_t> System<Derived>::get_eigenenergies() const {
      56        1898 :     if (hamiltonian_requires_construction) {
      57           0 :         construct_hamiltonian();
      58           0 :         hamiltonian_requires_construction = false;
      59             :     }
      60        1898 :     if (!this->is_diagonal()) {
      61           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
      62             :     }
      63        2888 :     return matrix.diagonal().real();
      64             : }
      65             : 
      66             : template <typename Derived>
      67             : const Eigen::SparseMatrix<typename System<Derived>::scalar_t, Eigen::RowMajor> &
      68         156 : System<Derived>::get_matrix() const {
      69         156 :     if (hamiltonian_requires_construction) {
      70         115 :         construct_hamiltonian();
      71         115 :         hamiltonian_requires_construction = false;
      72             :     }
      73         156 :     return matrix;
      74             : }
      75             : 
      76             : template <typename Derived>
      77             : const Transformation<typename System<Derived>::scalar_t> &
      78           0 : System<Derived>::get_transformation() const {
      79           0 :     if (hamiltonian_requires_construction) {
      80           0 :         construct_hamiltonian();
      81           0 :         hamiltonian_requires_construction = false;
      82             :     }
      83           0 :     return basis->get_transformation();
      84             : }
      85             : 
      86             : template <typename Derived>
      87        1213 : Sorting System<Derived>::get_sorter(const std::vector<TransformationType> &labels) const {
      88        1213 :     if (hamiltonian_requires_construction) {
      89           0 :         construct_hamiltonian();
      90           0 :         hamiltonian_requires_construction = false;
      91             :     }
      92             : 
      93        1213 :     basis->perform_sorter_checks(labels);
      94             : 
      95        1198 :     auto it = std::find(labels.begin(), labels.end(), TransformationType::SORT_BY_ENERGY);
      96        1182 :     std::vector<TransformationType> before_energy(labels.begin(), it);
      97        1195 :     bool contains_energy = (it != labels.end());
      98        1193 :     std::vector<TransformationType> after_energy(contains_energy ? it + 1 : labels.end(),
      99             :                                                  labels.end());
     100             : 
     101        1193 :     Sorting transformation;
     102        1195 :     transformation.matrix.resize(matrix.rows());
     103        1193 :     transformation.matrix.setIdentity();
     104             : 
     105        1206 :     if (!before_energy.empty()) {
     106         578 :         basis->get_sorter_without_checks(before_energy, transformation);
     107             :     }
     108             : 
     109        1203 :     if (contains_energy) {
     110         626 :         std::vector<real_t> energies_of_states;
     111         626 :         energies_of_states.reserve(matrix.rows());
     112       30209 :         for (int i = 0; i < matrix.rows(); ++i) {
     113       29625 :             energies_of_states.push_back(std::real(matrix.coeff(i, i)));
     114             :         }
     115             : 
     116         629 :         std::stable_sort(
     117         629 :             transformation.matrix.indices().data(),
     118         607 :             transformation.matrix.indices().data() + transformation.matrix.indices().size(),
     119      126861 :             [&](int i, int j) { return energies_of_states[i] < energies_of_states[j]; });
     120             : 
     121         626 :         transformation.transformation_type.push_back(TransformationType::SORT_BY_ENERGY);
     122         626 :     }
     123             : 
     124        1202 :     if (!after_energy.empty()) {
     125           0 :         basis->get_sorter_without_checks(after_energy, transformation);
     126             :     }
     127             : 
     128        1193 :     if (labels != transformation.transformation_type) {
     129           0 :         throw std::invalid_argument("The states could not be sorted by all the requested labels.");
     130             :     }
     131             : 
     132        2380 :     return transformation;
     133        1194 : }
     134             : 
     135             : template <typename Derived>
     136             : std::vector<IndicesOfBlock>
     137         628 : System<Derived>::get_indices_of_blocks(const std::vector<TransformationType> &labels) const {
     138         628 :     if (hamiltonian_requires_construction) {
     139           0 :         construct_hamiltonian();
     140           0 :         hamiltonian_requires_construction = false;
     141             :     }
     142             : 
     143         628 :     basis->perform_sorter_checks(labels);
     144             : 
     145         621 :     std::set<TransformationType> unique_labels(labels.begin(), labels.end());
     146         626 :     basis->perform_blocks_checks(unique_labels);
     147             : 
     148         628 :     auto it = unique_labels.find(TransformationType::SORT_BY_ENERGY);
     149         619 :     bool contains_energy = (it != unique_labels.end());
     150         620 :     if (contains_energy) {
     151           0 :         unique_labels.erase(it);
     152             :     }
     153             : 
     154         620 :     IndicesOfBlocksCreator blocks_creator({0, static_cast<size_t>(matrix.rows())});
     155             : 
     156         624 :     if (!unique_labels.empty()) {
     157         562 :         basis->get_indices_of_blocks_without_checks(unique_labels, blocks_creator);
     158             :     }
     159             : 
     160         630 :     if (contains_energy && matrix.rows() > 0) {
     161           0 :         scalar_t last_energy = std::real(matrix.coeff(0, 0));
     162           0 :         for (int i = 0; i < matrix.rows(); ++i) {
     163           0 :             if (std::real(matrix.coeff(i, i)) != last_energy) {
     164           0 :                 blocks_creator.add(i);
     165           0 :                 last_energy = std::real(matrix.coeff(i, i));
     166             :             }
     167             :         }
     168             :     }
     169             : 
     170        1254 :     return blocks_creator.create();
     171         627 : }
     172             : 
     173             : template <typename Derived>
     174           0 : System<Derived> &System<Derived>::transform(const Transformation<scalar_t> &transformation) {
     175           0 :     if (hamiltonian_requires_construction) {
     176           0 :         construct_hamiltonian();
     177           0 :         hamiltonian_requires_construction = false;
     178             :     }
     179             : 
     180           0 :     set_task_status("Applying transformation...");
     181           0 :     if (matrix.cols() != 0) {
     182           0 :         matrix = transformation.matrix.adjoint() * matrix * transformation.matrix;
     183           0 :         basis = basis->transformed(transformation);
     184             :     }
     185             : 
     186           0 :     hamiltonian_is_diagonal = false;
     187             : 
     188             :     // A transformed system might have lost its block-diagonalizability if the
     189             :     // transformation was not a sorting
     190           0 :     blockdiagonalizing_labels.clear();
     191             : 
     192           0 :     return *this;
     193             : }
     194             : 
     195             : template <typename Derived>
     196         625 : System<Derived> &System<Derived>::transform(const Sorting &transformation) {
     197         625 :     if (hamiltonian_requires_construction) {
     198           0 :         construct_hamiltonian();
     199           0 :         hamiltonian_requires_construction = false;
     200             :     }
     201             : 
     202         625 :     set_task_status("Applying transformation...");
     203         633 :     if (matrix.cols() != 0) {
     204         631 :         matrix = matrix.twistedBy(transformation.matrix.inverse());
     205         626 :         basis = basis->transformed(transformation);
     206             :     }
     207             : 
     208         630 :     return *this;
     209             : }
     210             : 
     211             : template <typename Derived>
     212         655 : System<Derived> &System<Derived>::diagonalize(const DiagonalizerInterface<scalar_t> &diagonalizer,
     213             :                                               std::optional<real_t> min_eigenenergy,
     214             :                                               std::optional<real_t> max_eigenenergy, double rtol,
     215             :                                               bool sort_by_energy) {
     216         655 :     set_task_status("Preparing Hamiltonian...");
     217             : 
     218         658 :     if (hamiltonian_requires_construction) {
     219         654 :         construct_hamiltonian();
     220         653 :         hamiltonian_requires_construction = false;
     221             :     }
     222             : 
     223         657 :     if (this->is_diagonal()) {
     224          27 :         if (sort_by_energy && !this->is_diagonal_and_sorted_by_energy()) {
     225           2 :             transform(get_sorter({TransformationType::SORT_BY_ENERGY}));
     226             :         }
     227          27 :         return *this;
     228             :     }
     229             : 
     230         619 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenvectors;
     231         625 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenenergies;
     232             : 
     233             :     // Sort the Hamiltonian according to the block structure
     234         629 :     if (!blockdiagonalizing_labels.empty()) {
     235         577 :         auto sorter = get_sorter(blockdiagonalizing_labels);
     236         576 :         matrix = matrix.twistedBy(sorter.matrix.inverse());
     237         579 :         basis = basis->transformed(sorter);
     238         574 :     }
     239             : 
     240             :     // Get the indices of the blocks
     241         628 :     auto blocks = get_indices_of_blocks(blockdiagonalizing_labels);
     242             : 
     243         627 :     assert((blockdiagonalizing_labels.empty() && blocks.size() == 1) ||
     244             :            !blockdiagonalizing_labels.empty());
     245             : 
     246         624 :     SPDLOG_DEBUG("Diagonalizing the Hamiltonian with {} blocks.", blocks.size());
     247             : 
     248             :     // Diagonalize the blocks in parallel
     249         631 :     std::vector<Eigen::VectorX<real_t>> eigenenergies_blocks(blocks.size());
     250         631 :     std::vector<Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>> eigenvectors_blocks(blocks.size());
     251         631 :     oneapi::tbb::parallel_for(
     252        1942 :         oneapi::tbb::blocked_range<size_t>(0, blocks.size()), [&](const auto &range) {
     253        2645 :             for (size_t idx = range.begin(); idx != range.end(); ++idx) {
     254        1326 :                 set_task_status("Diagonalizing Hamiltonian blocks...");
     255        4386 :                 auto eigensys = min_eigenenergy.has_value() || max_eigenenergy.has_value()
     256        2072 :                     ? diagonalizer.eigh(matrix.block(blocks[idx].start, blocks[idx].start,
     257         730 :                                                      blocks[idx].size(), blocks[idx].size()),
     258             :                                         min_eigenenergy, max_eigenenergy, rtol)
     259        3288 :                     : diagonalizer.eigh(matrix.block(blocks[idx].start, blocks[idx].start,
     260        1954 :                                                      blocks[idx].size(), blocks[idx].size()),
     261             :                                         rtol);
     262        1332 :                 eigenvectors_blocks[idx] = eigensys.eigenvectors;
     263        1331 :                 eigenenergies_blocks[idx] = eigensys.eigenvalues;
     264             :             }
     265             :         });
     266             : 
     267             :     // Get the number of non-zeros per row of the combined eigenvector matrix
     268         628 :     std::vector<Eigen::Index> non_zeros_per_inner_index;
     269         628 :     non_zeros_per_inner_index.reserve(matrix.rows());
     270         631 :     Eigen::Index num_rows = 0;
     271         631 :     Eigen::Index num_cols = 0;
     272        1969 :     for (const auto &matrix : eigenvectors_blocks) {
     273        1333 :         set_task_status("Collecting eigenvectors...");
     274       36247 :         for (int i = 0; i < matrix.outerSize(); ++i) {
     275       69860 :             non_zeros_per_inner_index.push_back(matrix.outerIndexPtr()[i + 1] -
     276       34927 :                                                 matrix.outerIndexPtr()[i]);
     277             :         }
     278        1312 :         num_rows += matrix.rows();
     279        1335 :         num_cols += matrix.cols();
     280             :     }
     281             : 
     282         630 :     assert(static_cast<size_t>(num_rows) == basis->get_number_of_states());
     283         629 :     assert(static_cast<size_t>(num_cols) <= basis->get_number_of_states());
     284             : 
     285         625 :     eigenvectors.resize(num_rows, num_cols);
     286         631 :     eigenenergies.resize(num_cols, num_cols);
     287             : 
     288         631 :     if (num_cols > 0) {
     289             :         // Get the combined eigenvector matrix (in case of an restricted energy range, it is not
     290             :         // square)
     291         629 :         eigenvectors.reserve(non_zeros_per_inner_index);
     292         627 :         Eigen::Index offset_rows = 0;
     293         627 :         Eigen::Index offset_cols = 0;
     294        1964 :         for (const auto &matrix : eigenvectors_blocks) {
     295        1334 :             set_task_status("Combining eigenvectors...");
     296       36370 :             for (Eigen::Index i = 0; i < matrix.outerSize(); ++i) {
     297       35038 :                 for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     298             :                          matrix, i);
     299      639748 :                      it; ++it) {
     300      605695 :                     eigenvectors.insert(it.row() + offset_rows, it.col() + offset_cols) =
     301      605152 :                         it.value();
     302             :                 }
     303             :             }
     304        1331 :             offset_rows += matrix.rows();
     305        1335 :             offset_cols += matrix.cols();
     306             :         }
     307         626 :         eigenvectors.makeCompressed();
     308             : 
     309         629 :         assert(
     310             :             eigenvectors.nonZeros() ==
     311             :             std::accumulate(non_zeros_per_inner_index.begin(), non_zeros_per_inner_index.end(), 0));
     312             : 
     313             :         // Get the combined eigenenergy matrix
     314         628 :         eigenenergies.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     315         627 :         Eigen::Index offset = 0;
     316        1966 :         for (const auto &matrix : eigenenergies_blocks) {
     317        1339 :             set_task_status("Combining eigenenergies...");
     318       31207 :             for (int i = 0; i < matrix.size(); ++i) {
     319       29862 :                 eigenenergies.insert(i + offset, i + offset) = matrix(i);
     320             :             }
     321        1315 :             offset += matrix.size();
     322             :         }
     323         629 :         eigenenergies.makeCompressed();
     324             : 
     325             :         // Fix phase ambiguity
     326         628 :         std::vector<scalar_t> map_col_to_max(num_cols, 0);
     327       35656 :         for (int row = 0; row < eigenvectors.outerSize(); ++row) {
     328       35024 :             set_task_status("Normalizing eigenvector phases...");
     329       35079 :             for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     330             :                      eigenvectors, row);
     331      636858 :                  it; ++it) {
     332      601074 :                 if (std::abs(it.value()) > std::abs(map_col_to_max[it.col()])) {
     333      139594 :                     map_col_to_max[it.col()] = it.value();
     334             :                 }
     335             :             }
     336             :         }
     337             : 
     338         626 :         Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> phase_matrix;
     339         629 :         phase_matrix.resize(num_cols, num_cols);
     340         629 :         phase_matrix.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     341       30464 :         for (int i = 0; i < num_cols; ++i) {
     342       29865 :             phase_matrix.insert(i, i) = std::abs(map_col_to_max[i]) / map_col_to_max[i];
     343             :         }
     344         599 :         phase_matrix.makeCompressed();
     345             : 
     346         629 :         set_task_status("Applying eigenvector phases...");
     347         629 :         eigenvectors = eigenvectors * phase_matrix;
     348         629 :     }
     349             : 
     350             :     // Store the diagonalized hamiltonian
     351         629 :     matrix = eigenenergies;
     352         628 :     basis = basis->transformed(eigenvectors);
     353             : 
     354         629 :     hamiltonian_is_diagonal = true;
     355         629 :     if (sort_by_energy) {
     356         629 :         transform(get_sorter({TransformationType::SORT_BY_ENERGY}));
     357             :     }
     358             : 
     359         627 :     return *this;
     360         627 : }
     361             : 
     362             : template <typename Derived>
     363        5078 : bool System<Derived>::is_diagonal() const {
     364        5078 :     if (hamiltonian_requires_construction) {
     365         373 :         construct_hamiltonian();
     366         373 :         hamiltonian_requires_construction = false;
     367             :     }
     368             : 
     369        5078 :     if (!hamiltonian_is_diagonal) {
     370        1102 :         real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
     371             : 
     372       29894 :         for (int row = 0; row < matrix.outerSize(); ++row) {
     373       29401 :             for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(matrix,
     374             :                                                                                            row);
     375       59393 :                  it; ++it) {
     376       30512 :                 if (it.row() != it.col() && std::abs(it.value()) > numerical_precision) {
     377         622 :                     return false;
     378             :                 }
     379             :             }
     380             :         }
     381             : 
     382         468 :         hamiltonian_is_diagonal = true;
     383             :     }
     384             : 
     385        4444 :     return true;
     386             : }
     387             : 
     388             : template <typename Derived>
     389        2079 : bool System<Derived>::is_diagonal_and_sorted_by_energy() const {
     390        2079 :     if (!this->is_diagonal()) {
     391           0 :         return false;
     392             :     }
     393             : 
     394        2079 :     Eigen::VectorX<real_t> eigenenergies = matrix.diagonal().real();
     395        2079 :     return std::is_sorted(eigenenergies.data(), eigenenergies.data() + eigenenergies.size());
     396        2079 : }
     397             : 
     398             : // Explicit instantiation
     399             : template class System<SystemAtom<double>>;
     400             : template class System<SystemAtom<std::complex<double>>>;
     401             : template class System<SystemPair<double>>;
     402             : template class System<SystemPair<std::complex<double>>>;
     403             : } // namespace pairinteraction