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/SorterType.hpp"
9 : #include "pairinteraction/interfaces/DiagonalizerInterface.hpp"
10 : #include "pairinteraction/ket/KetAtom.hpp"
11 : #include "pairinteraction/ket/KetPair.hpp"
12 : #include "pairinteraction/system/SystemAtom.hpp"
13 : #include "pairinteraction/system/SystemPair.hpp"
14 : #include "pairinteraction/utils/TaskControl.hpp"
15 : #include "pairinteraction/utils/eigen_assertion.hpp"
16 : #include "pairinteraction/utils/eigen_compat.hpp"
17 :
18 : #include <Eigen/SparseCore>
19 : #include <algorithm>
20 : #include <complex>
21 : #include <limits>
22 : #include <memory>
23 : #include <numeric>
24 : #include <oneapi/tbb.h>
25 : #include <optional>
26 : #include <spdlog/spdlog.h>
27 : #include <stdexcept>
28 : #include <vector>
29 :
30 : namespace pairinteraction {
31 : template <typename Derived>
32 1100 : System<Derived>::System(std::shared_ptr<const basis_t> basis)
33 1100 : : basis(std::move(basis)),
34 1100 : matrix(static_cast<Eigen::Index>(this->basis->get_number_of_states()),
35 3300 : static_cast<Eigen::Index>(this->basis->get_number_of_states())) {
36 : // The unperturbed Hamiltonian is constructed from the energies of the kets
37 150031 : for (const auto &ket : this->basis->get_kets()) {
38 148931 : if (static_cast<real_t>(ket->get_energy()) == std::numeric_limits<real_t>::max()) {
39 0 : throw std::invalid_argument("The kets of the basis do not have a well-defined energy.");
40 : }
41 : }
42 1100 : }
43 :
44 : template <typename Derived>
45 602 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_basis() const {
46 602 : if (hamiltonian_requires_construction) {
47 87 : construct_hamiltonian();
48 87 : hamiltonian_requires_construction = false;
49 : }
50 602 : return basis;
51 : }
52 :
53 : template <typename Derived>
54 1251 : std::shared_ptr<const typename System<Derived>::basis_t> System<Derived>::get_eigenbasis() const {
55 1251 : if (hamiltonian_requires_construction) {
56 134 : construct_hamiltonian();
57 134 : hamiltonian_requires_construction = false;
58 : }
59 1251 : if (!this->is_diagonal()) {
60 0 : throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
61 : }
62 1251 : return basis;
63 : }
64 :
65 : template <typename Derived>
66 931 : Eigen::VectorX<typename System<Derived>::real_t> System<Derived>::get_eigenenergies() const {
67 931 : if (hamiltonian_requires_construction) {
68 0 : construct_hamiltonian();
69 0 : hamiltonian_requires_construction = false;
70 : }
71 931 : if (!this->is_diagonal()) {
72 0 : throw std::runtime_error("The Hamiltonian has not been diagonalized yet.");
73 : }
74 1404 : return matrix.diagonal().real();
75 : }
76 :
77 : template <typename Derived>
78 : const Eigen::SparseMatrix<typename System<Derived>::scalar_t, Eigen::RowMajor> &
79 223 : System<Derived>::get_matrix() const {
80 223 : if (hamiltonian_requires_construction) {
81 177 : construct_hamiltonian();
82 177 : hamiltonian_requires_construction = false;
83 : }
84 223 : return matrix;
85 : }
86 :
87 : template <typename Derived>
88 : Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>
89 1242 : System<Derived>::get_sorter(const std::vector<SorterType> &labels) const {
90 : // The first sorter type in labels is supposed to be the primary sorting criterion.
91 1242 : if (hamiltonian_requires_construction) {
92 6 : construct_hamiltonian();
93 6 : hamiltonian_requires_construction = false;
94 : }
95 :
96 : // Drop repeated sorter types. This also ensures that energy occurs at most once.
97 1242 : std::vector<SorterType> unique_labels;
98 1242 : unique_labels.reserve(labels.size());
99 2486 : for (const auto &label : labels) {
100 1271 : if (std::find(unique_labels.begin(), unique_labels.end(), label) == unique_labels.end()) {
101 1241 : unique_labels.push_back(label);
102 : }
103 : }
104 :
105 1211 : auto it = std::find(unique_labels.begin(), unique_labels.end(), SorterType::ENERGY);
106 1215 : std::vector<SorterType> before_energy(unique_labels.begin(), it);
107 1216 : bool contains_energy = (it != unique_labels.end());
108 1221 : std::vector<SorterType> after_energy(contains_energy ? it + 1 : unique_labels.end(),
109 : unique_labels.end());
110 :
111 1229 : Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> sorter(matrix.rows());
112 1217 : sorter.setIdentity();
113 :
114 1252 : if (!after_energy.empty()) {
115 2 : basis->get_sorter_without_checks(after_energy, sorter);
116 : }
117 :
118 1238 : if (contains_energy) {
119 647 : std::vector<real_t> energies_of_states;
120 647 : energies_of_states.reserve(matrix.rows());
121 31087 : for (int i = 0; i < matrix.rows(); ++i) {
122 30476 : energies_of_states.push_back(std::real(matrix.coeff(i, i)));
123 : }
124 :
125 652 : std::stable_sort(
126 631 : sorter.indices().data(), sorter.indices().data() + sorter.indices().size(),
127 128900 : [&](int i, int j) { return energies_of_states[i] < energies_of_states[j]; });
128 664 : }
129 :
130 1248 : if (!before_energy.empty()) {
131 588 : basis->get_sorter_without_checks(before_energy, sorter);
132 : }
133 :
134 2451 : return sorter;
135 1236 : }
136 :
137 : template <typename Derived>
138 : std::vector<IndicesOfBlock>
139 657 : System<Derived>::get_indices_of_blocks(const std::vector<SorterType> &labels) const {
140 657 : if (hamiltonian_requires_construction) {
141 0 : construct_hamiltonian();
142 0 : hamiltonian_requires_construction = false;
143 : }
144 :
145 657 : std::set<SorterType> unique_labels(labels.begin(), labels.end());
146 649 : basis->perform_blocks_checks(unique_labels);
147 :
148 650 : IndicesOfBlocksCreator blocks_creator({0, static_cast<size_t>(matrix.rows())});
149 :
150 645 : if (!unique_labels.empty()) {
151 559 : basis->get_indices_of_blocks_without_checks(unique_labels, blocks_creator);
152 : }
153 :
154 1307 : return blocks_creator.create();
155 645 : }
156 :
157 : template <typename Derived>
158 : System<Derived> &
159 0 : System<Derived>::transform(const Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> &transformation) {
160 0 : if (hamiltonian_requires_construction) {
161 0 : construct_hamiltonian();
162 0 : hamiltonian_requires_construction = false;
163 : }
164 :
165 0 : set_task_status("Applying transformation...");
166 0 : if (matrix.cols() != 0) {
167 0 : matrix = transformation.adjoint() * matrix * transformation;
168 0 : basis = basis->transformed(transformation);
169 : }
170 :
171 0 : hamiltonian_is_diagonal = false;
172 :
173 : // A transformed system might have lost its block-diagonalizability if the
174 : // transformation was not a sorting
175 0 : blockdiagonalizing_labels.clear();
176 :
177 0 : return *this;
178 : }
179 :
180 : template <typename Derived>
181 : System<Derived> &
182 642 : System<Derived>::transform(const Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> &sorter) {
183 642 : if (hamiltonian_requires_construction) {
184 0 : construct_hamiltonian();
185 0 : hamiltonian_requires_construction = false;
186 : }
187 :
188 642 : set_task_status("Applying transformation...");
189 666 : if (matrix.cols() != 0) {
190 664 : matrix = matrix.twistedBy(sorter.inverse());
191 661 : basis = basis->transformed(sorter);
192 : }
193 :
194 662 : return *this;
195 : }
196 :
197 : template <typename Derived>
198 688 : System<Derived> &System<Derived>::diagonalize(const DiagonalizerInterface<scalar_t> &diagonalizer,
199 : std::optional<real_t> min_eigenenergy,
200 : std::optional<real_t> max_eigenenergy, double rtol,
201 : bool sort_by_energy) {
202 688 : set_task_status("Preparing Hamiltonian...");
203 :
204 699 : if (hamiltonian_requires_construction) {
205 695 : construct_hamiltonian();
206 695 : hamiltonian_requires_construction = false;
207 : }
208 :
209 699 : if (this->is_diagonal()) {
210 40 : if (sort_by_energy && !basis->is_sorted_by_energy()) {
211 2 : transform(get_sorter({SorterType::ENERGY}));
212 : }
213 40 : return *this;
214 : }
215 :
216 641 : Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenvectors;
217 632 : Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> eigenenergies;
218 :
219 : // Sort the Hamiltonian according to the block structure
220 652 : if (!blockdiagonalizing_labels.empty()) {
221 551 : auto sorter = get_sorter(blockdiagonalizing_labels);
222 574 : matrix = matrix.twistedBy(sorter.inverse());
223 590 : basis = basis->transformed(sorter);
224 577 : }
225 :
226 : // Get the indices of the blocks
227 648 : auto blocks = get_indices_of_blocks(blockdiagonalizing_labels);
228 :
229 648 : assert((blockdiagonalizing_labels.empty() && blocks.size() == 1) ||
230 : !blockdiagonalizing_labels.empty());
231 :
232 628 : SPDLOG_DEBUG("Diagonalizing the Hamiltonian with {} blocks.", blocks.size());
233 :
234 : // Diagonalize the blocks in parallel
235 659 : std::vector<Eigen::VectorX<real_t>> eigenenergies_blocks(blocks.size());
236 659 : std::vector<Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>> eigenvectors_blocks(blocks.size());
237 658 : oneapi::tbb::parallel_for(
238 2004 : oneapi::tbb::blocked_range<size_t>(0, blocks.size()), [&](const auto &range) {
239 2670 : for (size_t idx = range.begin(); idx != range.end(); ++idx) {
240 1347 : set_task_status("Diagonalizing Hamiltonian blocks...");
241 4457 : auto eigensys = min_eigenenergy.has_value() || max_eigenenergy.has_value()
242 2095 : ? diagonalizer.eigh(matrix.block(blocks[idx].start, blocks[idx].start,
243 724 : blocks[idx].size(), blocks[idx].size()),
244 : min_eigenenergy, max_eigenenergy, rtol)
245 3367 : : diagonalizer.eigh(matrix.block(blocks[idx].start, blocks[idx].start,
246 2016 : blocks[idx].size(), blocks[idx].size()),
247 : rtol);
248 1365 : eigenvectors_blocks[idx] = eigensys.eigenvectors;
249 1344 : eigenenergies_blocks[idx] = eigensys.eigenvalues;
250 : }
251 : });
252 :
253 : // Get the number of non-zeros per row of the combined eigenvector matrix
254 655 : std::vector<Eigen::Index> non_zeros_per_inner_index;
255 655 : non_zeros_per_inner_index.reserve(matrix.rows());
256 644 : Eigen::Index num_rows = 0;
257 644 : Eigen::Index num_cols = 0;
258 2017 : for (const auto &matrix : eigenvectors_blocks) {
259 1337 : set_task_status("Collecting eigenvectors...");
260 37964 : for (int i = 0; i < matrix.outerSize(); ++i) {
261 73259 : non_zeros_per_inner_index.push_back(matrix.outerIndexPtr()[i + 1] -
262 36628 : matrix.outerIndexPtr()[i]);
263 : }
264 1295 : num_rows += matrix.rows();
265 1372 : num_cols += matrix.cols();
266 : }
267 :
268 653 : assert(static_cast<size_t>(num_rows) == basis->get_number_of_states());
269 653 : assert(static_cast<size_t>(num_cols) <= basis->get_number_of_states());
270 :
271 652 : eigenvectors.resize(num_rows, num_cols);
272 650 : eigenenergies.resize(num_cols, num_cols);
273 :
274 653 : if (num_cols > 0) {
275 : // Get the combined eigenvector matrix (in case of an restricted energy range, it is not
276 : // square)
277 652 : eigenvectors.reserve(non_zeros_per_inner_index);
278 648 : Eigen::Index offset_rows = 0;
279 648 : Eigen::Index offset_cols = 0;
280 2019 : for (const auto &matrix : eigenvectors_blocks) {
281 1352 : set_task_status("Combining eigenvectors...");
282 38687 : for (Eigen::Index i = 0; i < matrix.outerSize(); ++i) {
283 37310 : for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
284 : matrix, i);
285 850983 : it; ++it) {
286 815003 : eigenvectors.insert(it.row() + offset_rows, it.col() + offset_cols) =
287 813726 : it.value();
288 : }
289 : }
290 1361 : offset_rows += matrix.rows();
291 1369 : offset_cols += matrix.cols();
292 : }
293 653 : eigenvectors.makeCompressed();
294 :
295 654 : assert(
296 : eigenvectors.nonZeros() ==
297 : std::accumulate(non_zeros_per_inner_index.begin(), non_zeros_per_inner_index.end(), 0));
298 :
299 : // Get the combined eigenenergy matrix
300 653 : eigenenergies.reserve(Eigen::VectorXi::Constant(num_cols, 1));
301 656 : Eigen::Index offset = 0;
302 2026 : for (const auto &matrix : eigenenergies_blocks) {
303 1346 : set_task_status("Combining eigenenergies...");
304 33159 : for (int i = 0; i < matrix.size(); ++i) {
305 31909 : eigenenergies.insert(i + offset, i + offset) = matrix(i);
306 : }
307 1328 : offset += matrix.size();
308 : }
309 654 : eigenenergies.makeCompressed();
310 :
311 : // Fix phase ambiguity
312 653 : std::vector<scalar_t> map_col_to_max(num_cols, 0);
313 38011 : for (int row = 0; row < eigenvectors.outerSize(); ++row) {
314 37358 : set_task_status("Normalizing eigenvector phases...");
315 37533 : for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
316 : eigenvectors, row);
317 849488 : it; ++it) {
318 809183 : if (std::abs(it.value()) > std::abs(map_col_to_max[it.col()])) {
319 164530 : map_col_to_max[it.col()] = it.value();
320 : }
321 : }
322 : }
323 :
324 643 : Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> phase_matrix;
325 657 : phase_matrix.resize(num_cols, num_cols);
326 657 : phase_matrix.reserve(Eigen::VectorXi::Constant(num_cols, 1));
327 32195 : for (int i = 0; i < num_cols; ++i) {
328 31553 : phase_matrix.insert(i, i) = std::abs(map_col_to_max[i]) / map_col_to_max[i];
329 : }
330 642 : phase_matrix.makeCompressed();
331 :
332 657 : set_task_status("Applying eigenvector phases...");
333 657 : eigenvectors = eigenvectors * phase_matrix;
334 657 : }
335 :
336 : // Store the diagonalized hamiltonian and assign the eigenenergies to the eigenstates
337 653 : matrix = eigenenergies;
338 648 : std::vector<real_t> energy_of_eigenstates;
339 648 : energy_of_eigenstates.reserve(num_cols);
340 1982 : for (const auto &eigenenergies_block : eigenenergies_blocks) {
341 1331 : energy_of_eigenstates.insert(energy_of_eigenstates.end(), eigenenergies_block.data(),
342 1338 : eigenenergies_block.data() + eigenenergies_block.size());
343 : }
344 639 : basis = basis->transformed(eigenvectors, energy_of_eigenstates);
345 :
346 644 : hamiltonian_is_diagonal = true;
347 644 : if (sort_by_energy) {
348 637 : transform(get_sorter({SorterType::ENERGY}));
349 : }
350 :
351 659 : return *this;
352 659 : }
353 :
354 : template <typename Derived>
355 2882 : bool System<Derived>::is_diagonal() const {
356 2882 : if (hamiltonian_requires_construction) {
357 0 : construct_hamiltonian();
358 0 : hamiltonian_requires_construction = false;
359 : }
360 :
361 2882 : if (!hamiltonian_is_diagonal) {
362 898 : real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
363 :
364 29656 : for (int row = 0; row < matrix.outerSize(); ++row) {
365 29403 : for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(matrix,
366 : row);
367 59329 : it; ++it) {
368 30280 : if (it.row() != it.col() && std::abs(it.value()) > numerical_precision) {
369 629 : return false;
370 : }
371 : }
372 : }
373 :
374 : // The states of the basis are eigenstates, so assign them the diagonal entries of the
375 : // Hamiltonian as energies
376 237 : Eigen::VectorX<real_t> diagonal = matrix.diagonal().real();
377 244 : Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> identity(matrix.rows(), matrix.cols());
378 244 : identity.setIdentity();
379 487 : basis = basis->transformed(
380 486 : identity, std::vector<real_t>(diagonal.data(), diagonal.data() + diagonal.size()));
381 :
382 244 : hamiltonian_is_diagonal = true;
383 244 : }
384 :
385 2225 : return true;
386 : }
387 :
388 : // Explicit instantiation
389 : template class System<SystemAtom<double>>;
390 : template class System<SystemAtom<std::complex<double>>>;
391 : template class System<SystemPair<double>>;
392 : template class System<SystemPair<std::complex<double>>>;
393 : } // namespace pairinteraction
|
