pairinteraction
A Rydberg Interaction Calculator
GreenTensor.cpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2025 Pairinteraction Developers
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
5
10
11#include <Eigen/Dense>
12#include <Eigen/Sparse>
13#include <complex>
14#include <map>
15#include <spdlog/spdlog.h>
16#include <unsupported/Eigen/Splines>
17
18namespace pairinteraction {
19template <typename Scalar>
20GreenTensor<Scalar>::ConstantEntry::ConstantEntry(int row, int col, Scalar val)
21 : row_(row), col_(col), val_(val) {}
22
23template <typename Scalar>
25 return val_;
26}
27
28template <typename Scalar>
30 return row_;
31}
32
33template <typename Scalar>
35 return col_;
36}
37
38template <typename Scalar>
40 Eigen::Spline<real_t, 1> real_spline,
41 Eigen::Spline<real_t, 1> imag_spline)
42 : row_(row), col_(col), real_spline(std::move(real_spline)),
43 imag_spline(std::move(imag_spline)) {}
44
45template <typename Scalar>
48 return {real_spline(omega)(0), imag_spline(omega)(0)};
49 } else {
50 return real_spline(omega)(0);
51 }
52}
53
54template <typename Scalar>
56 return row_;
57}
58
59template <typename Scalar>
61 return col_;
62}
63
64template <typename Scalar>
66 int kappa1, int kappa2, const Eigen::MatrixX<Scalar> &tensor_in_cartesian_coordinates) {
67
68 const real_t scale = tensor_in_cartesian_coordinates.norm();
69 const real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
70
71 Eigen::SparseMatrix<complex_t> tensor =
72 (spherical::get_transformator<complex_t>(kappa1) * tensor_in_cartesian_coordinates *
73 spherical::get_transformator<complex_t>(kappa2).adjoint())
74 .sparseView(1, numerical_precision);
75
76 std::vector<Entry> entries;
77 for (int k = 0; k < tensor.outerSize(); ++k) {
78 for (typename Eigen::SparseMatrix<complex_t>::InnerIterator it(tensor, k); it; ++it) {
80 entries.emplace_back(ConstantEntry(it.row(), it.col(), it.value()));
81 } else {
82 entries.emplace_back(ConstantEntry(it.row(), it.col(), it.value().real()));
83 assert(abs(it.value().imag()) < numerical_precision);
84 }
85 }
86 }
87 entries_map[{kappa1, kappa2}] = std::move(entries);
88}
89
90template <typename Scalar>
92 int kappa1, int kappa2,
93 const std::vector<Eigen::MatrixX<Scalar>> &tensors_in_cartesian_coordinates,
94 const std::vector<double> &omegas) {
95
96 if (tensors_in_cartesian_coordinates.size() != omegas.size()) {
97 throw std::invalid_argument("The number of tensors and omegas must match.");
98 }
99
100 if (tensors_in_cartesian_coordinates.size() < 4) {
101 throw std::invalid_argument(
102 "At least 4 tensors are required for the applied cubic spline interpolation.");
103 }
104
105 auto num_knots = static_cast<int>(omegas.size());
106 Eigen::Map<const Eigen::RowVectorXd> knots(omegas.data(), num_knots);
107
108 constexpr int spline_degree = 3; // cubic spline interpolation
109
110 // Temporary storage wih key = (row, col) and value = vector of one double per omega
111 std::map<std::pair<int, int>, std::pair<Eigen::RowVectorXd, Eigen::RowVectorXd>> temp_map;
112 for (int idx = 0; idx < num_knots; ++idx) {
113
114 const real_t scale = tensors_in_cartesian_coordinates[idx].norm();
115 const real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
116
117 Eigen::SparseMatrix<complex_t> tensor =
118 (spherical::get_transformator<complex_t>(kappa1) *
119 tensors_in_cartesian_coordinates[idx] *
120 spherical::get_transformator<complex_t>(kappa2).adjoint())
121 .sparseView(1, numerical_precision);
122
123 for (int k = 0; k < tensor.outerSize(); ++k) {
124 for (typename Eigen::SparseMatrix<complex_t>::InnerIterator it(tensor, k); it; ++it) {
125 std::pair<int, int> key{it.row(), it.col()};
126 auto &[vec_real, vec_imag] =
127 temp_map
128 .try_emplace(key, Eigen::RowVectorXd::Zero(num_knots),
129 Eigen::RowVectorXd::Zero(num_knots))
130 .first->second;
131 vec_real(idx) = it.value().real();
133 vec_imag(idx) = it.value().imag();
134 } else {
135 assert(abs(it.value().imag()) < numerical_precision);
136 }
137 }
138 }
139 }
140
141 // Set the green tensor entries with spline interpolation
142 std::vector<Entry> entries;
143 entries.reserve(temp_map.size());
144 for (const auto &[key, value] : temp_map) {
145 const auto &[vec_real, vec_imag] = value;
146 const auto &[row, col] = key;
147
148 Eigen::Spline<real_t, 1> real_spline =
149 Eigen::SplineFitting<Eigen::Spline<real_t, 1>>::Interpolate(vec_real, spline_degree,
150 knots);
151
152 Eigen::Spline<real_t, 1> imag_spline;
154 imag_spline = Eigen::SplineFitting<Eigen::Spline<real_t, 1>>::Interpolate(
155 vec_imag, spline_degree, knots);
156 }
157
158 entries.emplace_back(
159 OmegaDependentEntry(row, col, std::move(real_spline), std::move(imag_spline)));
160 }
161 entries_map[{kappa1, kappa2}] = std::move(entries);
162}
163
164template <typename Scalar>
165const std::vector<typename GreenTensor<Scalar>::Entry> &
166GreenTensor<Scalar>::get_spherical_entries(int kappa1, int kappa2) const {
167 if (auto it = entries_map.find({kappa1, kappa2}); it != entries_map.end()) {
168 return it->second;
169 }
170 static const std::vector<Entry> empty_entries;
171 return empty_entries;
172}
173
174// Explicit instantiations
175template class GreenTensor<double>;
177} // namespace pairinteraction
typename traits::NumTraits< Scalar >::real_t real_t
Definition: GreenTensor.hpp:23
void create_entries_from_cartesian(int kappa1, int kappa2, const Eigen::MatrixX< Scalar > &tensor_in_cartesian_coordinates)
Definition: GreenTensor.cpp:65
const std::vector< Entry > & get_spherical_entries(int kappa1, int kappa2) const
Matrix< Type, Dynamic, Dynamic > MatrixX
Helper struct to extract types from a numerical type.
Definition: traits.hpp:35