Coverage for gpaw/response/mft.py: 99%
218 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1from __future__ import annotations
3from abc import abstractmethod
4import warnings
6import numpy as np
8from gpaw.typing import Vector
9from gpaw.response import ResponseGroundStateAdaptable, ResponseContextInput
10from gpaw.response.frequencies import ComplexFrequencyDescriptor
11from gpaw.response.chiks import (ChiKSCalculator, RealAxisWarning,
12 get_smat_components, smat,
13 regularize_intraband_transitions)
14from gpaw.response.localft import LocalFTCalculator, add_LSDA_Wxc
15from gpaw.response.site_kernels import SiteKernels
16from gpaw.response.site_data import AtomicSites
17from gpaw.response.pair_integrator import PairFunction, PairFunctionIntegrator
18from gpaw.response.pair_transitions import PairTransitions
19from gpaw.response.matrix_elements import (SitePairDensityCalculator,
20 SiteZeemanPairEnergyCalculator,
21 SiteSpinPairEnergyCalculator)
23from ase.units import Hartree
26class IsotropicExchangeCalculator:
27 r"""Calculator class for the Heisenberg exchange constants
29 _ 2
30 J^ab(q) = - ‾‾ B^(xc†) K^(a†)(q) χ_KS^('+-)(q) K^b(q) B^(xc) (1)
31 V0
33 calculated for an isotropic system in a plane wave representation using
34 the magnetic force theorem within second order perturbation theory, see
35 [J. Phys.: Condens. Matter 35 (2023) 105802].
37 Entering the formula for the isotropic exchange constant at wave vector q
38 between sublattice a and b is the unit cell volume V0, the functional
39 derivative of the (LDA) exchange-correlation energy with respect to the
40 magnitude of the magnetization B^(xc), the sublattice site kernels K^a(q)
41 and K^b(q) as well as the reactive part of the static transverse magnetic
42 susceptibility of the Kohn-Sham system χ_KS^('+-)(q).
44 NB: To achieve numerical stability of the plane-wave implementation, we
45 use instead the following expression to calculate exchange parameters:
47 ˷ 2
48 J^ab(q) = - ‾‾ W_xc^(z†) K^(a†)(q) χ_KS^('+-)(q) K^b(q) W_xc^z (2)
49 V0
51 We do this since B^(xc)(r) = -|W_xc^z(r)| is nonanalytic in points of space
52 where the spin-polarization changes sign, why it is problematic to evaluate
53 Eq. (1) numerically within a plane-wave representation.
54 If the site partitionings only include spin-polarization of the same sign,
55 Eqs. (1) and (2) should yield identical exchange parameters, but for
56 antiferromagnetically aligned sites, the coupling constants differ by a
57 sign.
59 The site kernels encode the partitioning of real space into sites of the
60 Heisenberg model. This is not a uniquely defined procedure, why the user
61 has to define them externally through the SiteKernels interface."""
63 def __init__(self,
64 chiks_calc: ChiKSCalculator,
65 localft_calc: LocalFTCalculator):
66 """Construct the IsotropicExchangeCalculator object."""
67 # Check that chiks has the assumed properties
68 assumed_props = dict(
69 gammacentered=True,
70 nblocks=1
71 )
72 for key, item in assumed_props.items():
73 assert getattr(chiks_calc, key) == item, \
74 f'Expected chiks.{key} == {item}. '\
75 f'Got: {getattr(chiks_calc, key)}'
77 self.chiks_calc = chiks_calc
78 self.context = chiks_calc.context
80 # Check assumed properties of the LocalFTCalculator
81 assert localft_calc.context is self.context
82 assert localft_calc.gs is chiks_calc.gs
83 self.localft_calc = localft_calc
85 # W_xc^z buffer
86 self._Wxc_G = None
88 # χ_KS^('+-) buffer
89 self._chiksr = None
91 def __call__(self, q_c, site_kernels: SiteKernels, txt=None):
92 """Calculate the isotropic exchange constants for a given wavevector.
94 Parameters
95 ----------
96 q_c : nd.array
97 Wave vector q in relative coordinates
98 site_kernels : SiteKernels
99 Site kernels instance defining the magnetic sites of the crystal
100 txt : str
101 Separate file to store the chiks calculation output in (optional).
102 If not supplied, the output will be written to the standard text
103 output location specified when initializing chiks.
105 Returns
106 -------
107 J_abp : nd.array (dtype=complex)
108 Isotropic Heisenberg exchange constants between magnetic sites a
109 and b for all the site partitions p given by the site_kernels.
110 """
111 # Get ingredients
112 Wxc_G = self.get_Wxc()
113 chiksr = self.get_chiksr(q_c, txt=txt)
114 qpd, chiksr_GG = chiksr.qpd, chiksr.array[0] # array = chiksr_zGG
115 V0 = qpd.gd.volume
117 # Allocate an array for the exchange constants
118 nsites = site_kernels.nsites
119 J_pab = np.empty(site_kernels.shape + (nsites,), dtype=complex)
121 # Compute exchange coupling
122 for J_ab, K_aGG in zip(J_pab, site_kernels.calculate(qpd)):
123 for a in range(nsites):
124 for b in range(nsites):
125 J = np.conj(Wxc_G) @ np.conj(K_aGG[a]).T @ chiksr_GG \
126 @ K_aGG[b] @ Wxc_G
127 J_ab[a, b] = - 2. * J / V0
129 # Transpose to have the partitions index last
130 J_abp = np.transpose(J_pab, (1, 2, 0))
132 return J_abp * Hartree # Convert from Hartree to eV
134 def get_Wxc(self):
135 """Get B^(xc)_G from buffer."""
136 if self._Wxc_G is None: # Calculate if buffer is empty
137 self._Wxc_G = self._calculate_Wxc()
139 return self._Wxc_G
141 def _calculate_Wxc(self):
142 """Calculate the Fourier transform W_xc^z(G)."""
143 # Create a plane wave descriptor encoding the plane wave basis. Input
144 # q_c is arbitrary, since we are assuming that chiks.gammacentered == 1
145 qpd0 = self.chiks_calc.get_pw_descriptor([0., 0., 0.])
147 return self.localft_calc(qpd0, add_LSDA_Wxc)
149 def get_chiksr(self, q_c, txt=None):
150 """Get χ_KS^('+-)(q) from buffer."""
151 q_c = np.asarray(q_c)
153 # Calculate if buffer is empty or a new q-point is given
154 if self._chiksr is None or not np.allclose(q_c, self._chiksr.q_c):
155 self._chiksr = self._calculate_chiksr(q_c, txt=txt)
157 return self._chiksr
159 def _calculate_chiksr(self, q_c, txt=None):
160 r"""Use the ChiKSCalculator to calculate the reactive part of the
161 static Kohn-Sham susceptibility χ_KS^('+-)(q).
163 First, the dynamic Kohn-Sham susceptibility
165 __ __
166 1 \ \ f_nk↑ - f_mk+q↓
167 χ_KS,GG'^+-(q,ω+iη) = ‾ / / ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
168 V ‾‾ ‾‾ ħω - (ε_mk+q↓ - ε_nk↑) + iħη
169 k n,m
170 x n_nk↑,mk+q↓(G+q) n_mk+q↓,nk↑(-G'-q)
172 is calculated in the static limit ω=0 and without broadening η=0. Then,
173 the reactive part (see [PRB 103, 245110 (2021)]) is extracted:
175 1
176 χ_KS,GG'^(+-')(q,z) = ‾ [χ_KS,GG'^+-(q,z) + χ_KS,-G'-G^-+(-q,-z*)].
177 2
178 """
179 # Initiate new output file, if supplied
180 if txt is not None:
181 self.context.new_txt_and_timer(txt)
183 # Even though the Heisenberg exchange constants are difficult to
184 # converge for metals, it does not really help to add finite broadening
185 # of the susceptibility. Therefore, we bite the sour apple and always
186 # evaluate the χ_KS on the real axis.
187 zd = ComplexFrequencyDescriptor.from_array([0. + 0.j])
188 with warnings.catch_warnings():
189 warnings.simplefilter('ignore', category=RealAxisWarning)
190 chiks = self.chiks_calc.calculate('+-', q_c, zd)
191 if np.allclose(q_c, 0.):
192 chiks.symmetrize_reciprocity()
194 # Take the reactive part
195 chiksr = chiks.copy_reactive_part()
197 return chiksr
200def calculate_single_particle_site_magnetization(
201 gs: ResponseGroundStateAdaptable,
202 sites: AtomicSites,
203 context: ResponseContextInput = '-'):
204 """Calculate the single-particle site magnetization.
206 Returns
207 -------
208 sp_magmom_ap : np.ndarray
209 Magnetic moment in μB of site a under partitioning p, calculated based
210 on a single-particle sum rule.
211 """
212 single_particle_calc = SingleParticleSiteMagnetizationCalculator(
213 gs, sites, context=context)
214 site_magnetization = single_particle_calc()
215 return site_magnetization.array
218def calculate_single_particle_site_zeeman_energy(
219 gs: ResponseGroundStateAdaptable,
220 sites: AtomicSites,
221 context: ResponseContextInput = '-'):
222 """Calculate the single-particle site Zeeman energy.
224 Returns
225 -------
226 sp_EZ_ap : np.ndarray
227 Local Zeeman energy in eV of site a under partitioning p, calculated
228 based on a single-particle sum rule.
229 """
230 single_particle_calc = SingleParticleSiteZeemanEnergyCalculator(
231 gs, sites, context=context)
232 site_zeeman_energy = single_particle_calc()
233 return site_zeeman_energy.array * Hartree # Ha -> eV
236def calculate_pair_site_magnetization(
237 gs: ResponseGroundStateAdaptable,
238 sites: AtomicSites,
239 context: ResponseContextInput = '-',
240 q_c=[0., 0., 0.],
241 nbands: int | None = None,
242 nblocks: int | str = 1):
243 """Calculate the pair site magnetization.
245 Parameters
246 ----------
247 q_c : Vector
248 q-vector to evaluate the pair site magnetization for.
249 nbands : int or None
250 Number of bands to include in the band summation of the pair site
251 magnetization. If nbands is None, it includes all bands.
252 nblocks : int or str
253 The workload is parallelized over k-points and band+spin transitions.
254 The latter is divided into nblocks, integrating nprocessors / nblocks
255 k-points at a time.
257 Returns
258 -------
259 magmom_abp : np.ndarray
260 Pair magnetization in μB of site a and b under partitioning p,
261 calculated based on a two-particle sum rule.
262 """
263 two_particle_calc = TwoParticleSiteMagnetizationCalculator(
264 gs, sites, context=context, nbands=nbands, nblocks=nblocks)
265 pair_site_magnetization = two_particle_calc(q_c)
266 return pair_site_magnetization.array
269def calculate_pair_site_zeeman_energy(
270 gs: ResponseGroundStateAdaptable,
271 sites: AtomicSites,
272 context: ResponseContextInput = '-',
273 q_c=[0., 0., 0.],
274 nbands: int | None = None,
275 nblocks: int | str = 1):
276 """Calculate the pair site Zeeman energy.
278 Parameters
279 ----------
280 q_c : Vector
281 q-vector to evaluate the pair site Zeeman energy for.
282 nbands : int or None
283 Number of bands to include in the band summation of the pair site
284 Zeeman energy. If nbands is None, it includes all bands.
285 nblocks : int or str
286 The workload is parallelized over k-points and band+spin transitions.
287 The latter is divided into nblocks, integrating nprocessors / nblocks
288 k-points at a time.
290 Returns
291 -------
292 EZ_abp : np.ndarray
293 Local pair Zeeman energy in eV of site a and b under partitioning p,
294 calculated based on a two-particle sum rule.
295 """
296 two_particle_calc = TwoParticleSiteZeemanEnergyCalculator(
297 gs, sites, context=context, nbands=nbands, nblocks=nblocks)
298 pair_site_zeeman_energy = two_particle_calc(q_c)
299 return pair_site_zeeman_energy.array * Hartree # Ha -> eV
302def calculate_exchange_parameters(
303 gs: ResponseGroundStateAdaptable,
304 sites: AtomicSites,
305 q_c: Vector,
306 context: ResponseContextInput = '-',
307 nbands: int | None = None,
308 nblocks: int | str = 1):
309 """Calculate the Heisenberg exchange parameters.
311 Parameters
312 ----------
313 q_c : Vector
314 q-vector to evaluate the pair site Zeeman energy for.
315 nbands : int or None
316 Number of bands to include in the band summation of the pair site
317 Zeeman energy. If nbands is None, it includes all bands.
318 nblocks : int or str
319 The workload is parallelized over k-points and band+spin transitions.
320 The latter is divided into nblocks, integrating nprocessors / nblocks
321 k-points at a time.
323 Returns
324 -------
325 J_abp : np.ndarray
326 Heisenberg exchange parameters in eV of sites a and b under
327 partitioning p.
328 """
329 heisenberg_calc = HeisenbergExchangeCalculator(
330 gs, sites, context=context, nbands=nbands, nblocks=nblocks)
331 heisenberg_exchange = heisenberg_calc(q_c)
332 return heisenberg_exchange.array
335class SiteFunction(PairFunction):
336 r"""Data object for single-particle site functions f_a.
338 A single-particle site function is understood as any function that can be
339 constructed as a sum over the system eigenstates
340 __
341 \ a
342 f_a = / f
343 ‾‾ α
344 α
346 with site dependent weights f^a_α representing some projection onto a local
347 (atomic) site.
348 """
349 def __init__(self, sites: AtomicSites):
350 self.sites = sites
351 super().__init__(q_c=[0., 0., 0.]) # no crystal momentum transfer
353 @property
354 def shape(self):
355 return self.sites.shape
357 def zeros(self):
358 return np.zeros(self.shape, dtype=float)
361class SingleParticleSiteSumRuleCalculator(PairFunctionIntegrator):
362 r"""Calculator for single-particle site sum rules.
364 For any site matrix element f^a_(nks,n'k's') of the Kohn-Sham system, one
365 may define a single-particle site sum rule by its weighted trace
366 __ __
367 1 \ \
368 f_a^μ = ‾‾‾ / / σ^μ_ss f_nks f^a_(nks,nks)
369 N_k ‾‾ ‾‾
370 k n,s
372 where μ∊{0,z}.
373 """
375 def __init__(self, gs, sites, context='-'):
376 super().__init__(gs, context, qsymmetry=False)
377 self.transitions = self.get_band_and_spin_transitions()
378 # Set up calculator for the f^a matrix element
379 self.sites = sites
380 self.matrix_element_calc = self.create_matrix_element_calculator()
382 @abstractmethod
383 def create_matrix_element_calculator(self):
384 """Create the desired site matrix element calculator."""
386 @abstractmethod
387 def get_pauli_matrix(self):
388 """Get the desired Pauli matrix σ^μ."""
390 def get_band_and_spin_transitions(self):
391 """Set up all intraband transitions (n,s)->(n,s)."""
392 nocc2 = self.gs.nocc2
393 n_n = list(range(nocc2))
394 n_t = np.array(n_n + n_n)
395 s_t = np.array([0] * nocc2 + [1] * nocc2)
396 return PairTransitions(n1_t=n_t, n2_t=n_t, s1_t=s_t, s2_t=s_t)
398 def __call__(self):
399 site_function = SiteFunction(sites=self.sites)
400 self._integrate(site_function, self.transitions)
401 return site_function
403 def add_integrand(self, kptpair, weight, site_function):
404 r"""Add the integrand of the outer k-point integral.
406 The integrand is given by (see gpaw.response.pair_integrator)
407 __
408 \
409 (...)_k = V0 / σ^μ_ss f_nks f^a_(nks,nks)
410 ‾‾
411 n,s
412 """
413 # Calculate matrix elements
414 site_matrix_element = self.matrix_element_calc(
415 kptpair, site_function.q_c)
416 assert site_matrix_element.tblocks.blockcomm.size == 1
417 f_tap = site_matrix_element.get_global_array()
419 # Since we only use diagonal site matrix elements, corresponding
420 # to the expectation value of the real functions Θ(r∊Ω_ap) and f(r),
421 # f^a_(nks,nks) = <ψ_nks|Θ(r∊Ω_ap)f(r)|ψ_nks>,
422 # the matrix elements are real
423 assert np.allclose(f_tap.imag, 0.)
424 f_tap = f_tap.real
426 # Calculate Pauli matrix factors and multiply the occupations
427 sigma = self.get_pauli_matrix()
428 sigma_t = sigma[kptpair.transitions.s1_t, kptpair.transitions.s2_t]
429 f_t = kptpair.get_all(kptpair.ikpt1.f_myt)
430 sigmaf_t = sigma_t * f_t
432 # Calculate and add integrand
433 site_function.array[:] += self.gs.volume * weight * np.einsum(
434 't, tap -> ap', sigmaf_t, f_tap)
437class SingleParticleSiteMagnetizationCalculator(
438 SingleParticleSiteSumRuleCalculator):
439 r"""Calculator for the single-particle site magnetization sum rule.
441 The site magnetization is calculated from the site pair density:
442 __ __
443 1 \ \
444 n_a^z = ‾‾‾ / / σ^z_ss f_nks n^a_(nks,nks)
445 N_k ‾‾ ‾‾
446 k n,s
447 """
448 def create_matrix_element_calculator(self):
449 return SitePairDensityCalculator(self.gs, self.context, self.sites)
451 def get_pauli_matrix(self):
452 return smat('z')
455class SingleParticleSiteZeemanEnergyCalculator(
456 SingleParticleSiteMagnetizationCalculator):
457 r"""Calculator for the single-particle site Zeeman energy sum rule.
458 __ __
459 1 \ \
460 E_a^Z = ‾‾‾ / / σ^z_ss f_nks E^(Z,a)_(nks,nks)
461 N_k ‾‾ ‾‾
462 k n,s
463 """
464 def create_matrix_element_calculator(self):
465 return SiteZeemanPairEnergyCalculator(
466 self.gs, self.context, self.sites, rshewmin=1e-8)
469class SitePairFunction(PairFunction):
470 r"""Data object for site pair functions.
472 A site pair function is understood as any function that can be written on
473 the form of a pair function,
474 __
475 \ ab
476 pf_ab(q) = / pf δ_{q,q_{α',α}}
477 ‾‾ αα'
478 α,α'
480 with site-dependent pair function weights pf^(ab)_{αα'}.
482 Typically, the site pair function will be related to a more general lattice
483 periodic pair function pf(r,r') = pf(r+R,r'+R), which can be written in
484 terms of its lattice Fourier transform
486 V0 /
487 pf(r,r') = ‾‾‾‾‾‾ | dq pf(r,r',q)
488 (2π)^D /
489 BZ
490 where
491 __
492 \ iq⋅R
493 pf(r,r',q) = / e pf(r,r'+R)
494 ‾‾
495 R
497 The site-projected lattice Fourier transform then constitutes a site pair
498 function:
500 //
501 pf_ab(q) = || drdr' Θ(r∊Ω_a) pf(r,r',q) Θ(r'∊Ω_b)
502 //
503 """
504 def __init__(self, q_c: Vector, sites: AtomicSites):
505 self.sites = sites
506 super().__init__(q_c)
508 @property
509 def shape(self):
510 nsites = len(self.sites)
511 npartitions = self.sites.npartitions
512 return nsites, nsites, npartitions
514 def zeros(self):
515 return np.zeros(self.shape, dtype=complex)
518class SitePairFunctionCalculator(PairFunctionIntegrator):
519 r"""Calculator for site-projected pair functions.
521 In the Kohn-Sham system, site-projected pair functions are constructed
522 straight-forwardly as a sum over Kohn-Sham eigenstate transitions,
523 __ __ __
524 1 \ \ \ /
525 pf_ab(q) = ‾‾‾ / / / | σ^μ_ss' σ^ν_s's w_(ε_nks,ε_n'k+qs')
526 N_k ‾‾ ‾‾ ‾‾ \ \
527 k n,n' s,s' × f^a_(nks,n'k+qs') g^b_(n'k+qs',nks) |
528 /
530 summing up the site-projected matrix elements f^a and f^b, weighted by
531 Pauli-like 2x2 spin-matrices σ^μ and σ^ν and some function
532 w_(ε_nks,ε_n'k+qs') of the Kohn-Sham eigenvalues.
533 """
534 def __init__(self,
535 gs: ResponseGroundStateAdaptable,
536 sites: AtomicSites,
537 context: ResponseContextInput = '-',
538 nbands: int | None = None,
539 nblocks: int | str = 1):
540 """Construct the two-particle site sum rule calculator."""
541 super().__init__(gs, context,
542 # Disable q-symmetry for now. To enable it, we need
543 # to implement site pair function symmetrization.
544 qsymmetry=False,
545 nblocks=nblocks)
546 self.nbands = nbands
547 self.bandsummation = 'double'
548 self.transitions = self.get_band_and_spin_transitions()
550 # Set up calculators for the f^a and g^b matrix elements
551 self.sites = sites
552 mecalc1, mecalc2 = self.create_matrix_element_calculators()
553 self.matrix_element_calc1 = mecalc1
554 self.matrix_element_calc2 = mecalc2
556 @abstractmethod
557 def create_matrix_element_calculators(self):
558 """Create the desired site matrix element calculators."""
560 @abstractmethod
561 def get_spincomponent(self):
562 """Define how to rotate the spins via the spin component (μν)."""
564 @abstractmethod
565 def calculate_eigenvalue_dependent_weights(self, kptpair):
566 """Calculate w_(ε_nks,ε_n'k+qs') for band and spin transitions myt."""
568 def get_band_and_spin_transitions(self):
569 return super().get_band_and_spin_transitions(
570 self.get_spincomponent(),
571 nbands=self.nbands, bandsummation=self.bandsummation)
573 def __call__(self, q_c):
574 """Calculate the site pair function for a given wave vector q_c."""
575 self.context.print(self.get_info_string(q_c))
576 site_pair_function = SitePairFunction(q_c, self.sites)
577 self._integrate(site_pair_function, self.transitions)
578 return site_pair_function
580 def add_integrand(self, kptpair, weight, site_pair_function):
581 r"""Add the site pair function integrand of the outer k-point integral.
583 The integrand is given by (see gpaw.response.pair_integrator)
584 __ __
585 \ \ /
586 (...)_k = V0 / / | σ^μ_ss' σ^ν_s's w_(ε_nks,ε_n'k+qs')
587 ‾‾ ‾‾ \ \
588 n,n' s,s' × f^a_(nks,n'k+qs') g^b_(n'k+qs',nks) |
589 /
591 where V0 is the cell volume.
592 """
593 # Calculate the product of site matrix elements
594 q_c = site_pair_function.q_c
595 matrix_element1 = self.matrix_element_calc1(kptpair, q_c)
596 if self.matrix_element_calc2 is self.matrix_element_calc1:
597 matrix_element2 = matrix_element1
598 else:
599 matrix_element2 = self.matrix_element_calc2(kptpair, q_c)
600 f_mytap = matrix_element1.local_array_view
601 g_mytap = matrix_element2.local_array_view
602 fgcc_mytabp = f_mytap[:, :, np.newaxis] * g_mytap.conj()[:, np.newaxis]
604 # Sum over local transitions, weighted by the spin matrices and
605 # eigenvalue-dependent weights
606 scomps_myt = get_smat_components(
607 self.get_spincomponent(), *kptpair.get_local_spin_indices())
608 weps_myt = self.calculate_eigenvalue_dependent_weights(kptpair)
609 x_myt = scomps_myt * weps_myt # σ^μ_ss' σ^ν_s's w_(ε_nks,ε_n'k+qs')
610 integrand_abp = np.einsum('t, tabp -> abp', x_myt, fgcc_mytabp)
611 # Sum over distributed transitions
612 kptpair.tblocks.blockcomm.sum(integrand_abp)
613 # Add integrand to output array
614 site_pair_function.array[:] += self.gs.volume * weight * integrand_abp
616 def get_info_string(self, q_c):
617 """Get information about the calculation"""
618 info_list = ['',
619 'Calculating site pair function with:'
620 f' q_c: [{q_c[0]}, {q_c[1]}, {q_c[2]}]',
621 self.get_band_and_transitions_info_string(
622 self.nbands, len(self.transitions)),
623 '',
624 self.get_basic_info_string()]
625 return '\n'.join(info_list)
628class TwoParticleSiteSumRuleCalculator(SitePairFunctionCalculator):
629 r"""Calculator for two-particle site sum rules.
631 For any set of site matrix elements f^a and g^b, one may define a two-
632 particle site sum rule,
633 __ __ __
634 1 \ \ \ /
635 ̄x_ab(q) = ‾‾‾ / / / | σ^μ_ss' σ^ν_s's (f_nks - f_n'k+qs')
636 N_k ‾‾ ‾‾ ‾‾ \ \
637 k n,n' s,s' × f^a_(nks,n'k+qs') g^b_(n'k+qs',nks) |
638 /
640 that is, with eigenvalue-dependent weights
642 w_(ε_nks,ε_n'k+qs') = f_nks - f_n'k+qs'
643 """
644 @staticmethod
645 def calculate_eigenvalue_dependent_weights(kptpair):
646 return kptpair.ikpt1.f_myt - kptpair.ikpt2.f_myt # df_myt
649class TwoParticleSiteMagnetizationCalculator(TwoParticleSiteSumRuleCalculator):
650 r"""Calculator for the two-particle site magnetization sum rule.
652 The site magnetization can be calculated from the site pair densities via
653 the following sum rule [publication in preparation]:
654 __ __
655 1 \ \
656 ̄n_ab^z(q) = ‾‾‾ / / (f_nk↑ - f_mk+q↓) n^a_(nk↑,mk+q↓) n^b_(mk+q↓,nk↑)
657 N_k ‾‾ ‾‾
658 k n,m
660 = δ_(a,b) n_a^z
662 This is directly related to the sum rule of the χ^(+-) spin component of
663 the four-component susceptibility tensor.
664 """
665 def create_matrix_element_calculators(self):
666 site_pair_density_calc = SitePairDensityCalculator(
667 self.gs, self.context, self.sites)
668 return site_pair_density_calc, site_pair_density_calc
670 def get_spincomponent(self):
671 return '+-'
674class TwoParticleSiteZeemanEnergyCalculator(
675 TwoParticleSiteMagnetizationCalculator):
676 r"""Calculator for the two-particle site Zeeman energy sum rule.
678 The site Zeeman energy can be calculated from the site pair density and
679 site Zeeman pair energy via the following sum rule [publication in
680 preparation]:
681 __ __
682 ˍ 1 \ \ /
683 E_ab^Z(q) = ‾‾‾ / / | (f_nk↑ - f_mk+q↓)
684 N_k ‾‾ ‾‾ \ \
685 k n,m × E^(Z,a)_(nk↑,mk+q↓) n^b_(mk+q↓,nk↑) |
686 /
687 = δ_(a,b) E_a^Z
688 """
689 def create_matrix_element_calculators(self):
690 site_zeeman_pair_energy_calc = SiteZeemanPairEnergyCalculator(
691 self.gs, self.context, self.sites, rshewmin=1e-8)
692 site_pair_density_calc = SitePairDensityCalculator(
693 self.gs, self.context, self.sites)
694 return site_zeeman_pair_energy_calc, site_pair_density_calc
697class HeisenbergExchangeCalculator(SitePairFunctionCalculator):
698 r"""Calculator for the site-projected Heisenberg exchange.
700 The Heisenberg exchange parameters can be calculated as a function of the
701 wave vector q, by projecting the exchange field J(r,r') onto a series of
702 magnetic sites. The site pair function which follows is given by
703 __ __ /
704 _ 2 \ \ | f_nk↑ - f_mk+q↓
705 J_ab(q) = - ‾‾‾ / / | ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
706 N_k ‾‾ ‾‾ \ ε_nk↑ - ε_mk+q↓ \
707 k n,m |
708 × d^(xc,a)_(nk↑,mk+q↓) d^(xc,b)_(mk+q↓,nk↑) |
709 /
710 where d^(xc,a) is the site spin pair energy, see [publication in
711 preparation] and [J. Phys.: Condens. Matter 35 (2023) 105802].
712 """
713 def __call__(self, q_c):
714 out = super().__call__(q_c)
715 if np.allclose(q_c, 0.):
716 # Symmetrize reciprocity [J^ab(q)]^*=J^ab(-q)
717 J_abp = out.array
718 out.array[:] = (J_abp + J_abp.conj()) / 2.
719 out.array *= Hartree # Ha -> eV
720 return out
722 def create_matrix_element_calculators(self):
723 mcalc = SiteSpinPairEnergyCalculator(
724 self.gs, self.context, self.sites, rshewmin=1e-8)
725 return mcalc, mcalc
727 def get_spincomponent(self):
728 return '+-'
730 @staticmethod
731 def calculate_eigenvalue_dependent_weights(kptpair):
732 """Calculate the eigenvalue-dependent weights.
734 Calculates
736 f_nks - f_mk+qs' f_mk+qs' - f_nks
737 ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ = ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
738 ε_nks - ε_mk+qs' ε_mk+qs' - ε_nks
740 weighted by a prefactor of -2.
741 """
742 nom_myt = kptpair.df_myt # df = (f_n'k's' - f_nks)
743 denom_myt = kptpair.deps_myt # dε = (ε_n'k's' - ε_nks)
744 regularize_intraband_transitions(denom_myt, kptpair)
745 return -2 * nom_myt / denom_myt