Coverage for gpaw/response/goldstone.py: 98%
138 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1from __future__ import annotations
3from abc import abstractmethod
4from functools import partial
5from typing import TYPE_CHECKING
7import numpy as np
8from scipy.optimize import minimize
10from gpaw.response.dyson import HXCScaling, DysonEquation, DysonEquations
12if TYPE_CHECKING:
13 from gpaw.response.chiks import SelfEnhancementCalculator
16class GoldstoneScaling(HXCScaling):
17 """Scale the Dyson equation to fulfill a Goldstone condition."""
19 def _calculate_scaling(self, dyson_equations):
20 """Calculate scaling coefficient λ."""
21 self.check_descriptors(dyson_equations)
23 # Find the frequency to determine the scaling from and identify where
24 # the Dyson equation in question is distributed
25 wgs = self.find_goldstone_frequency(
26 dyson_equations.zd.omega_w)
27 wblocks = dyson_equations.zblocks
28 rgs, mywgs = wblocks.find_global_index(wgs)
30 # Let the rank which holds the Goldstone frequency find and broadcast λ
31 lambdbuf = np.empty(1, dtype=float)
32 if wblocks.blockcomm.rank == rgs:
33 lambdbuf[:] = self.find_goldstone_scaling(dyson_equations[mywgs])
34 wblocks.blockcomm.broadcast(lambdbuf, rgs)
35 lambd = lambdbuf[0]
37 return lambd
39 @staticmethod
40 def check_descriptors(dyson_equations):
41 if not (dyson_equations.qpd.optical_limit and
42 dyson_equations.spincomponent in ['+-', '-+']):
43 raise ValueError(
44 'The Goldstone condition only applies to χ^(+-)(q=0).')
46 @abstractmethod
47 def find_goldstone_frequency(self, omega_w):
48 """Determine frequency index for the Goldstone condition."""
50 @abstractmethod
51 def find_goldstone_scaling(self, dyson_equation: DysonEquation) -> float:
52 """Calculate the Goldstone scaling parameter λ."""
55class FMGoldstoneScaling(GoldstoneScaling):
56 """Fulfil ferromagnetic Goldstone condition."""
58 @staticmethod
59 def find_goldstone_frequency(omega_w):
60 """Ferromagnetic Goldstone condition is based on χ^(+-)(ω=0)."""
61 wgs = np.abs(omega_w).argmin()
62 assert abs(omega_w[wgs]) < 1.e-8, \
63 "Frequency grid needs to include ω=0."
64 return wgs
66 @staticmethod
67 def find_goldstone_scaling(dyson_equation):
68 return find_fm_goldstone_scaling(dyson_equation)
71class NewFMGoldstoneScaling(FMGoldstoneScaling):
72 """Fulfil Goldstone condition by maximizing a^(+-)(ω=0)."""
74 def __init__(self,
75 lambd: float | None = None,
76 m_G: np.ndarray | None = None):
77 """Construct the scaling object.
79 If the λ-parameter hasn't yet been calculated, the (normalized)
80 spin-polarization |m> is needed in order to extract the acoustic
81 magnon mode lineshape, a^(+-)(ω=0) = -Im[<m|χ^(+-)(ω=0)|m>]/π.
82 """
83 super().__init__(lambd=lambd)
84 self.m_G = m_G
86 @classmethod
87 def from_xi_calculator(cls, xi_calc: SelfEnhancementCalculator):
88 """Construct scaling object with |m> consistent with a xi_calc."""
89 return cls(m_G=cls.calculate_m(xi_calc))
91 @staticmethod
92 def calculate_m(xi_calc: SelfEnhancementCalculator):
93 """Calculate the normalized spin-polarization |m>."""
94 from gpaw.response.localft import (LocalFTCalculator,
95 add_spin_polarization)
96 localft_calc = LocalFTCalculator.from_rshe_parameters(
97 xi_calc.gs, xi_calc.context,
98 rshelmax=xi_calc.rshelmax, rshewmin=xi_calc.rshewmin)
99 qpd = xi_calc.get_pw_descriptor(q_c=[0., 0., 0.])
100 nz_G = localft_calc(qpd, add_spin_polarization)
101 return nz_G / np.linalg.norm(nz_G)
103 def find_goldstone_scaling(self, dyson_equation):
104 assert self.m_G is not None, \
105 'Please supply spin-polarization to calculate λ'
107 def acoustic_antispectrum(lambd):
108 """Calculate -a^(+-)(ω=0)."""
109 return - calculate_acoustic_spectrum(
110 lambd, dyson_equation, self.m_G)
112 # Maximize a^(+-)(ω=0)
113 res = minimize(acoustic_antispectrum, x0=[1.], bounds=[(0.1, 10.)])
114 assert res.success
115 return res.x[0]
118class RefinedFMGoldstoneScaling(HXCScaling):
119 """Ensures that a^(+-)(ω) has a maximum in ω=0."""
121 def __init__(self,
122 lambd: float | None = None,
123 base_scaling: NewFMGoldstoneScaling | None = None):
124 """Construct the scaling object.
126 If the λ-parameter hasn't yet been calculated, we use a base scaling
127 class, which calculates λ approximately (and gives us access to |m>),
128 thus providing a starting point for the refinement.
129 """
130 super().__init__(lambd=lambd)
131 self._base_scaling = base_scaling
133 @classmethod
134 def from_xi_calculator(cls, xi_calc: SelfEnhancementCalculator):
135 return cls(
136 base_scaling=NewFMGoldstoneScaling.from_xi_calculator(xi_calc))
138 @property
139 def m_G(self):
140 assert self._base_scaling is not None
141 assert self._base_scaling.m_G is not None
142 return self._base_scaling.m_G
144 def _calculate_scaling(self, dyson_equations: DysonEquations) -> float:
145 """Calculate the scaling coefficient λ."""
146 # First we calculate the base scaling based on a^(+-)(ω=0)
147 assert self._base_scaling is not None
148 self._base_scaling.calculate_scaling(dyson_equations)
149 base_lambd = self._base_scaling.lambd
151 # Secondly, we extract the spectral peak position by performing a
152 # parabolic fit to the five points with lowest |ω|.
153 omega_W = dyson_equations.zd.omega_w
154 wblocks = dyson_equations.zblocks
155 fiveW_w = np.argpartition(np.abs(omega_W), 5)[:5]
156 omega_w = omega_W[fiveW_w]
158 def near_acoustic_spectrum(lambd):
159 a_w = np.empty(5, dtype=float)
160 for w, W in enumerate(fiveW_w):
161 wrank, myw = wblocks.find_global_index(W)
162 if wblocks.blockcomm.rank == wrank:
163 a_w[w] = calculate_acoustic_spectrum(
164 lambd, dyson_equations[myw], self.m_G)
165 wblocks.blockcomm.broadcast(a_w[w:w + 1], wrank)
166 return a_w
168 def acoustic_magnon_frequency(lambd):
169 a_w = near_acoustic_spectrum(lambd)
170 a, b, c = np.polyfit(omega_w, a_w, 2)
171 return -b / (2 * a)
173 # Lastly, we minimize the (absolute) peak frequency |ω_0| to obtain the
174 # refined λ. To do so efficiently, we define a (hyperbolic) cost
175 # function which is linear in |ω_0| in the meV range, but parabolic in
176 # the μeV range, such that derivatives remain smooth at the minimum.
178 def cost_function(lambd):
179 return np.sqrt(5e-5 + acoustic_magnon_frequency(lambd)**2)
181 res = minimize(cost_function, x0=[base_lambd],
182 bounds=[(base_lambd * 0.975, base_lambd * 1.025)])
183 assert res.success, res
184 return res.x[0]
187class AFMGoldstoneScaling(GoldstoneScaling):
188 """Fulfil antiferromagnetic Goldstone condition."""
190 @staticmethod
191 def find_goldstone_frequency(omega_w):
192 """Antiferromagnetic Goldstone condition is based on ω->0^+."""
193 # Set ω<=0. to np.inf
194 omega1_w = np.where(omega_w < 1.e-8, np.inf, omega_w)
195 # Sort for the two smallest positive frequencies
196 omega2_w = np.partition(omega1_w, 1)
197 # Find original index of second smallest positive frequency
198 wgs = np.abs(omega_w - omega2_w[1]).argmin()
199 return wgs
201 @staticmethod
202 def find_goldstone_scaling(dyson_equation):
203 return find_afm_goldstone_scaling(dyson_equation)
206def find_fm_goldstone_scaling(dyson_equation):
207 """Find goldstone scaling of the kernel by ensuring that the
208 macroscopic inverse enhancement function has a root in (q=0, omega=0)."""
209 fnct = partial(calculate_macroscopic_kappa,
210 dyson_equation=dyson_equation)
212 def is_converged(kappaM):
213 return abs(kappaM) < 1e-7
215 return find_root(fnct, is_converged)
218def find_afm_goldstone_scaling(dyson_equation):
219 """Find goldstone scaling of the kernel by ensuring that the
220 macroscopic magnon spectrum vanishes at q=0. for finite frequencies."""
221 fnct = partial(calculate_macroscopic_spectrum,
222 dyson_equation=dyson_equation)
224 def is_converged(SM):
225 # We want the macroscopic spectrum to be strictly positive
226 return 0. < SM < 1.e-7
228 return find_root(fnct, is_converged)
231def find_root(fnct, is_converged):
232 """Find the root f(λ)=0, where the scaling parameter λ~1.
234 |f(λ)| is minimized iteratively, assuming that f(λ) is continuous and
235 monotonically decreasing with λ for λ∊]0.1, 10[.
236 """
237 lambd = 1. # initial guess for the scaling parameter
238 value = fnct(lambd)
239 lambd_incr = 0.1 * np.sign(value) # Increase λ to decrease f(λ)
240 while not is_converged(value) or abs(lambd_incr) > 1.e-7:
241 # Update λ
242 lambd += lambd_incr
243 if lambd <= 0.1 or lambd >= 10.:
244 raise Exception(f'Found an invalid λ-value of {lambd:.4f}')
245 # Update value and refine increment, if we have passed f(λ)=0
246 value = fnct(lambd)
247 if np.sign(value) != np.sign(lambd_incr):
248 lambd_incr *= -0.2
249 return lambd
252def calculate_macroscopic_kappa(lambd, dyson_equation):
253 """Invert dyson equation and calculate the inverse enhancement function."""
254 chi_GG = dyson_equation.invert(lambd=lambd)
255 return (dyson_equation.chiks_GG[0, 0] / chi_GG[0, 0]).real
258def calculate_macroscopic_spectrum(lambd, dyson_equation):
259 """Invert dyson equation and extract the macroscopic spectrum."""
260 chi_GG = dyson_equation.invert(lambd=lambd)
261 return - chi_GG[0, 0].imag / np.pi
264def calculate_acoustic_spectrum(lambd, dyson_equation, m_G):
265 """Invert the dyson equation and extract the acoustic spectrum."""
266 chi_GG = dyson_equation.invert(lambd=lambd)
267 chi_projection = np.conj(m_G) @ chi_GG @ m_G
268 return - chi_projection.imag / np.pi