Coverage for gpaw/hyperfine.py: 97%
191 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1"""Hyperfine parameters.
3See:
5 First-principles calculations of defects in oxygen-deficient
6 silica exposed to hydrogen
8 Peter E. Blöchl
10 Phys. Rev. B 62, 6158 – Published 1 September 2000
12 https://doi.org/10.1103/PhysRevB.62.6158
14"""
15import argparse
16from math import pi
17from typing import List
19import ase.units as units
20import numpy as np
21from scipy.integrate import simpson
23from gpaw import get_scipy_version
24from gpaw.atom.aeatom import Channel
25from gpaw.atom.configurations import configurations
26from gpaw.atom.radialgd import RadialGridDescriptor
27from gpaw.calculator import GPAW
28from gpaw.gaunt import gaunt
29from gpaw.grid_descriptor import GridDescriptor
30from gpaw.pw.descriptor import PWDescriptor
31from gpaw.setup import Setup
32from gpaw.typing import Array1D, Array2D, Array3D
33from gpaw.utilities import unpack_density
34from gpaw.xc.functional import XCFunctional
36# Fine-structure constant: (~1/137)
37alpha = 0.5 * units._mu0 * units._c * units._e**2 / units._hplanck
39G_FACTOR_E = 2.00231930436256
42def hyperfine_parameters(calc: GPAW,
43 exclude_core=False) -> Array3D:
44 r"""Calculate isotropic and anisotropic hyperfine coupling parameters.
46 One tensor (:math:`A_{ij}`) per atom is returned in eV units.
47 In Hartree atomic units, we have the isotropic part
48 :math:`a = \text{Tr}(\mathbf{A}) / 3`:
50 .. math::
52 a = \frac{2 \alpha^2 g_e m_e}{3 m_p}
53 \int \delta_T(\mathbf{r}) \rho_s(\mathbf{r}) d\mathbf{r},
55 and the anisotropic part :math:`\mathbf{A} - a`:
57 .. math::
59 \frac{\alpha^2 g_e m_e}{4 \pi m_p}
60 \int \frac{3 r_i r_j - \delta_{ij} r^2}{r^5}
61 \rho_s(\mathbf{r}) d\mathbf{r}.
63 Remember to multiply each tensor by the g-factors of the nuclei.
65 Use ``exclude_core=True`` to exclude contribution from "frozen" core.
66 """
67 dens = calc.density
68 nt_sR = dens.nt_sG
69 A_avv = smooth_part(
70 nt_sR[0] - nt_sR[1],
71 dens.gd,
72 calc.atoms.get_scaled_positions())
74 D_asp = calc.density.D_asp
75 for a, D_sp in D_asp.items():
76 density_sii = unpack_density(D_sp)
77 setup = calc.wfs.setups[a]
79 A_vv = paw_correction(density_sii,
80 setup,
81 calc.hamiltonian.xc,
82 exclude_core)
84 A_avv[a] += A_vv
86 A_avv *= pi * alpha**2 * G_FACTOR_E * units._me / units._mp * units.Ha
88 return A_avv
91def smooth_part(spin_density_R: Array3D,
92 gd: GridDescriptor,
93 spos_ac: Array2D,
94 ecut: float = None) -> Array3D:
95 """Contribution from pseudo spin-density."""
96 pd = PWDescriptor(ecut, gd)
97 spin_density_G = pd.fft(spin_density_R)
98 G_Gv = pd.get_reciprocal_vectors(add_q=False)
99 # eiGR_aG = np.exp(-1j * spos_ac.dot(gd.cell_cv).dot(G_Gv.T))
100 eiGR_aG = np.exp(-1j * spos_ac @ gd.cell_cv @ G_Gv.T)
102 # Isotropic term:
103 W1_a = pd.integrate(spin_density_G, eiGR_aG) / gd.dv * (2 / 3)
105 spin_density_G[0] = 0.0
106 G2_G = pd.G2_qG[0].copy()
107 G2_G[0] = 1.0
108 spin_density_G /= G2_G
110 # Anisotropic term:
111 W_vva = np.empty((3, 3, len(spos_ac)))
112 for v1 in range(3):
113 for v2 in range(3):
114 W_a = pd.integrate(G_Gv[:, v1] * G_Gv[:, v2] * spin_density_G,
115 eiGR_aG)
116 W_vva[v1, v2] = -W_a / gd.dv
118 W_a = np.trace(W_vva) / 3
119 for v in range(3):
120 W_vva[v, v] -= W_a
121 W_vva[v, v] += W1_a
123 return W_vva.transpose((2, 0, 1))
126# Normalization constants for xy, yz, 3z^2-r^2, xz, x^2-y^2:
127Y2_m = (np.array([15 / 4,
128 15 / 4,
129 5 / 16,
130 15 / 4,
131 15 / 16]) / pi)**0.5
132# Second derivatives:
133Y2_mvv = np.array([[[0, 1, 0],
134 [1, 0, 0],
135 [0, 0, 0]],
136 [[0, 0, 0],
137 [0, 0, 1],
138 [0, 1, 0]],
139 [[-2, 0, 0],
140 [0, -2, 0],
141 [0, 0, 4]],
142 [[0, 0, 1],
143 [0, 0, 0],
144 [1, 0, 0]],
145 [[2, 0, 0],
146 [0, -2, 0],
147 [0, 0, 0]]])
150def paw_correction(density_sii: Array3D,
151 setup: Setup,
152 xc: XCFunctional = None,
153 exclude_core: bool = False) -> Array2D:
154 """Corrections from 1-center expansions of spin-density."""
155 # Spherical part:
156 spin_density_ii = density_sii[0] - density_sii[1]
157 D0_jj = expand(spin_density_ii, setup.l_j, l=0)[0]
159 phit_jg = np.array(setup.data.phit_jg)
160 phi_jg = np.array(setup.data.phi_jg)
162 rgd = setup.rgd
164 # Spin-density from pseudo density:
165 nt0 = phit_jg[:, 0].dot(D0_jj).dot(phit_jg[:, 0]) / (4 * pi)**0.5
167 # All-electron contribution diverges as r^-beta and must be integrated
168 # over a small region of size rT:
169 n0_g = np.einsum('ab, ag, bg -> g', D0_jj, phi_jg, phi_jg) / (4 * pi)**0.5
170 if not exclude_core and setup.Nc > 0 and xc is not None:
171 n0_g += core_contribution(density_sii, setup, xc)
173 beta = 2 * (1 - (1 - (setup.Z * alpha)**2)**0.5)
174 rT = setup.Z * alpha**2
175 n0 = integrate(n0_g, rgd, rT, beta)
177 W1 = (n0 - nt0) * 2 / 3 # isotropic result
179 # Now the anisotropic part from the l=2 part of the density:
180 D2_mjj = expand(spin_density_ii, setup.l_j, 2)
181 dn2_mg = np.einsum('mab, ag, bg -> mg', D2_mjj, phi_jg, phi_jg)
182 dn2_mg -= np.einsum('mab, ag, bg -> mg', D2_mjj, phit_jg, phit_jg)
183 A_m = dn2_mg[:, 1:].dot(rgd.dr_g[1:] / rgd.r_g[1:]) / 5
184 A_m *= Y2_m
185 # W_vv: Array2D = Y2_mvv.T.dot(A_m) # type: ignore
186 W_vv = Y2_mvv.T @ A_m
187 W = np.trace(W_vv) / 3
188 for v in range(3):
189 W_vv[v, v] -= W
190 W_vv[v, v] += W1
192 return W_vv
195def expand(D_ii: Array2D,
196 l_j: List[int],
197 l: int) -> Array3D:
198 """Get expansion coefficients."""
199 G_LLm = gaunt(lmax=2)[:, :, l**2:(l + 1)**2]
200 D_mjj = np.empty((2 * l + 1, len(l_j), len(l_j)))
201 i1a = 0
202 for j1, l1 in enumerate(l_j):
203 i1b = i1a + 2 * l1 + 1
204 i2a = 0
205 for j2, l2 in enumerate(l_j):
206 i2b = i2a + 2 * l2 + 1
207 D_mjj[:, j1, j2] = np.einsum('ab, abm -> m',
208 D_ii[i1a:i1b, i2a:i2b],
209 G_LLm[l1**2:(l1 + 1)**2,
210 l2**2:(l2 + 1)**2])
211 i2a = i2b
212 i1a = i1b
213 return D_mjj
216def delta(r: Array1D, rT: float) -> Array1D:
217 """Extended delta function."""
218 return 2 / rT / (1 + 2 * r / rT)**2
221def integrate(n0_g: Array1D,
222 rgd: RadialGridDescriptor,
223 rT: float,
224 beta: float) -> float:
225 """Take care of r^-beta divergence."""
226 r_g = rgd.r_g
227 a_i = np.polyfit(r_g[1:5], n0_g[1:5] * r_g[1:5]**beta, 3)
228 r4 = r_g[4]
229 dr = rT / 50
230 n = max(int(r4 / dr / 2) * 2 + 1, 3)
231 r_j = np.linspace(0, r4, n)
233 b_i = np.arange(3, -1, -1) + 1 - beta
234 d0 = 2 / rT # delta(0, rT)
235 d1 = -8 / rT**2
236 n0 = a_i.dot(d0 * r4**b_i / b_i + d1 * r4**(b_i + 1) / (b_i + 1))
238 d_j = delta(r_j, rT) - (d0 + d1 * r_j)
240 head_j = d_j * np.polyval(a_i, r_j)
241 head_j[1:] *= r_j[1:]**-beta
242 n0 += simpson(head_j, x=r_j)
244 tail_g = n0_g[4:] * delta(r_g[4:], rT)
245 if get_scipy_version() >= [1, 11]:
246 n0 += simpson(tail_g, x=r_g[4:])
247 else:
248 n0 += simpson(tail_g, x=r_g[4:], even='first')
250 return n0
253def core_contribution(density_sii: Array3D,
254 setup: Setup,
255 xc: XCFunctional) -> Array1D:
256 """Calculate spin-density from "frozen" core."""
257 # Spherical part:
258 D_sjj = [expand(density_ii, setup.l_j, 0)[0]
259 for density_ii in density_sii]
260 phi_jg = np.array(setup.data.phi_jg)
261 rgd = setup.rgd
263 # Densities with frozen core:
264 n_sg = np.einsum('ag, sab, bg -> sg',
265 phi_jg, D_sjj, phi_jg) / (4 * pi)**0.5
266 n_sg += setup.data.nc_g * (0.5 / (4 * pi)**0.5)
268 # Potential:
269 v_sg = np.zeros_like(n_sg)
270 xc.calculate_spherical(rgd, n_sg, v_sg)
271 vr_sg = v_sg * rgd.r_g
272 vr_sg -= setup.Z
273 vr_sg += rgd.poisson(n_sg.sum(axis=0))
275 # Find first bound s-state included in PAW potential:
276 for n, l in zip(setup.n_j, setup.l_j):
277 if l == 0:
278 assert n > 0
279 n0 = n
280 break
281 else:
282 assert False, (setup.n_j, setup.l_j)
284 # Initial guesses for core s-states:
285 eigs = [e for n, l, f, e in configurations[setup.symbol][1]
286 if n < n0 and l == 0]
288 # Solve spherical scalar-relativistic Schrödinger equation:
289 core_spin_density_g = rgd.zeros()
290 sign = 1.0
291 for vr_g in vr_sg:
292 channel = Channel(l=0, f_n=[1] * len(eigs))
293 channel.e_n = eigs
294 channel.phi_ng = rgd.empty(len(eigs))
295 channel.solve2(vr_g, rgd=rgd, scalar_relativistic=True, Z=setup.Z)
296 assert channel.solve2ok
297 core_spin_density_g += sign * channel.calculate_density()
298 sign = -1.0
300 return core_spin_density_g
303# From https://en.wikipedia.org/wiki/Gyromagnetic_ratio
304# Units: MHz/T
305gyromagnetic_ratios = {'H': (1, 42.577478518),
306 'He': (3, -32.434),
307 'Li': (7, 16.546),
308 'C': (13, 10.7084),
309 'N': (14, 3.077),
310 'O': (17, -5.772),
311 'F': (19, 40.052),
312 'Na': (23, 11.262),
313 'Al': (27, 11.103),
314 'Si': (29, -8.465),
315 'P': (31, 17.235),
316 'Fe': (57, 1.382),
317 'Cu': (63, 11.319),
318 'Zn': (67, 2.669),
319 'Xe': (129, -11.777)}
322def main(argv: List[str] = None) -> None:
323 """Command-line interface."""
324 parser = argparse.ArgumentParser(
325 prog='python3 -m gpaw.hyperfine',
326 description='Calculate hyperfine parameters.')
327 add = parser.add_argument
328 add('file', metavar='input-file',
329 help='GPW-file (with or without wave functions).')
330 add('-g', '--g-factors',
331 help='G-factors. Example: "-g H:5.6,O:-0.76".')
332 add('-u', '--units', default='ueV', choices=['ueV', 'MHz'],
333 help='Units. Must be "ueV" (micro-eV, default) or "MHz".')
334 add('-x', '--exclude-core', action='store_true')
335 add('-d', '--diagonalize', action='store_true',
336 help='Show eigenvalues of tensor.')
338 args = parser.parse_intermixed_args(argv)
340 calc = GPAW(args.file)
341 atoms = calc.get_atoms()
343 symbols = atoms.symbols
344 magmoms = atoms.get_magnetic_moments()
345 total_magmom = atoms.get_magnetic_moment()
347 g_factors = {symbol: ratio * 1e6 * 4 * pi * units._mp / units._e
348 for symbol, (n, ratio) in gyromagnetic_ratios.items()}
350 if args.g_factors:
351 for symbol, value in (part.split(':')
352 for part in args.g_factors.split(',')):
353 g_factors[symbol] = float(value)
355 if args.units == 'ueV':
356 scale = 1e6
357 unit = 'μeV'
358 else:
359 scale = units._e / units._hplanck * 1e-6
360 unit = 'MHz'
362 A_avv = hyperfine_parameters(calc, args.exclude_core)
364 print('Hyperfine coupling paramters '
365 f'in {unit}:\n')
367 if args.diagonalize:
368 columns = ['1.', '2.', '3.']
369 else:
370 columns = ['xx', 'yy', 'zz', 'yz', 'xz', 'xy']
371 print(' atom magmom ', ' '.join(columns))
373 used = {}
374 for a, A_vv in enumerate(A_avv):
375 symbol = symbols[a]
376 magmom = magmoms[a]
377 g_factor = g_factors.get(symbol, 1.0)
378 used[symbol] = g_factor
379 A_vv *= g_factor * scale
380 if args.diagonalize:
381 numbers = np.linalg.eigvalsh(A_vv).tolist()
382 else:
383 numbers = [A_vv[0, 0], A_vv[1, 1], A_vv[2, 2],
384 A_vv[1, 2], A_vv[0, 2], A_vv[0, 1]]
385 print(f'{a:3} {symbol:>2} {magmom:6.3f}',
386 ''.join(f'{x:9.2f}' for x in numbers))
388 print('\nCore correction',
389 'NOT included!' if args.exclude_core else 'included')
390 print(f'Total magnetic moment: {total_magmom:.3f}')
392 print('\nG-factors used:')
393 for symbol, g in used.items():
394 print(f'{symbol:2} {g:10.3f}')
397if __name__ == '__main__':
398 main()