Coverage for gpaw/test/hyperfine/test_hyperfine.py: 97%
126 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1from math import pi
3import numpy as np
4import pytest
5from scipy.special import expn
6import ase.units as units
8from gpaw.grid_descriptor import GridDescriptor
9from gpaw.hyperfine import (hyperfine_parameters, paw_correction, smooth_part,
10 integrate, alpha, G_FACTOR_E, core_contribution)
11from gpaw import GPAW
12from gpaw.atom.aeatom import AllElectronAtom
13from gpaw.atom.radialgd import RadialGridDescriptor
14from gpaw.lfc import LFC
15from gpaw.setup import create_setup
16from gpaw.xc import XC
19@pytest.mark.serial
20def test_thomson_integral():
21 dr = 0.0001
22 rgd = RadialGridDescriptor.new('equidistant', dr, 16000)
23 beta = 0.0
24 rT = 0.001
25 n0_g = np.exp(-rgd.r_g)
27 n0 = integrate(n0_g, rgd, rT, beta)
28 ref = expn(2, rT / 2) * np.exp(rT / 2)
29 assert n0 == pytest.approx(ref, 1e-5)
32@pytest.mark.serial
33def test_thomson_integral2():
34 dr = 0.0001
35 rgd = RadialGridDescriptor.new('equidistant', dr, 66000)
36 rT = 0.001
37 beta = 0.001
38 n0_g = rgd.zeros()
39 n0_g[1:] = rgd.r_g[1:]**-beta
40 n0_g[0] = np.nan
42 n0 = integrate(n0_g, rgd, rT, beta)
43 ref = (2 / rT)**beta * np.pi * rT / np.sin(np.pi * rT)
44 assert n0 == pytest.approx(ref, 1e-4)
47class Setup:
48 def __init__(self):
49 self.rgd = RadialGridDescriptor.new('equidistant', 0.01, 400)
50 self.data = Data(self.rgd)
51 self.l_j = [0, 1]
52 self.Z = 1
53 self.Nc = 0
56class Data:
57 def __init__(self, rgd: RadialGridDescriptor):
58 self.phit_jg = [np.exp(-(rgd.r_g / 0.5)**2),
59 np.exp(-(rgd.r_g / 0.5)**2) * rgd.r_g]
60 self.phi_jg = [rgd.zeros(), rgd.zeros()]
63@pytest.fixture
64def things():
65 setup = Setup()
67 n = 40
68 L = n * 0.1
69 gd = GridDescriptor((n, n, n), (L, L, L), (1, 1, 1))
71 splines = []
72 for j, phit_g in enumerate(setup.data.phit_jg):
73 splines.append(setup.rgd.spline(phit_g, 1.5, l=j, points=101))
75 lfc = LFC(gd, [splines])
77 spos_ac = np.zeros((1, 3)) + 0.5
78 lfc.set_positions(spos_ac)
80 return gd, lfc, setup, spos_ac
83@pytest.mark.serial
84def test_gaussian(things):
85 gd, lfc, setup, spos_ac = things
86 spin_density_R = gd.zeros()
87 lfc.add(spin_density_R, {0: np.array([1.0, 0, 0, 0])})
88 spin_density_R *= spin_density_R
90 W_avv = smooth_part(spin_density_R, gd, spos_ac)
91 print(W_avv)
92 assert abs(W_avv[0] - np.eye(3) * W_avv[0, 0, 0]).max() < 1e-7
94 density_sii = np.zeros((2, 4, 4))
95 density_sii[0, 0, 0] = 1.0
96 W1_vv = paw_correction(density_sii, setup)
97 print(W1_vv)
98 assert abs(W_avv[0] + W1_vv).max() < 1e-7
101@pytest.mark.serial
102def test_gaussian2(things):
103 gd, lfc, setup, spos_ac = things
104 spin_density_R = gd.zeros()
105 lfc.add(spin_density_R, {0: np.array([0, 0, 1.0, 0])})
106 spin_density2_R = gd.zeros()
107 lfc.add(spin_density2_R, {0: np.array([0, 1.0, 0, 0])})
108 spin_density_R = spin_density_R * spin_density2_R
110 W_avv = smooth_part(spin_density_R, gd, spos_ac)
111 print(W_avv)
112 assert abs(W_avv[0] - np.array([[0, 0, 0],
113 [0, 0, 1],
114 [0, 1, 0]]) * W_avv[0, 1, 2]).max() < 1e-7
116 density_sii = np.zeros((2, 4, 4))
117 density_sii[0, 2, 1] = 1.0
118 W1_vv = paw_correction(density_sii, setup)
119 print(W1_vv)
120 assert abs(W_avv[0] + W1_vv).max() < 1e-6
123G_FACTOR_PROTON = 5.586
126@pytest.mark.serial
127def test_h(gpw_files):
128 calc = GPAW(gpw_files['h_pw'])
129 A_vv = hyperfine_parameters(calc)[0] * G_FACTOR_PROTON
130 print(A_vv)
132 energy = (2 / 3 * alpha**2 * units.Ha * units._me / units._mp *
133 G_FACTOR_E * G_FACTOR_PROTON) # in eV
134 frequency = energy * units._e / units._hplanck # Hz
135 wavelength = units._c / frequency # meters
136 assert wavelength == pytest.approx(0.211, abs=0.001)
138 energy *= 0.94 # density at nucleus is slightly below 1/pi for PBE
139 print(energy)
140 assert abs(A_vv - np.eye(3) * energy).max() < 1e-7
143def thomson():
144 """Analytic integrals for testing."""
145 from sympy import var, integrate, oo, E, expint
146 x, a, b = var('x, a, b')
147 print(integrate(E**(-b * x) / (1 + x)**2, (x, 0, oo)))
148 print(expint(2, 1.0))
151def test_hyper_core():
152 """Test polarization of "frozen" nitrogen-1s core."""
153 # Calculate 1s spin-density from 2s and 2p polarizarion:
154 setup = create_setup('N')
155 xc = XC('LDA')
156 D_sii = np.zeros((2, 13, 13))
157 D_sii[:, 0, 0] = 1.0
158 D_sii[0, 1:4, 1:4] = np.eye(3)
159 spin_density1_g = core_contribution(D_sii, setup, xc)
161 assert setup.rgd.integrate(spin_density1_g) == pytest.approx(0, abs=1e-11)
163 # Do all-electron calculation:
164 aea = AllElectronAtom('N', spinpol=True)
165 aea.initialize()
166 aea.run()
167 aea.scalar_relativistic = True
168 aea.refine()
170 # 1s:
171 spin_density2_g = (aea.channels[0].phi_ng[0]**2 -
172 aea.channels[1].phi_ng[0]**2) / (4 * pi)
173 # 2s:
174 spin_density3_g = (aea.channels[0].phi_ng[1]**2 -
175 aea.channels[1].phi_ng[1]**2) / (4 * pi)
177 assert aea.rgd.integrate(spin_density2_g) == pytest.approx(0, abs=1e-11)
178 assert aea.rgd.integrate(spin_density3_g) == pytest.approx(0, abs=1e-11)
180 if 0:
181 import matplotlib.pyplot as plt
182 plt.plot(setup.rgd.r_g, spin_density1_g)
183 plt.plot(aea.rgd.r_g, spin_density2_g)
184 plt.show()
186 assert spin_density1_g[0] == pytest.approx(spin_density2_g[0], abs=0.04)