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        1230 : System<Derived>::System(std::shared_ptr<const basis_t> basis)
      29        1230 :     : basis(std::move(basis)),
      30        1230 :       matrix(static_cast<Eigen::Index>(this->basis->get_number_of_states()),
      31        3690 :              static_cast<Eigen::Index>(this->basis->get_number_of_states())) {}
      32             : 
      33             : template <typename Derived>
      34        1471 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_basis() const {
      35        1471 :     if (hamiltonian_requires_construction) {
      36          75 :         construct_hamiltonian();
      37          75 :         hamiltonian_requires_construction = false;
      38             :     }
      39        1471 :     return basis;
      40             : }
      41             : 
      42             : template <typename Derived>
      43         426 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_eigenbasis() const {
      44         426 :     if (hamiltonian_requires_construction) {
      45           2 :         construct_hamiltonian();
      46           2 :         hamiltonian_requires_construction = false;
      47             :     }
      48         426 :     if (!this->is_diagonal()) {
      49           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
      50             :     }
      51         426 :     return basis;
      52             : }
      53             : 
      54             : template <typename Derived>
      55        1874 : Eigen::VectorX<typename System<Derived>::real_t> System<Derived>::get_eigenenergies() const {
      56        1874 :     if (hamiltonian_requires_construction) {
      57          36 :         construct_hamiltonian();
      58          36 :         hamiltonian_requires_construction = false;
      59             :     }
      60        1874 :     if (!this->is_diagonal()) {
      61           0 :         throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
      62             :     }
      63        2942 :     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             : Transformation<typename System<Derived>::scalar_t>
      88           0 : System<Derived>::get_rotator(real_t alpha, real_t beta, real_t gamma) const {
      89           0 :     if (hamiltonian_requires_construction) {
      90           0 :         construct_hamiltonian();
      91           0 :         hamiltonian_requires_construction = false;
      92             :     }
      93           0 :     return basis->get_rotator(alpha, beta, gamma);
      94             : }
      95             : 
      96             : template <typename Derived>
      97        2130 : Sorting System<Derived>::get_sorter(const std::vector<TransformationType> &labels) const {
      98        2130 :     if (hamiltonian_requires_construction) {
      99           0 :         construct_hamiltonian();
     100           0 :         hamiltonian_requires_construction = false;
     101             :     }
     102             : 
     103        2130 :     basis->perform_sorter_checks(labels);
     104             : 
     105        2130 :     auto it = std::find(labels.begin(), labels.end(), TransformationType::SORT_BY_ENERGY);
     106        2130 :     std::vector<TransformationType> before_energy(labels.begin(), it);
     107        2130 :     bool contains_energy = (it != labels.end());
     108        2130 :     std::vector<TransformationType> after_energy(contains_energy ? it + 1 : labels.end(),
     109             :                                                  labels.end());
     110             : 
     111        2130 :     Sorting transformation;
     112        2130 :     transformation.matrix.resize(matrix.rows());
     113        2130 :     transformation.matrix.setIdentity();
     114             : 
     115        2130 :     if (!before_energy.empty()) {
     116         578 :         basis->get_sorter_without_checks(before_energy, transformation);
     117             :     }
     118             : 
     119        2130 :     if (contains_energy) {
     120        1552 :         std::vector<real_t> energies_of_states;
     121        1552 :         energies_of_states.reserve(matrix.rows());
     122      146992 :         for (int i = 0; i < matrix.rows(); ++i) {
     123      145440 :             energies_of_states.push_back(std::real(matrix.coeff(i, i)));
     124             :         }
     125             : 
     126        1552 :         std::stable_sort(
     127        1552 :             transformation.matrix.indices().data(),
     128        1552 :             transformation.matrix.indices().data() + transformation.matrix.indices().size(),
     129      664514 :             [&](int i, int j) { return energies_of_states[i] < energies_of_states[j]; });
     130             : 
     131        1552 :         transformation.transformation_type.push_back(TransformationType::SORT_BY_ENERGY);
     132        1552 :     }
     133             : 
     134        2130 :     if (!after_energy.empty()) {
     135           0 :         basis->get_sorter_without_checks(after_energy, transformation);
     136             :     }
     137             : 
     138        2130 :     if (labels != transformation.transformation_type) {
     139           0 :         throw std::invalid_argument("The states could not be sorted by all the requested labels.");
     140             :     }
     141             : 
     142        4260 :     return transformation;
     143        2130 : }
     144             : 
     145             : template <typename Derived>
     146             : std::vector<IndicesOfBlock>
     147         629 : System<Derived>::get_indices_of_blocks(const std::vector<TransformationType> &labels) const {
     148         629 :     if (hamiltonian_requires_construction) {
     149           0 :         construct_hamiltonian();
     150           0 :         hamiltonian_requires_construction = false;
     151             :     }
     152             : 
     153         629 :     basis->perform_sorter_checks(labels);
     154             : 
     155         629 :     std::set<TransformationType> unique_labels(labels.begin(), labels.end());
     156         629 :     basis->perform_blocks_checks(unique_labels);
     157             : 
     158         629 :     auto it = unique_labels.find(TransformationType::SORT_BY_ENERGY);
     159         629 :     bool contains_energy = (it != unique_labels.end());
     160         629 :     if (contains_energy) {
     161           0 :         unique_labels.erase(it);
     162             :     }
     163             : 
     164         629 :     IndicesOfBlocksCreator blocks_creator({0, static_cast<size_t>(matrix.rows())});
     165             : 
     166         629 :     if (!unique_labels.empty()) {
     167         578 :         basis->get_indices_of_blocks_without_checks(unique_labels, blocks_creator);
     168             :     }
     169             : 
     170         629 :     if (contains_energy && matrix.rows() > 0) {
     171           0 :         scalar_t last_energy = std::real(matrix.coeff(0, 0));
     172           0 :         for (int i = 0; i < matrix.rows(); ++i) {
     173           0 :             if (std::real(matrix.coeff(i, i)) != last_energy) {
     174           0 :                 blocks_creator.add(i);
     175           0 :                 last_energy = std::real(matrix.coeff(i, i));
     176             :             }
     177             :         }
     178             :     }
     179             : 
     180        1258 :     return blocks_creator.create();
     181         629 : }
     182             : 
     183             : template <typename Derived>
     184           0 : System<Derived> &System<Derived>::transform(const Transformation<scalar_t> &transformation) {
     185           0 :     if (hamiltonian_requires_construction) {
     186           0 :         construct_hamiltonian();
     187           0 :         hamiltonian_requires_construction = false;
     188             :     }
     189             : 
     190           0 :     set_task_status("Applying transformation...");
     191           0 :     if (matrix.cols() != 0) {
     192           0 :         matrix = transformation.matrix.adjoint() * matrix * transformation.matrix;
     193           0 :         basis = basis->transformed(transformation);
     194             :     }
     195             : 
     196           0 :     hamiltonian_is_diagonal = false;
     197             : 
     198             :     // A transformed system might have lost its block-diagonalizability if the
     199             :     // transformation was not a sorting
     200           0 :     blockdiagonalizing_labels.clear();
     201             : 
     202           0 :     return *this;
     203             : }
     204             : 
     205             : template <typename Derived>
     206        1552 : System<Derived> &System<Derived>::transform(const Sorting &transformation) {
     207        1552 :     if (hamiltonian_requires_construction) {
     208           0 :         construct_hamiltonian();
     209           0 :         hamiltonian_requires_construction = false;
     210             :     }
     211             : 
     212        1552 :     set_task_status("Applying transformation...");
     213        1552 :     if (matrix.cols() != 0) {
     214        1552 :         matrix = matrix.twistedBy(transformation.matrix.inverse());
     215        1552 :         basis = basis->transformed(transformation);
     216             :     }
     217             : 
     218        1552 :     return *this;
     219             : }
     220             : 
     221             : template <typename Derived>
     222         656 : System<Derived> &System<Derived>::diagonalize(const DiagonalizerInterface<scalar_t> &diagonalizer,
     223             :                                               std::optional<real_t> min_eigenenergy,
     224             :                                               std::optional<real_t> max_eigenenergy, double rtol) {
     225         656 :     set_task_status("Preparing Hamiltonian...");
     226             : 
     227         656 :     if (hamiltonian_requires_construction) {
     228         642 :         construct_hamiltonian();
     229         642 :         hamiltonian_requires_construction = false;
     230             :     }
     231             : 
     232         656 :     if (this->is_diagonal()) {
     233          27 :         return *this;
     234             :     }
     235             : 
     236         629 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenvectors;
     237         629 :     Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenenergies;
     238             : 
     239             :     // Sort the Hamiltonian according to the block structure
     240         629 :     if (!blockdiagonalizing_labels.empty()) {
     241         578 :         auto sorter = get_sorter(blockdiagonalizing_labels);
     242         578 :         matrix = matrix.twistedBy(sorter.matrix.inverse());
     243         578 :         basis = basis->transformed(sorter);
     244         578 :     }
     245             : 
     246             :     // Get the indices of the blocks
     247         629 :     auto blocks = get_indices_of_blocks(blockdiagonalizing_labels);
     248             : 
     249         629 :     assert((blockdiagonalizing_labels.empty() && blocks.size() == 1) ||
     250             :            !blockdiagonalizing_labels.empty());
     251             : 
     252         629 :     SPDLOG_DEBUG("Diagonalizing the Hamiltonian with {} blocks.", blocks.size());
     253             : 
     254             :     // Diagonalize the blocks in parallel
     255         629 :     std::vector<Eigen::VectorX<real_t>> eigenenergies_blocks(blocks.size());
     256         629 :     std::vector<Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>> eigenvectors_blocks(blocks.size());
     257         629 :     oneapi::tbb::parallel_for(
     258        1969 :         oneapi::tbb::blocked_range<size_t>(0, blocks.size()), [&](const auto &range) {
     259        2680 :             for (size_t idx = range.begin(); idx != range.end(); ++idx) {
     260        1340 :                 set_task_status("Diagonalizing Hamiltonian blocks...");
     261        4385 :                 auto eigensys = min_eigenenergy.has_value() || max_eigenenergy.has_value()
     262        2070 :                     ? diagonalizer.eigh(matrix.block(blocks[idx].start, blocks[idx].start,
     263         730 :                                                      blocks[idx].size(), blocks[idx].size()),
     264             :                                         min_eigenenergy, max_eigenenergy, rtol)
     265        3290 :                     : diagonalizer.eigh(matrix.block(blocks[idx].start, blocks[idx].start,
     266        1950 :                                                      blocks[idx].size(), blocks[idx].size()),
     267             :                                         rtol);
     268        1340 :                 eigenvectors_blocks[idx] = eigensys.eigenvectors;
     269        1340 :                 eigenenergies_blocks[idx] = eigensys.eigenvalues;
     270             :             }
     271             :         });
     272             : 
     273             :     // Get the number of non-zeros per row of the combined eigenvector matrix
     274         629 :     std::vector<Eigen::Index> non_zeros_per_inner_index;
     275         629 :     non_zeros_per_inner_index.reserve(matrix.rows());
     276         629 :     Eigen::Index num_rows = 0;
     277         629 :     Eigen::Index num_cols = 0;
     278        1969 :     for (const auto &matrix : eigenvectors_blocks) {
     279        1340 :         set_task_status("Collecting eigenvectors...");
     280       37076 :         for (int i = 0; i < matrix.outerSize(); ++i) {
     281       71472 :             non_zeros_per_inner_index.push_back(matrix.outerIndexPtr()[i + 1] -
     282       35736 :                                                 matrix.outerIndexPtr()[i]);
     283             :         }
     284        1340 :         num_rows += matrix.rows();
     285        1340 :         num_cols += matrix.cols();
     286             :     }
     287             : 
     288         629 :     assert(static_cast<size_t>(num_rows) == basis->get_number_of_states());
     289         629 :     assert(static_cast<size_t>(num_cols) <= basis->get_number_of_states());
     290             : 
     291         629 :     eigenvectors.resize(num_rows, num_cols);
     292         629 :     eigenenergies.resize(num_cols, num_cols);
     293             : 
     294         629 :     if (num_cols > 0) {
     295             :         // Get the combined eigenvector matrix (in case of an restricted energy range, it is not
     296             :         // square)
     297         627 :         eigenvectors.reserve(non_zeros_per_inner_index);
     298         627 :         Eigen::Index offset_rows = 0;
     299         627 :         Eigen::Index offset_cols = 0;
     300        1965 :         for (const auto &matrix : eigenvectors_blocks) {
     301        1338 :             set_task_status("Combining eigenvectors...");
     302       36994 :             for (Eigen::Index i = 0; i < matrix.outerSize(); ++i) {
     303       35656 :                 for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     304             :                          matrix, i);
     305      660728 :                      it; ++it) {
     306      625072 :                     eigenvectors.insert(it.row() + offset_rows, it.col() + offset_cols) =
     307      625072 :                         it.value();
     308             :                 }
     309             :             }
     310        1338 :             offset_rows += matrix.rows();
     311        1338 :             offset_cols += matrix.cols();
     312             :         }
     313         627 :         eigenvectors.makeCompressed();
     314             : 
     315         627 :         assert(
     316             :             eigenvectors.nonZeros() ==
     317             :             std::accumulate(non_zeros_per_inner_index.begin(), non_zeros_per_inner_index.end(), 0));
     318             : 
     319             :         // Get the combined eigenenergy matrix
     320         627 :         eigenenergies.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     321         627 :         Eigen::Index offset = 0;
     322        1965 :         for (const auto &matrix : eigenenergies_blocks) {
     323        1338 :             set_task_status("Combining eigenenergies...");
     324       31976 :             for (int i = 0; i < matrix.size(); ++i) {
     325       30638 :                 eigenenergies.insert(i + offset, i + offset) = matrix(i);
     326             :             }
     327        1338 :             offset += matrix.size();
     328             :         }
     329         627 :         eigenenergies.makeCompressed();
     330             : 
     331             :         // Fix phase ambiguity
     332         627 :         std::vector<scalar_t> map_col_to_max(num_cols, 0);
     333       36283 :         for (int row = 0; row < eigenvectors.outerSize(); ++row) {
     334       35656 :             set_task_status("Normalizing eigenvector phases...");
     335       35656 :             for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
     336             :                      eigenvectors, row);
     337      660728 :                  it; ++it) {
     338      625072 :                 if (std::abs(it.value()) > std::abs(map_col_to_max[it.col()])) {
     339      143302 :                     map_col_to_max[it.col()] = it.value();
     340             :                 }
     341             :             }
     342             :         }
     343             : 
     344         627 :         Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> phase_matrix;
     345         627 :         phase_matrix.resize(num_cols, num_cols);
     346         627 :         phase_matrix.reserve(Eigen::VectorXi::Constant(num_cols, 1));
     347       31265 :         for (int i = 0; i < num_cols; ++i) {
     348       30638 :             phase_matrix.insert(i, i) = std::abs(map_col_to_max[i]) / map_col_to_max[i];
     349             :         }
     350         627 :         phase_matrix.makeCompressed();
     351             : 
     352         627 :         set_task_status("Applying eigenvector phases...");
     353         627 :         eigenvectors = eigenvectors * phase_matrix;
     354         627 :     }
     355             : 
     356             :     // Store the diagonalized hamiltonian
     357         629 :     matrix = eigenenergies;
     358         629 :     basis = basis->transformed(eigenvectors);
     359             : 
     360         629 :     hamiltonian_is_diagonal = true;
     361             : 
     362         629 :     return *this;
     363         629 : }
     364             : 
     365             : template <typename Derived>
     366        3950 : bool System<Derived>::is_diagonal() const {
     367        3950 :     if (hamiltonian_requires_construction) {
     368         359 :         construct_hamiltonian();
     369         359 :         hamiltonian_requires_construction = false;
     370             :     }
     371             : 
     372        3950 :     if (!hamiltonian_is_diagonal) {
     373        1114 :         real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
     374             : 
     375       39990 :         for (int row = 0; row < matrix.outerSize(); ++row) {
     376       39518 :             for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(matrix,
     377             :                                                                                            row);
     378       79057 :                  it; ++it) {
     379       40181 :                 if (it.row() != it.col() && std::abs(it.value()) > numerical_precision) {
     380         642 :                     return false;
     381             :                 }
     382             :             }
     383             :         }
     384             : 
     385         472 :         hamiltonian_is_diagonal = true;
     386             :     }
     387             : 
     388        3308 :     return true;
     389             : }
     390             : 
     391             : // Explicit instantiation
     392             : template class System<SystemAtom<double>>;
     393             : template class System<SystemAtom<std::complex<double>>>;
     394             : template class System<SystemPair<double>>;
     395             : template class System<SystemPair<std::complex<double>>>;
     396             : } // namespace pairinteraction