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