Coverage for gpaw/test/sphere/test_integrate.py: 100%
108 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1import pytest
3import numpy as np
5from ase.units import Bohr
7from gpaw import GPAW
9from gpaw.sphere.integrate import (integrate_lebedev, radial_trapz,
10 radial_truncation_function,
11 periodic_truncation_function,
12 spherical_truncation_function_collection,
13 default_spherical_drcut,
14 find_volume_conserving_lambd)
17def generate_analytical_integrals():
18 """Return pair of functions that (1) evaluates a given f(r) (2) integrates
20 r_cut
21 /
22 | r^2 dr f(r)
23 /
24 0
26 analytically together with a relative tolerance.
27 """
29 def linear(r):
30 f = 2. - r
31 f[r > 2.] = 0.
32 return f
34 def integrate_linear(rcut):
35 if rcut > 2.:
36 return 4. / 3.
37 else:
38 return (2 / 3. - rcut / 4.) * rcut**3.
40 def gaussian(r):
41 return np.exp(-r**2 / 4)
43 def integrate_gaussian(rcut):
44 from scipy.special import erf
45 return 2 * np.sqrt(np.pi) * erf(rcut / 2) \
46 - 2 * rcut * np.exp(-rcut**2 / 4)
48 def lorentzian(r):
49 return 1 / (r**2 + 4)
51 def integrate_lorentzian(rcut):
52 return rcut - 2. * np.arctan(rcut / 2)
54 return [
55 # Exact numerical integration with linear interpolation
56 (linear, integrate_linear, 1e-8),
57 # Approximate numerical integration
58 (gaussian, integrate_gaussian, 1e-5),
59 (lorentzian, integrate_lorentzian, 1e-5)
60 ]
63@pytest.mark.parametrize('analytical_integral',
64 generate_analytical_integrals())
65def test_radial_trapz(analytical_integral):
66 func, integrate, rtol = analytical_integral
67 # Create radial grid and evaluate the function
68 r_g = np.linspace(0., 3.0, 250)
69 f_g = func(r_g)
71 # Compute the radial integral for a number of different cutoffs:
72 for ncut in np.arange(100, 260, 10):
73 rcut = r_g[ncut - 1]
74 ref = integrate(rcut)
75 result = radial_trapz(f_g[:ncut], r_g[:ncut])
76 assert abs(result - ref) <= 1e-8 + rtol * ref
79@pytest.mark.parametrize('rc', np.linspace(0.5, 4.5, 9))
80def test_smooth_truncation_function(rc):
81 # Define radial grid
82 r_g = np.linspace(0., 5.0, 501)
83 # Calculate spherical volume corresponding to the cutoff
84 ref = 4 * np.pi * rc**3. / 3.
85 # Test drc from grid spacing x2 to grid spacing x50
86 for drc in np.linspace(0.02, 0.5, 10):
87 theta_g = radial_truncation_function(r_g, rc, drc)
88 # Calculate spherical volume with truncation function
89 vol = 4 * np.pi * radial_trapz(theta_g, r_g)
90 assert abs(vol - ref) <= 1e-8 + 1e-4 * ref
93def test_fe_augmentation_sphere(gpw_files):
94 # Extract the spherical grid information from the iron fixture
95 calc = GPAW(gpw_files['fe_pw'], txt=None)
96 setup = calc.setups[0]
97 rgd = setup.rgd
98 Y_nL = setup.xc_correction.Y_nL
100 # Create a function which is unity over the entire grid
101 f_ng = np.ones((Y_nL.shape[0], rgd.N))
103 # Integrate f(r) with different cutoffs, to check that the volume is
104 # correctly recovered
105 r_g = rgd.r_g
106 for rcut in np.linspace(0.1 / Bohr, np.max(r_g), 100):
107 ref = 4 * np.pi * rcut**3. / 3.
109 # Do standard numerical truncation
110 # Integrate angular components, then radial
111 f_g = integrate_lebedev(f_ng)
112 vol = rgd.integrate_trapz(f_g, rcut=rcut)
113 assert abs(vol - ref) <= 1e-8 + 1e-6 * ref
114 # Integrate radial components, then angular
115 f_n = rgd.integrate_trapz(f_ng, rcut=rcut)
116 vol = integrate_lebedev(f_n)
117 assert abs(vol - ref) <= 1e-8 + 1e-6 * ref
119 # Integrate f(r) θ(r<rc) using a smooth truncation function
120 if rcut > np.max(setup.rcut_j):
121 # This method relies on a sufficiently dense grid sampling to be
122 # accurate, so we only test values inside the augmentation sphere
123 continue
124 theta_g = radial_truncation_function(r_g, rcut)
125 ft_ng = f_ng * theta_g[np.newaxis]
126 # Integrate angular components, then radial
127 ft_g = integrate_lebedev(ft_ng)
128 vol = rgd.integrate_trapz(ft_g)
129 assert abs(vol - ref) <= 1e-8 + 1e-4 * ref
130 # Integrate radial components, then angular
131 ft_n = rgd.integrate_trapz(ft_ng)
132 vol = integrate_lebedev(ft_n)
133 assert abs(vol - ref) <= 1e-8 + 1e-4 * ref
136def test_fe_periodic_truncation_function(gpw_files):
137 # Extract the grid information from the iron fixture
138 calc = GPAW(gpw_files['fe_pw'], txt=None)
139 finegd = calc.density.finegd
140 spos_ac = calc.spos_ac
142 # Integrate θ(r) with different cutoffs, to check that the sphere volume
143 # is correctly recovered
144 a = 2.867 # lattice constant in Å
145 rcut_max = 2 * a / (3 * Bohr) # 2a / 3 in Bohr
146 # Get default drcut corresponding to the coarse real-space grid
147 drcut = default_spherical_drcut(calc.density.gd)
148 for rcut in np.linspace(rcut_max / 6, rcut_max, 13):
149 ref = 4 * np.pi * rcut**3. / 3.
151 # Optimize λ-parameter, generate θ(r) and integrate
152 lambd = find_volume_conserving_lambd(rcut, drcut)
153 theta_r = periodic_truncation_function(finegd, spos_ac[0],
154 rcut, drcut, lambd)
155 vol = finegd.integrate(theta_r)
156 assert abs(vol - ref) <= 1e-8 + 1e-2 * ref
158 # Make sure that the difference between coarse and fine grid drcut is not
159 # too large
160 finedrcut_vol = finegd.integrate(
161 periodic_truncation_function(finegd, spos_ac[0], rcut))
162 assert abs(vol - finedrcut_vol) / finedrcut_vol < 1e-3
163 # Make sure that we get a different value numerically, if we change drcut
164 diff_vol = finegd.integrate(
165 periodic_truncation_function(finegd, spos_ac[0], rcut,
166 drcut=drcut * 1.5))
167 assert abs(vol - diff_vol) > 1e-8
168 # Test that the actual value of the integral changes, if we use a different
169 # λ-parameter
170 diff_vol = finegd.integrate(
171 periodic_truncation_function(finegd, spos_ac[0], rcut,
172 lambd=0.75))
173 assert abs(vol - diff_vol) > 1e-2 * ref
176def test_co_spherical_truncation_function_collection(gpw_files):
177 # Extract grid information from the cobalt fixture
178 calc = GPAW(gpw_files['co_pw'], txt=None)
179 finegd = calc.density.finegd
180 spos_ac = calc.spos_ac
181 drcut = default_spherical_drcut(calc.density.gd)
183 # Generate collection of spherical truncation functions with varrying rcut
184 a = 2.5071
185 c = 4.0695
186 nn_dist = min(a, np.sqrt(a**2 / 3 + c**2 / 4))
187 rcut_j = np.linspace(drcut, 2 * nn_dist / 3, 13)
188 rcut_aj = [rcut_j, rcut_j]
189 stfc = spherical_truncation_function_collection(finegd, spos_ac, rcut_aj,
190 drcut=drcut)
192 # Integrate collection of spherical truncation functions
193 ones_r = finegd.empty()
194 ones_r[:] = 1.
195 vol_aj = {0: np.empty(len(rcut_j)), 1: np.empty(len(rcut_j))}
196 stfc.integrate(ones_r, vol_aj)
198 # Check that the integrated volume matches the spherical volume and an
199 # analogous manual integration
200 for a, spos_c in enumerate(spos_ac):
201 for rcut, vol in zip(rcut_j, vol_aj[a]):
202 ref = 4 * np.pi * rcut**3. / 3.
203 assert abs(vol - ref) <= 1e-8 + 1e-2 * ref
205 # "Manual" integration
206 theta_r = periodic_truncation_function(finegd, spos_c, rcut, drcut)
207 manual_vol = finegd.integrate(theta_r)
208 assert abs(vol - manual_vol) <= 1e-8 + 1e-6 * ref