18#include <spdlog/spdlog.h>
21template <
typename Scalar>
25template <
typename Scalar>
27 this->hamiltonian_requires_construction =
true;
30 throw std::invalid_argument(
31 "The field must not have a y-component if the scalar type is real.");
34 electric_field = field;
39template <
typename Scalar>
41 this->hamiltonian_requires_construction =
true;
44 throw std::invalid_argument(
45 "The field must not have a y-component if the scalar type is real.");
48 magnetic_field = field;
53template <
typename Scalar>
55 this->hamiltonian_requires_construction =
true;
56 diamagnetism_enabled = enable;
60template <
typename Scalar>
63 this->hamiltonian_requires_construction =
true;
66 throw std::invalid_argument(
67 "The distance vector must not have a y-component if the scalar type is real.");
70 ion_distance_vector = vector;
75template <
typename Scalar>
77 this->hamiltonian_requires_construction =
true;
82template <
typename Scalar>
84 this->hamiltonian_requires_construction =
true;
85 if (value < 2 || value > 3) {
86 throw std::invalid_argument(
"The order of the Rydberg-ion interaction must be 2 or 3");
88 ion_interaction_order = value;
92template <
typename Scalar>
94 auto basis = this->hamiltonian->
get_basis();
98 this->hamiltonian_is_diagonal =
true;
99 bool sort_by_quantum_number_f =
true;
100 bool sort_by_quantum_number_m =
true;
101 bool sort_by_parity =
true;
106 real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
108 real_t typical_magnetic_dipole = 1e2;
109 real_t typical_electric_dipole = 1e4;
110 real_t typical_electric_quadrupole = 1e8;
115 Eigen::Map<const Eigen::Vector3<real_t>> vector_map(ion_distance_vector.data(),
116 ion_distance_vector.size());
117 real_t distance = vector_map.norm();
118 SPDLOG_DEBUG(
"Distance to the ion: {}", distance);
121 if (std::isfinite(distance) && ion_interaction_order >= 2) {
125 spherical::get_transformator<Scalar>(1) * vector_map / distance;
126 vector_dipole_order = vector_dipole_order.conjugate() / std::pow(distance, 2);
128 for (
int q = -1; q <= 1; ++q) {
129 if (std::abs(vector_dipole_order[q + 1]) * typical_electric_dipole >
130 numerical_precision) {
131 *this->hamiltonian -= ion_charge * vector_dipole_order[q + 1] *
133 this->hamiltonian_is_diagonal =
false;
134 sort_by_quantum_number_f =
false;
135 sort_by_parity =
false;
136 sort_by_quantum_number_m &= (q == 0);
143 if (std::isfinite(distance) && ion_interaction_order >= 3) {
148 Eigen::KroneckerProduct(vector_map / distance, vector_map / distance);
149 vector_quadrupole_order = 3 * vector_quadrupole_order.conjugate() / std::pow(distance, 3);
151 for (
int q = -2; q <= 2; ++q) {
152 if (std::abs(vector_quadrupole_order[q + 2]) * typical_electric_quadrupole >
153 numerical_precision) {
154 *this->hamiltonian -= ion_charge * vector_quadrupole_order[q + 2] *
156 this->hamiltonian_is_diagonal =
false;
157 sort_by_quantum_number_f =
false;
158 sort_by_quantum_number_m &= (q == 0);
167 Eigen::Map<const Eigen::Vector3<real_t>>(electric_field.data(), electric_field.size());
168 electric_field_spherical = electric_field_spherical.conjugate();
172 Eigen::Map<const Eigen::Vector3<real_t>>(magnetic_field.data(), magnetic_field.size());
173 magnetic_field_spherical = magnetic_field_spherical.conjugate();
179 for (
int q = -1; q <= 1; ++q) {
180 if (std::abs(electric_field_spherical[q + 1]) * typical_electric_dipole >
181 numerical_precision) {
182 *this->hamiltonian -= electric_field_spherical[q + 1] *
184 this->hamiltonian_is_diagonal =
false;
185 sort_by_quantum_number_f =
false;
186 sort_by_parity =
false;
187 sort_by_quantum_number_m &= (q == 0);
195 for (
int q = -1; q <= 1; ++q) {
196 if (std::abs(magnetic_field_spherical[q + 1]) * typical_magnetic_dipole >
197 numerical_precision) {
198 *this->hamiltonian -= magnetic_field_spherical[q + 1] *
200 this->hamiltonian_is_diagonal =
false;
201 sort_by_quantum_number_f =
false;
202 sort_by_quantum_number_m &= (q == 0);
213 if (diamagnetism_enabled) {
214 if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
215 numerical_precision) {
216 *this->hamiltonian +=
static_cast<real_t>(1 / 12.) *
217 static_cast<Scalar
>(std::pow(magnetic_field_spherical[1], 2)) *
219 *this->hamiltonian -=
static_cast<real_t>(1 / 12.) *
220 static_cast<Scalar
>(std::pow(magnetic_field_spherical[1], 2)) *
222 this->hamiltonian_is_diagonal =
false;
223 sort_by_quantum_number_f =
false;
225 if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
226 numerical_precision &&
227 std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
228 numerical_precision) {
229 *this->hamiltonian -=
static_cast<real_t>(2 / 12.) * magnetic_field_spherical[0] *
230 magnetic_field_spherical[2] *
232 *this->hamiltonian -=
static_cast<real_t>(1 / 12.) * magnetic_field_spherical[0] *
233 magnetic_field_spherical[2] *
235 this->hamiltonian_is_diagonal =
false;
236 sort_by_quantum_number_f =
false;
238 if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
239 numerical_precision &&
240 std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
241 numerical_precision) {
242 *this->hamiltonian -=
static_cast<real_t>(std::sqrt(3.0) / 12.) *
243 magnetic_field_spherical[1] * magnetic_field_spherical[2] *
245 this->hamiltonian_is_diagonal =
false;
246 sort_by_quantum_number_f =
false;
247 sort_by_quantum_number_m =
false;
249 if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
250 numerical_precision &&
251 std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
252 numerical_precision) {
253 *this->hamiltonian -=
static_cast<real_t>(std::sqrt(3.0) / 12.) *
254 magnetic_field_spherical[1] * magnetic_field_spherical[0] *
256 this->hamiltonian_is_diagonal =
false;
257 sort_by_quantum_number_f =
false;
258 sort_by_quantum_number_m =
false;
260 if (std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
261 numerical_precision) {
262 *this->hamiltonian -=
static_cast<real_t>(std::sqrt(1.5) / 12.) *
263 static_cast<Scalar
>(std::pow(magnetic_field_spherical[2], 2)) *
265 this->hamiltonian_is_diagonal =
false;
266 sort_by_quantum_number_f =
false;
267 sort_by_quantum_number_m =
false;
269 if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
270 numerical_precision) {
271 *this->hamiltonian -=
static_cast<real_t>(std::sqrt(1.5) / 12.) *
272 static_cast<Scalar
>(std::pow(magnetic_field_spherical[0], 2)) *
274 this->hamiltonian_is_diagonal =
false;
275 sort_by_quantum_number_f =
false;
276 sort_by_quantum_number_m =
false;
281 this->blockdiagonalizing_labels.clear();
282 if (sort_by_quantum_number_f) {
285 if (sort_by_quantum_number_m) {
288 if (sort_by_parity) {
294template class SystemAtom<double>;
295template class SystemAtom<std::complex<double>>;
Type & set_ion_charge(real_t charge)
typename traits::CrtpTraits< Type >::real_t real_t
Type & set_diamagnetism_enabled(bool enable)
SystemAtom(std::shared_ptr< const basis_t > basis)
Type & set_electric_field(const std::array< real_t, 3 > &field)
Type & set_magnetic_field(const std::array< real_t, 3 > &field)
Type & set_ion_distance_vector(const std::array< real_t, 3 > &vector)
Type & set_ion_interaction_order(int value)
std::shared_ptr< const basis_t > get_basis() const
Matrix< Type, Size, 1 > Vector
Matrix< Type, 3, 1 > Vector3
Matrix< Type, Dynamic, 1 > VectorX
@ ELECTRIC_QUADRUPOLE_ZERO
@ SORT_BY_QUANTUM_NUMBER_M
@ SORT_BY_QUANTUM_NUMBER_F
Helper struct to extract types from a numerical type.