Line data Source code
1 : # SPDX-FileCopyrightText: 2025 PairInteraction Developers
2 : # SPDX-License-Identifier: LGPL-3.0-or-later
3 1 : from __future__ import annotations
4 :
5 1 : import logging
6 1 : from collections.abc import Iterable
7 1 : from typing import TYPE_CHECKING, overload
8 :
9 1 : import numpy as np
10 1 : from scipy import sparse
11 1 : from scipy.sparse import csr_matrix
12 1 : from typing_extensions import deprecated
13 :
14 1 : import pairinteraction as pi_complex
15 1 : import pairinteraction.real as pi_real
16 1 : from pairinteraction.perturbative.perturbation_theory import calculate_perturbative_hamiltonian
17 1 : from pairinteraction.units import QuantityArray, QuantityScalar, ureg
18 :
19 : if TYPE_CHECKING:
20 : from collections.abc import Sequence
21 :
22 : from scipy.sparse import csr_matrix
23 :
24 : from pairinteraction.ket import KetAtomTuple, KetPairLike
25 : from pairinteraction.units import NDArray, PintArray, PintFloat
26 :
27 : SystemPair = pi_real.SystemPair | pi_complex.SystemPair
28 :
29 1 : logger = logging.getLogger(__name__)
30 :
31 :
32 : @overload
33 : def get_effective_hamiltonian_from_system(
34 : ket_tuple_list: Sequence[KetPairLike],
35 : system_pair: SystemPair,
36 : order: int = 2,
37 : required_overlap: float = 0.9,
38 : *,
39 : return_only_specified_order: bool = False,
40 : unit: None = None,
41 : ) -> tuple[PintArray, csr_matrix]: ...
42 :
43 :
44 : @overload
45 : def get_effective_hamiltonian_from_system(
46 : ket_tuple_list: Sequence[KetPairLike],
47 : system_pair: SystemPair,
48 : order: int = 2,
49 : required_overlap: float = 0.9,
50 : return_only_specified_order: bool = False,
51 : *,
52 : unit: str,
53 : ) -> tuple[NDArray, csr_matrix]: ...
54 :
55 :
56 1 : @deprecated("Use EffectiveSystemPair(ket_tuples).get_effective_hamiltonian() instead.")
57 1 : def get_effective_hamiltonian_from_system(
58 : ket_tuple_list: Sequence[KetPairLike],
59 : system_pair: SystemPair,
60 : order: int = 2,
61 : required_overlap: float = 0.9,
62 : return_only_specified_order: bool = False,
63 : unit: str | None = None,
64 : ) -> tuple[NDArray | PintArray, csr_matrix]:
65 : r"""Get the perturbative Hamiltonian at a desired order in Rayleigh-Schrödinger perturbation theory.
66 :
67 : This function takes a list of tuples of ket states, which forms the basis of the model space in which the effective
68 : Hamiltonian is calculated. The whole Hamiltonian is taken from a pair system.
69 : The Hamiltonian of the pair system is assumed to be diagonal in the unperturbed Hamiltonian, all off-diagonal
70 : elements are assumed to belong to the perturbative term.
71 : The function also checks for resonances between all states and states in the model space.
72 :
73 : Args:
74 : ket_tuple_list: List of all pair states that span up the model space.
75 : The effective Hamiltonian is calculated for these states.
76 : system_pair: Two-Atom-System, diagonal in the basis of the unperturbed Hamiltonian.
77 : order: Order up to which the perturbation theory is expanded. Support up to third order.
78 : Default is second order.
79 : required_overlap: If set, the code checks for validity of a perturbative treatment.
80 : Error is thrown if the perturbed eigenstate has less overlap than this value with the
81 : unperturbed eigenstate.
82 : return_only_specified_order: If True, the returned effective Hamiltonian will only contain the specified order.
83 : Default is False, which returns the sum of all orders up to the specified order.
84 : unit: The unit to which to convert the result. Default None will return a pint quantity.
85 :
86 : Returns:
87 : - Effective Hamiltonian as a :math:`m \times m` matrix, where m is the length of `ket_tuple_list`.
88 : - Eigenvectors in perturbation theory due to interaction with states out of the model space, returned as
89 : a sparse matrix in compressed row format. Each row represents the corresponding eigenvector.
90 :
91 : Raises:
92 : ValueError: If a resonance between a state in the model space and a state not in the model space occurs.
93 :
94 : """
95 1 : if np.isinf(system_pair.get_distance().magnitude):
96 0 : raise ValueError(
97 : "Pair system is initialized without a distance. "
98 : "Please set a distance for calculating an effective Hamiltonian."
99 : )
100 :
101 1 : model_inds = _get_model_inds(ket_tuple_list, system_pair)
102 1 : h_au = system_pair.get_hamiltonian().to_base_units().magnitude # Hamiltonian in atomic units
103 1 : h_eff_dict_au, eigvec_perturb = calculate_perturbative_hamiltonian(h_au, model_inds, order)
104 1 : if not 0 <= required_overlap <= 1:
105 0 : raise ValueError("Required overlap has to be a positive real number between zero and one.")
106 1 : if required_overlap > 0:
107 1 : _check_for_resonances(model_inds, eigvec_perturb, system_pair, required_overlap)
108 :
109 1 : if return_only_specified_order:
110 1 : h_eff_au: NDArray = h_eff_dict_au[order]
111 : else:
112 1 : h_eff_au = sum(h_eff for h_eff in h_eff_dict_au.values()) # type: ignore [assignment]
113 :
114 1 : h_eff = QuantityArray.convert_au_to_user(h_eff_au, "energy", unit)
115 1 : return h_eff, eigvec_perturb
116 :
117 :
118 : @overload
119 : def get_c3_from_system(
120 : ket_tuple_list: Sequence[KetPairLike], system_pair: SystemPair, *, unit: None = None
121 : ) -> PintFloat: ...
122 :
123 :
124 : @overload
125 : def get_c3_from_system(ket_tuple_list: Sequence[KetPairLike], system_pair: SystemPair, unit: str) -> float: ...
126 :
127 :
128 1 : @deprecated("Use C3(ket1, ket2).get() instead.")
129 1 : def get_c3_from_system(
130 : ket_tuple_list: Sequence[KetPairLike], system_pair: SystemPair, unit: str | None = None
131 : ) -> float | PintFloat:
132 : r"""Calculate the :math:`C_3` coefficient for a list of two 2-tuples of single-atom ket states.
133 :
134 : This function calculates the :math:`C_3` coefficient in the desired unit. The input is a list of two 2-tuples of
135 : single-atom ket states. We use the convention :math:`\Delta E = \frac{C_3}{r^3}`.
136 :
137 : Args:
138 : ket_tuple_list: The input as a list of tuples of two states [(a,b),(c,d)],
139 : the :math:`C_3` coefficient is calculated for (a,b)->(c,d).
140 : If there are not exactly two tuples in the list, a ValueError is raised.
141 : system_pair: The pair system that is used for the calculation.
142 : unit: The unit to which to convert the result. Default None will return a pint quantity.
143 :
144 : Returns:
145 : The :math:`C_3` coefficient with its unit.
146 :
147 : Raises:
148 : ValueError: If a list of not exactly two tuples of single-atom states is given.
149 :
150 : """
151 1 : if len(ket_tuple_list) != 2:
152 0 : raise ValueError("C3 coefficient can be calculated only between two 2-atom states.")
153 :
154 1 : r = system_pair.get_distance()
155 1 : if np.isinf(r.magnitude):
156 0 : logger.warning(
157 : "Pair system is initialized without a distance. "
158 : "Calculating the C3 coefficient at a distance vector of [0, 0, 20] mum."
159 : )
160 0 : old_distance_vector = system_pair.get_distance_vector()
161 0 : system_pair.set_distance_vector([0, 0, 20], "micrometer")
162 0 : c3 = get_c3_from_system(ket_tuple_list, system_pair, unit=unit)
163 0 : system_pair.set_distance_vector(old_distance_vector)
164 0 : return c3
165 :
166 1 : h_eff, _ = get_effective_hamiltonian_from_system(ket_tuple_list, system_pair, order=1)
167 1 : c3_pint = h_eff[0, 1] * r**3 # type: ignore [index] # PintArray does not know it can be indexed
168 1 : return QuantityScalar.from_pint(c3_pint, "c3").to_pint_or_unit(unit)
169 :
170 :
171 : @overload
172 : def get_c6_from_system(ket_tuple: KetPairLike, system_pair: SystemPair, *, unit: None = None) -> PintFloat: ...
173 :
174 :
175 : @overload
176 : def get_c6_from_system(ket_tuple: KetPairLike, system_pair: SystemPair, unit: str) -> float: ...
177 :
178 :
179 1 : @deprecated("Use C6(ket1, ket2).get() instead.")
180 1 : def get_c6_from_system(ket_tuple: KetPairLike, system_pair: SystemPair, unit: str | None = None) -> float | PintFloat:
181 : r"""Calculate the :math:`C_6` coefficient for a given tuple of ket states.
182 :
183 : This function calculates the :math:`C_6` coefficient in the desired unit. The input is a 2-tuple of single-atom ket
184 : states.
185 :
186 : Args:
187 : ket_tuple: The input is a tuple repeating the same single-atom state in the format (a,a).
188 : If a tuple with not exactly two identical states is given, a ValueError is raised.
189 : system_pair: The pair system that is used for the calculation.
190 : unit: The unit to which to convert the result. Default None will return a pint quantity.
191 :
192 : Returns:
193 : The :math:`C_6` coefficient. If a unit is specified, the value in this unit is returned.
194 :
195 : Raises:
196 : ValueError: If a tuple with more than two single-atom states is given.
197 :
198 : """
199 1 : if isinstance(ket_tuple, Iterable):
200 1 : if len(ket_tuple) != 2:
201 0 : raise ValueError("C6 coefficient can be calculated only for a single 2-atom state.")
202 1 : if ket_tuple[0].species == ket_tuple[1].species and ket_tuple[0] != ket_tuple[1]:
203 0 : raise ValueError(
204 : "If you want to calculate 2nd order perturbations of two different states a and b, "
205 : "please use the get_effective_hamiltonian_from_system([(a,b), (b,a)], system_pair) function."
206 : )
207 :
208 1 : r = system_pair.get_distance()
209 1 : if np.isinf(r.magnitude):
210 0 : logger.warning(
211 : "Pair system is initialized without a distance. "
212 : "Calculating the C6 coefficient at a distance vector of [0, 0, 20] mum."
213 : )
214 0 : old_distance_vector = system_pair.get_distance_vector()
215 0 : system_pair.set_distance_vector([0, 0, 20], "micrometer")
216 0 : c6 = get_c6_from_system(ket_tuple, system_pair, unit=unit)
217 0 : system_pair.set_distance_vector(old_distance_vector)
218 0 : return c6
219 :
220 1 : h_eff, _ = get_effective_hamiltonian_from_system(
221 : [ket_tuple], system_pair, order=2, return_only_specified_order=True
222 : )
223 1 : c6_pint = h_eff[0, 0] * r**6 # type: ignore [index] # PintArray does not know it can be indexed
224 1 : return QuantityScalar.from_pint(c6_pint, "c6").to_pint_or_unit(unit)
225 :
226 :
227 1 : def _get_model_inds(ket_tuple_list: Sequence[KetPairLike], system_pair: SystemPair) -> list[int]:
228 : """Get the indices of all ket tuples in the basis of pair system.
229 :
230 : This function takes a list of 2-tuples of ket states, and a pair system holding the entire basis.
231 : It returns an array of indices of the pair system basis in the order of the tuple list.
232 :
233 : Args:
234 : ket_tuple_list: List of all pair states that span up the model space.
235 : system_pair: Two-Atom-System, diagonal in the basis of the unperturbed Hamiltonian.
236 :
237 : Returns:
238 : List of indices corresponding to the states that span up the model space.
239 :
240 : """
241 1 : model_inds = []
242 1 : for kets in ket_tuple_list:
243 1 : overlap = system_pair.basis.get_overlaps(kets)
244 1 : index = np.argmax(overlap)
245 1 : if overlap[index] == 0:
246 0 : raise ValueError(f"The pairstate {kets} is not part of the basis of the pair system.")
247 1 : if overlap[index] < 0.5:
248 0 : raise ValueError(f"The pairstate {kets} cannot be identified uniquely (max overlap: {overlap[index]}).")
249 1 : model_inds.append(int(index))
250 1 : return model_inds
251 :
252 :
253 1 : def _check_for_resonances(
254 : model_inds: list[int],
255 : eigvec_perturb: csr_matrix,
256 : system_pair: SystemPair,
257 : required_overlap: float,
258 : ) -> None:
259 : r"""Check for resonance between the states in the model space and other states.
260 :
261 : This function takes the perturbed eigenvectors of the perturbation theory as an input.
262 : If the overlap of the perturbed eigenstate with its corresponding unperturbed state are too small,
263 : this function raises an error, as perturbation theory breaks down.
264 : In this case, it also prints all states with a relevant admixture that should therefore be also included in the
265 : model space, to allow perturbation theory.
266 :
267 : Args:
268 : model_inds: List of indices corresponding to the states that span up the model space.
269 : eigvec_perturb: Sparse representation of the perturbed eigenstates in the desired order of
270 : perturbation theory. Each row corresponds to the eigestate according to `state model indices.`
271 : system_pair: Two-Atom-System, diagonal in the basis of the unperturbed Hamiltonian.
272 : order: Order up to which the perturbation theory is expanded. Support up to third order.
273 : Default is second order.
274 : required_overlap: If set, the code checks for validity of a perturbative treatment.
275 : Error is thrown if the perturbed eigenstate has less overlap than this value with the unperturbed eigenstate
276 :
277 : Returns:
278 : Effective Hamiltonian as a :math:`m \times m` matrix, where m is the length of `ket_tuple_list`
279 : Eigenvectors in perturbation theory due to interaction with states out of the model space,
280 : returned as a sparse matrix in compressed row format. Each row represent the corresponding eigenvector
281 :
282 : Raises:
283 : ValueError: If a resonance between a state in the model space and a state not in the model space occurs.
284 :
285 : """
286 1 : overlaps = (eigvec_perturb.multiply(eigvec_perturb.conj())).real
287 1 : error_flag = False
288 1 : for i, j in zip(range(len(model_inds)), model_inds, strict=True):
289 1 : vector_norm = sparse.linalg.norm(overlaps[i, :])
290 1 : overlap = overlaps[i, j] / vector_norm
291 1 : if overlap >= required_overlap:
292 1 : continue
293 0 : error_flag = True
294 0 : print_above_admixture = (1 - required_overlap) * 0.05
295 0 : indices = sparse.find(overlaps[i, :] >= print_above_admixture * vector_norm)[1]
296 0 : logger.error(
297 : "The state %s has resonances with the following states, please consider adding them to your model space:",
298 : system_pair.basis.get_ket(j),
299 : )
300 0 : for index in indices:
301 0 : if index == j:
302 0 : continue
303 0 : admixture = 1 if np.isinf(overlaps[i, index]) else overlaps[i, index] / vector_norm
304 0 : logger.error(" - %s with admixture %.3f", system_pair.basis.get_ket(int(index)), admixture)
305 1 : if error_flag:
306 0 : raise ValueError(
307 : "Error. Perturbative Calculation not possible due to resonances. "
308 : "Add more states to the model space or adapt your required overlap."
309 : )
310 :
311 :
312 1 : @deprecated("Use EffectiveSystemPair(ket_tuples) instead to create a system for perturbative calculations.")
313 1 : def create_system_for_perturbative( # noqa: C901, PLR0912, PLR0915
314 : ket_tuple_list: Sequence[KetAtomTuple],
315 : electric_field: PintArray | None = None,
316 : magnetic_field: PintArray | None = None,
317 : distance_vector: PintArray | None = None,
318 : multipole_order: int = 3,
319 : with_diamagnetism: bool = False,
320 : perturbation_order: int = 2,
321 : number_of_considered_pair_kets: int = 2_000,
322 : ) -> SystemPair:
323 : r"""Create a good estimate for a system in which to perform perturbative calculations.
324 :
325 : This function takes a list of 2-tuples of ket states and creates a pair system holding a larger basis.
326 : The parameters of the basis are adjusted by the electric and magnetic field vectors of the system, as well
327 : as the distance and the multipole-order of the interaction. For higher-order perturbation theory, larger
328 : systems can be considered. Diamagnetism can be considered as well.
329 :
330 : Args:
331 : ket_tuple_list: List of all pair states that span up the model space. The system is created such that
332 : the effective Hamiltonian of the model system can be calculated accurately at a later stage.
333 : electric_field: Electric field in the system.
334 : magnetic_field: Magnetic field in the system.
335 : distance_vector: Distance vector between the atoms.
336 : multipole_order: Multipole-order of the interaction. Default is 3 (dipole-dipole).
337 : with_diamagnetism: True if diamagnetic term should be considered. Default is False.
338 : perturbation_order: Order of perturbative calculation the system shall be used for. Default is 2.
339 : number_of_considered_pair_kets: Number of pair kets that are considered in the system. Default is 2000.
340 :
341 : Returns:
342 : Pair system that can be used for perturbative calculations.
343 :
344 : """
345 1 : electric_field = electric_field if electric_field is not None else ureg.Quantity([0, 0, 0], "V/cm")
346 1 : magnetic_field = magnetic_field if magnetic_field is not None else ureg.Quantity([0, 0, 0], "G")
347 :
348 1 : is_real = isinstance(ket_tuple_list[0][0], pi_real.KetAtom)
349 1 : pi = pi_real if is_real else pi_complex
350 1 : are_fields_along_z = all(x == 0 for x in [*magnetic_field[:2], *electric_field[:2]]) # type: ignore [index]
351 :
352 1 : system_atoms: list[pi_real.SystemAtom | pi_complex.SystemAtom] = []
353 :
354 1 : delta_n = 7
355 1 : delta_l = perturbation_order * (multipole_order - 2)
356 1 : for i in range(2):
357 1 : kets = [ket_tuple[i] for ket_tuple in ket_tuple_list]
358 1 : nlfm = np.transpose([[ket.n, ket.l, ket.f, ket.m] for ket in kets])
359 1 : n_range = (int(np.min(nlfm[0])) - delta_n, int(np.max(nlfm[0])) + delta_n)
360 1 : l_range = (np.min(nlfm[1]) - delta_l, np.max(nlfm[1]) + delta_l)
361 1 : if any(ket.is_calculated_with_mqdt for ket in kets):
362 : # for mqdt we increase delta_l by 1 to take into account the variance ...
363 0 : l_range = (np.min(nlfm[1]) - delta_l - 1, np.max(nlfm[1]) + delta_l + 1)
364 1 : m_range = None
365 1 : if are_fields_along_z:
366 1 : m_range = (np.min(nlfm[3]) - delta_l, np.max(nlfm[3]) + delta_l)
367 1 : basis = pi.BasisAtom(kets[0].species, n=n_range, l=l_range, m=m_range)
368 1 : system = pi.SystemAtom(basis)
369 1 : system.set_diamagnetism_enabled(with_diamagnetism)
370 1 : system.set_magnetic_field(magnetic_field)
371 1 : system.set_electric_field(electric_field)
372 1 : system_atoms.append(system)
373 :
374 1 : pi.diagonalize(system_atoms)
375 :
376 1 : pair_energies_au = [
377 : sum(
378 : system.get_corresponding_energy(ket).to_base_units().magnitude
379 : for system, ket in zip(system_atoms, ket_tuple, strict=True)
380 : )
381 : for ket_tuple in ket_tuple_list
382 : ]
383 :
384 1 : def get_basis_pair(delta_energy_au: float) -> pi_real.BasisPair | pi_complex.BasisPair:
385 1 : return pi.BasisPair( # type: ignore [no-any-return]
386 : system_atoms,
387 : energy=(min(pair_energies_au) - delta_energy_au, max(pair_energies_au) + delta_energy_au),
388 : energy_unit="hartree",
389 : )
390 :
391 1 : mhz_au = QuantityScalar.convert_user_to_au(1, "MHz", "energy")
392 1 : delta_energy_au = mhz_au
393 1 : min_delta, max_delta = None, None
394 :
395 : # make a bisect search to get a sensible basis size between:
396 : # number_of_considered_pair_kets and 1.2 * number_of_considered_pair_kets
397 1 : while delta_energy_au < 1: # stop if delta_energy_au is 1 and the basis is still very small
398 1 : basis_pair = get_basis_pair(delta_energy_au)
399 1 : if basis_pair.number_of_kets < number_of_considered_pair_kets:
400 1 : min_delta = delta_energy_au
401 1 : if max_delta is None:
402 1 : delta_energy_au *= 2
403 : else:
404 0 : delta_energy_au = (delta_energy_au + max_delta) / 2
405 1 : elif basis_pair.number_of_kets > number_of_considered_pair_kets * 1.2:
406 1 : max_delta = delta_energy_au
407 1 : if min_delta is None:
408 0 : delta_energy_au /= 2
409 : else:
410 1 : delta_energy_au = (delta_energy_au + min_delta) / 2
411 : else:
412 1 : break
413 1 : if max_delta is not None and min_delta is not None and max_delta - min_delta < mhz_au:
414 0 : break
415 :
416 1 : logger.debug("The pair basis for the perturbative calculations consists of %d kets.", basis_pair.number_of_kets)
417 :
418 1 : system_pair = pi.SystemPair(basis_pair)
419 1 : if distance_vector is not None:
420 1 : system_pair.set_distance_vector(distance_vector)
421 1 : system_pair.set_interaction_order(multipole_order)
422 1 : return system_pair # type: ignore [no-any-return]
|
