Line data Source code
1 : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/basis/BasisPair.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/basis/BasisPairCreator.hpp"
8 : #include "pairinteraction/database/Database.hpp"
9 : #include "pairinteraction/ket/KetAtom.hpp"
10 : #include "pairinteraction/ket/KetPair.hpp"
11 : #include "pairinteraction/system/SystemAtom.hpp"
12 : #include "pairinteraction/utils/Range.hpp"
13 : #include "pairinteraction/utils/tensor.hpp"
14 :
15 : #include <cassert>
16 : #include <memory>
17 : #include <oneapi/tbb.h>
18 : #include <vector>
19 :
20 : namespace pairinteraction {
21 : template <typename Scalar>
22 103 : BasisPair<Scalar>::BasisPair(Private /*unused*/, ketvec_t &&kets,
23 : map_range_t &&map_range_of_state_index2,
24 : map_indices_t &&state_indices_to_ket_index,
25 : std::shared_ptr<const BasisAtom<Scalar>> basis1,
26 : std::shared_ptr<const BasisAtom<Scalar>> basis2)
27 103 : : Basis<BasisPair<Scalar>>(std::move(kets)),
28 103 : map_range_of_state_index2(std::move(map_range_of_state_index2)),
29 103 : state_indices_to_ket_index(std::move(state_indices_to_ket_index)), basis1(std::move(basis1)),
30 206 : basis2(std::move(basis2)) {}
31 :
32 : template <typename Scalar>
33 : const typename BasisPair<Scalar>::range_t &
34 3263918 : BasisPair<Scalar>::get_index_range(size_t state_index1) const {
35 3263918 : return map_range_of_state_index2.at(state_index1);
36 : }
37 :
38 : template <typename Scalar>
39 108 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis1() const {
40 108 : return basis1;
41 : }
42 :
43 : template <typename Scalar>
44 108 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis2() const {
45 108 : return basis2;
46 : }
47 :
48 : template <typename Scalar>
49 4063506 : int BasisPair<Scalar>::get_ket_index_from_tuple(size_t state_index1, size_t state_index2) const {
50 4063506 : if (state_indices_to_ket_index.count({state_index1, state_index2}) == 0) {
51 280669 : return -1;
52 : }
53 3782837 : return state_indices_to_ket_index.at({state_index1, state_index2});
54 : }
55 :
56 : template <typename Scalar>
57 : Eigen::VectorX<Scalar>
58 34 : BasisPair<Scalar>::get_amplitudes(std::shared_ptr<const KetAtom> ket1,
59 : std::shared_ptr<const KetAtom> ket2) const {
60 34 : return get_amplitudes(basis1->get_canonical_state_from_ket(ket1),
61 34 : basis2->get_canonical_state_from_ket(ket2))
62 102 : .transpose();
63 : }
64 :
65 : template <typename Scalar>
66 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
67 34 : BasisPair<Scalar>::get_amplitudes(std::shared_ptr<const BasisAtom<Scalar>> other1,
68 : std::shared_ptr<const BasisAtom<Scalar>> other2) const {
69 68 : if (other1->get_id_of_kets() != basis1->get_id_of_kets() ||
70 34 : other2->get_id_of_kets() != basis2->get_id_of_kets()) {
71 0 : throw std::invalid_argument("The other objects must be expressed using the same kets.");
72 : }
73 :
74 34 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
75 :
76 34 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> coefficients1 =
77 34 : basis1->get_coefficients().adjoint() * other1->get_coefficients();
78 34 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> coefficients2 =
79 34 : basis2->get_coefficients().adjoint() * other2->get_coefficients();
80 :
81 34 : oneapi::tbb::concurrent_vector<Eigen::Triplet<Scalar>> triplets;
82 :
83 : // Loop over the rows of the first coefficient matrix
84 34 : oneapi::tbb::parallel_for(
85 34 : oneapi::tbb::blocked_range<Eigen::Index>(0, coefficients1.outerSize()),
86 4861 : [&](const auto &range) {
87 10281 : for (Eigen::Index row1 = range.begin(); row1 != range.end(); ++row1) {
88 :
89 5454 : const auto &range_row2 = this->get_index_range(row1);
90 :
91 : // Loop over the rows of the second coefficient matrix that are energetically
92 : // allowed
93 54596 : for (auto row2 = static_cast<Eigen::Index>(range_row2.min());
94 54596 : row2 < static_cast<Eigen::Index>(range_row2.max()); ++row2) {
95 :
96 49142 : Eigen::Index row = get_ket_index_from_tuple(row1, row2);
97 49142 : if (row < 0) {
98 709 : continue;
99 : }
100 :
101 : // Loop over the non-zero column elements of the first coefficient matrix
102 49306 : for (typename Eigen::SparseMatrix<Scalar, Eigen::RowMajor>::InnerIterator it1(
103 : coefficients1, row1);
104 49306 : it1; ++it1) {
105 :
106 873 : Eigen::Index col1 = it1.col();
107 873 : Scalar value1 = it1.value();
108 :
109 : // Loop over the non-zero column elements of the second coefficient matrix
110 924 : for (typename Eigen::SparseMatrix<Scalar, Eigen::RowMajor>::InnerIterator
111 873 : it2(coefficients2, row2);
112 924 : it2; ++it2) {
113 :
114 51 : Eigen::Index col2 = it2.col();
115 51 : Scalar value2 = it2.value();
116 51 : Eigen::Index col = col1 * coefficients2.cols() + col2;
117 :
118 : // Store the entry
119 32 : Scalar value = value1 * value2;
120 51 : if (std::abs(value) > numerical_precision) {
121 51 : triplets.emplace_back(row, col, value);
122 : }
123 : }
124 : }
125 : }
126 : }
127 : });
128 :
129 : // Construct the combined matrix from the triplets
130 34 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> matrix(this->get_number_of_kets(),
131 34 : other1->get_number_of_states() *
132 34 : other2->get_number_of_states());
133 34 : matrix.setFromTriplets(triplets.begin(), triplets.end());
134 34 : matrix.makeCompressed();
135 :
136 68 : return matrix.adjoint() * this->get_coefficients();
137 34 : }
138 :
139 : template <typename Scalar>
140 : Eigen::VectorX<typename BasisPair<Scalar>::real_t>
141 34 : BasisPair<Scalar>::get_overlaps(std::shared_ptr<const KetAtom> ket1,
142 : std::shared_ptr<const KetAtom> ket2) const {
143 34 : return get_amplitudes(ket1, ket2).cwiseAbs2();
144 : }
145 :
146 : template <typename Scalar>
147 : Eigen::SparseMatrix<typename BasisPair<Scalar>::real_t, Eigen::RowMajor>
148 0 : BasisPair<Scalar>::get_overlaps(std::shared_ptr<const BasisAtom<Scalar>> other1,
149 : std::shared_ptr<const BasisAtom<Scalar>> other2) const {
150 0 : return get_amplitudes(other1, other2).cwiseAbs2();
151 : }
152 :
153 : template <typename Scalar>
154 0 : Eigen::VectorX<Scalar> BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const ket_t> /*ket*/,
155 : OperatorType /*type*/,
156 : int /*q*/) const {
157 0 : throw std::invalid_argument("It is required to specify one operator for each atom.");
158 : }
159 :
160 : template <typename Scalar>
161 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
162 0 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const Type> /*final*/, OperatorType /*type*/,
163 : int /*q*/) const {
164 0 : throw std::invalid_argument("It is required to specify one operator for each atom.");
165 : }
166 :
167 : template <typename Scalar>
168 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
169 11 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const Type> final, OperatorType type1,
170 : OperatorType type2, int q1, int q2) const {
171 11 : auto initial1 = this->get_basis1();
172 11 : auto initial2 = this->get_basis2();
173 11 : auto final1 = final->get_basis1();
174 11 : auto final2 = final->get_basis2();
175 :
176 11 : auto matrix_elements1 =
177 11 : initial1->get_database().get_matrix_elements(initial1, final1, type1, q1);
178 11 : auto matrix_elements2 =
179 11 : initial2->get_database().get_matrix_elements(initial2, final2, type2, q2);
180 11 : auto matrix_elements = utils::calculate_tensor_product(this->shared_from_this(), final,
181 : matrix_elements1, matrix_elements2);
182 11 : assert(static_cast<size_t>(matrix_elements.rows()) == final->get_number_of_kets());
183 11 : assert(static_cast<size_t>(matrix_elements.cols()) == this->get_number_of_kets());
184 :
185 22 : return final->get_coefficients().adjoint() * matrix_elements * this->get_coefficients();
186 11 : }
187 :
188 : template <typename Scalar>
189 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
190 1 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const BasisAtom<Scalar>> final1,
191 : std::shared_ptr<const BasisAtom<Scalar>> final2,
192 : OperatorType type1, OperatorType type2, int q1,
193 : int q2) const {
194 : // Construct a pair basis with the two single-atom bases
195 1 : auto system1 = SystemAtom<Scalar>(final1);
196 1 : auto system2 = SystemAtom<Scalar>(final2);
197 1 : auto final = BasisPairCreator<Scalar>().add(system1).add(system2).create();
198 1 : assert(final->get_number_of_states() ==
199 : final1->get_number_of_states() * final2->get_number_of_states());
200 1 : assert(final->get_number_of_kets() ==
201 : final1->get_number_of_states() * final2->get_number_of_states());
202 :
203 2 : return this->get_matrix_elements(final, type1, type2, q1, q2);
204 1 : }
205 :
206 : template <typename Scalar>
207 : Eigen::VectorX<Scalar>
208 1 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const ket_t> ket, OperatorType type1,
209 : OperatorType type2, int q1, int q2) const {
210 : // Construct a pair basis containing only the pair ket
211 1 : auto final = this->get_canonical_state_from_ket(ket);
212 1 : assert(final->get_number_of_states() == 1);
213 :
214 3 : return this->get_matrix_elements(final, type1, type2, q1, q2).row(0);
215 1 : }
216 :
217 : template <typename Scalar>
218 : Eigen::VectorX<Scalar>
219 1 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const KetAtom> ket1,
220 : std::shared_ptr<const KetAtom> ket2, OperatorType type1,
221 : OperatorType type2, int q1, int q2) const {
222 : // Construct a pair basis with the two single-atom kets
223 1 : auto final1 = this->get_basis1()->get_canonical_state_from_ket(ket1);
224 1 : auto final2 = this->get_basis2()->get_canonical_state_from_ket(ket2);
225 1 : auto system1 = SystemAtom<Scalar>(final1);
226 1 : auto system2 = SystemAtom<Scalar>(final2);
227 1 : auto final = BasisPairCreator<Scalar>().add(system1).add(system2).create();
228 1 : assert(final->get_number_of_states() == 1);
229 1 : assert(final->get_number_of_kets() == 1);
230 :
231 3 : return this->get_matrix_elements(final, type1, type2, q1, q2).row(0);
232 1 : }
233 :
234 : // Explicit instantiations
235 : template class BasisPair<double>;
236 : template class BasisPair<std::complex<double>>;
237 : } // namespace pairinteraction
|
