Coverage for gpaw/test/response/test_localft.py: 100%
128 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1"""Test functionality to compute Fourier Transforms with PAW corrections"""
3# General modules
4import numpy as np
5import pytest
7# Script modules
8from ase.units import Ha
10from gpaw import GPAW
11import gpaw.mpi as mpi
12from gpaw.pw.descriptor import PWDescriptor
13from gpaw.kpt_descriptor import KPointDescriptor
14from gpaw.grid_descriptor import GridDescriptor
15from gpaw.lfc import LFC
16from gpaw.atom.radialgd import AERadialGridDescriptor
18from gpaw.response import ResponseGroundStateAdapter, ResponseContext
19from gpaw.response.localft import (LocalFTCalculator, MicroSetup,
20 add_total_density, add_LSDA_Wxc)
21from gpaw.response.pair_functions import get_pw_coordinates
22from gpaw.test.response.test_site_kernels import get_pw_descriptor
25# ---------- Test parametrization ---------- #
27# 1s orbital radii
28testa_a = np.linspace(0.5, 1.5, 10) # a.u.
31def ae_1s_density(r_g, a=1.0):
32 """Construct the radial dependence of the density from a 1s orbital on the
33 radial grid r_g."""
34 assert np.all(r_g >= 0)
35 prefactor = 1 / (np.pi * a**3.)
36 n_g = prefactor * np.exp(-2. * r_g / a)
38 return n_g
41def ae_1s_density_plane_waves(pd, R_v, a=1.0):
42 """Calculate the plane-wave components of the density from a 1s
43 orbital centered at a given position analytically."""
44 # List of all plane waves
45 G_Gv = np.array([pd.G_Qv[Q] for Q in pd.Q_qG[0]])
46 Gnorm_G = np.linalg.norm(G_Gv, axis=1)
48 position_prefactor_G = np.exp(-1.j * np.dot(G_Gv, R_v))
49 atomcentered_n_G = 1 / (1 + (Gnorm_G * a / 2.)**2.)**2.
51 n_G = position_prefactor_G * atomcentered_n_G
53 return n_G
56# ---------- Actual tests ---------- #
58@pytest.mark.response
59@pytest.mark.parametrize("a", testa_a)
60def test_localft_grid_calculator(a):
61 """Test that the LocalGridFTCalculator is able to correctly Fourier
62 transform the all-electron density of an 1s orbital."""
63 # ---------- Inputs ---------- #
65 # Real-space grid
66 lc_to_a_ratio = 10 # lattice constant to orbital radii
67 N_grid_points = 50**3
69 # Plane-wave cutoff
70 relative_ecut = 20 # eV relative to a=1
72 # Test tolerance
73 rtol = 1e-3
75 # ---------- Script ---------- #
77 # Set up atomic position at the center of the unit cell
78 lattice_constant = lc_to_a_ratio * a # a.u.
79 R_v = np.array([lattice_constant, lattice_constant,
80 lattice_constant]) / 2. # Place atom at the center
82 # Set up grid descriptor
83 cell_cv = np.array([[lattice_constant, 0., 0.],
84 [0., lattice_constant, 0.],
85 [0., 0., lattice_constant]])
86 N_c = np.array([int(N_grid_points**(1 / 3.))] * 3)
87 gd = GridDescriptor(N_c, cell_cv=cell_cv, comm=mpi.serial_comm)
89 # Set up plane-wave descriptor
90 qd = KPointDescriptor(np.array([[0., 0., 0.]]))
91 pd = PWDescriptor(relative_ecut / Ha / a**2., gd, complex, qd)
93 # Calculate the plane-wave components analytically
94 ntest_G = ae_1s_density_plane_waves(pd, R_v, a=a)
96 # Calculate the atomic radius at all grid points
97 r_vR = gd.get_grid_point_coordinates()
98 r_R = np.linalg.norm(r_vR - R_v[:, np.newaxis, np.newaxis, np.newaxis],
99 axis=0)
101 # Calculate the all-electron density on the real-space grid
102 n_sR = np.array([ae_1s_density(r_R, a=a)])
104 # Initialize the LocalGridFTCalculator with an empty ground state adapter
105 gs = EmptyGSAdapter() # hack to pass isinstance in constructor
106 context = ResponseContext()
107 localft_calc = LocalFTCalculator.from_rshe_parameters(gs, context,
108 rshelmax=None)
110 # Compute the plane-wave components numerically
111 n_G = localft_calc._calculate(pd, n_sR, add_total_density)
113 # Check validity of results
114 assert np.allclose(n_G, ntest_G, rtol=rtol)
117@pytest.mark.response
118@pytest.mark.parametrize("a", testa_a)
119def test_localft_paw_engine(a):
120 """Test that the LocalPAWFTEngine is able to correctly Fourier
121 transform the all-electron density of an 1s orbital."""
122 # ---------- Inputs ---------- #
124 # Real-space grid
125 lc_to_a_ratio = 10 # lattice constant to orbital radii
126 N_grid_points = 50**3
128 # Radial grid (using standard parameters from Li)
129 rgd_a = 0.0023570226039551583
130 rgd_b = 0.0004528985507246377
131 rcut = 2.0 # a.u.
133 # Plane-wave cutoff
134 relative_ecut = 20 # eV relative to a=1
136 # Settings for the expansion in real spherical harmonics
137 rshe_params_p = [{},
138 {'rshelmax': 0}, # test that only l=0 contributes
139 {'rshewmin': 1e-8}] # test coefficient filter
141 # Test tolerance
142 rtol = 5e-4
144 # ---------- Script ---------- #
146 # Set up atomic position at the center of the unit cell
147 lattice_constant = lc_to_a_ratio * a # a.u.
148 R_v = np.array([lattice_constant, lattice_constant,
149 lattice_constant]) / 2.
150 pos_ac = np.array([[0.5, 0.5, 0.5]]) # Relative atomic positions
152 # Set up grid descriptor
153 cell_cv = np.array([[lattice_constant, 0., 0.],
154 [0., lattice_constant, 0.],
155 [0., 0., lattice_constant]])
156 N_c = np.array([int(N_grid_points**(1 / 3.))] * 3)
157 gd = GridDescriptor(N_c, cell_cv=cell_cv, comm=mpi.serial_comm)
159 # Set up radial grid descriptor extending all the way to the edge of the
160 # unit cell
161 redge = np.sqrt(3) * lattice_constant / 2. # center-corner distance
162 Ng = int(np.floor(redge / (rgd_a + rgd_b * redge)) + 1)
163 rgd = AERadialGridDescriptor(rgd_a, rgd_b, N=Ng)
165 # Set up plane-wave descriptor
166 qd = KPointDescriptor(np.array([[0., 0., 0.]]))
167 pd = PWDescriptor(relative_ecut / Ha / a**2., gd, complex, qd)
169 # Calculate the plane-wave components analytically
170 ntest_G = ae_1s_density_plane_waves(pd, R_v, a=a)
172 # Calculate the pseudo and ae densities on the radial grid
173 n_g = ae_1s_density(rgd.r_g, a=a)
174 gcut = rgd.floor(rcut)
175 nt_g, _ = rgd.pseudize(n_g, gcut)
177 # Set up pseudo and ae densities on the Lebedev quadrature
178 from gpaw.sphere.lebedev import Y_nL
179 Y_nL = Y_nL[:, :9] # include only s, p and d
180 nL = Y_nL.shape[1]
181 n_sLg = np.zeros((1, nL, Ng), dtype=float)
182 nt_sLg = np.zeros((1, nL, Ng), dtype=float)
183 # 1s <=> (l,m) = (0,0) <=> L = 0
184 n_sLg[0, 0, :] += np.sqrt(4. * np.pi) * n_g # Y_0 = 1 / sqrt(4pi)
185 nt_sLg[0, 0, :] += np.sqrt(4. * np.pi) * nt_g
187 # Calculate the pseudo density on the real-space grid
188 # ------------------------------------------------- #
189 # Generate splines on for the pseudo density on the radial grid
190 spline = rgd.spline(nt_g, l=0, rcut=redge)
191 # Use the LocalizedFunctionsCollection to generate pseudo density
192 # on the cubic real space grid
193 nt_R = gd.zeros()
194 lfc = LFC(gd, [[spline]])
195 lfc.set_positions(pos_ac)
196 lfc.add(nt_R, c_axi=np.sqrt(4. * np.pi)) # Y_0 = 1 / sqrt(4pi)
197 nt_sR = np.array([nt_R])
199 # Create MicroSetup(s)
200 micro_setup = MicroSetup(rgd, Y_nL, n_sLg, nt_sLg)
201 micro_setups = [micro_setup]
203 gs = EmptyGSAdapter() # hack to pass isinstance in constructor
204 context = ResponseContext()
205 for rshe_params in rshe_params_p:
206 # Initialize the LocalPAWFTCalculator with an empty gs adapter
207 localft_calc = LocalFTCalculator.from_rshe_parameters(gs, context,
208 **rshe_params)
210 # Compute the plane-wave components numerically
211 n_G = localft_calc.engine.calculate(pd, nt_sR, [R_v], micro_setups,
212 add_total_density)
214 # Check validity of numerical results
215 assert np.allclose(n_G, ntest_G, rtol=rtol)
218@pytest.mark.response
219def test_Fe_bxc(gpw_files):
220 """Test the symmetry relation
222 W_xc^z(G)^* = W_xc^z(-G)
224 for a real life system with d-electrons (bcc-Fe)."""
225 # ---------- Inputs ---------- #
227 # Bxc calculation
228 ecut = 100
230 # ---------- Script ---------- #
232 # Bxc calculation
234 # Set up calculator and plane-wave descriptor
235 calc = GPAW(gpw_files['fe_pw'], parallel=dict(domain=1))
236 atoms = calc.atoms
237 gs = ResponseGroundStateAdapter(calc)
238 context = ResponseContext()
239 localft_calc = LocalFTCalculator.from_rshe_parameters(gs, context)
240 pd0 = get_pw_descriptor(atoms, calc, [0., 0., 0.],
241 ecut=ecut,
242 gammacentered=True)
244 Wxc_G = localft_calc(pd0, add_LSDA_Wxc)
246 # Part 3: Check symmetry relation
247 G1_G, G2_G = get_inversion_pairs(pd0)
249 assert np.allclose(np.conj(Wxc_G[G1_G]), Wxc_G[G2_G])
252# ---------- Test functionality ---------- #
255class EmptyGSAdapter(ResponseGroundStateAdapter):
256 # Make an empty subclass to pass isinstance in constructor
257 # In a future where the response code has been liberated from GPAW
258 # calculator objects, the
259 # >>> assert isinstance(gs, ResponseGroundStateAdapter)
260 # statements can be deleted, making this class redundant.
262 def __init__(self):
263 pass
266def get_inversion_pairs(pd0):
267 """Get all pairs of G-indices which correspond to inverted reciprocal
268 lattice vectors G and -G."""
269 G_Gc = get_pw_coordinates(pd0)
271 G1_G = []
272 G2_G = []
273 paired_indices = []
274 for G1, G1_c in enumerate(G_Gc):
275 if G1 in paired_indices:
276 continue # Already paired
278 for G2, G2_c in enumerate(G_Gc):
279 if np.all(G2_c == -G1_c):
280 G1_G.append(G1)
281 G2_G.append(G2)
282 paired_indices += [G1, G2]
283 break
285 assert len(np.unique(paired_indices)) == len(G_Gc)
287 return G1_G, G2_G