Coverage for gpaw/test/response/test_site_data.py: 100%
108 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import pytest
3import numpy as np
5from ase.spacegroup import crystal
6from ase.units import Bohr
8from gpaw import GPAW
9from gpaw.sphere.integrate import integrate_lebedev
11from gpaw.response import ResponseGroundStateAdapter
12from gpaw.response.site_data import (AtomicSites, AtomicSiteData,
13 calculate_site_magnetization,
14 calculate_site_zeeman_energy,
15 get_site_radii_range,
16 maximize_site_magnetization)
17from gpaw.response.localft import add_spin_polarization
20@pytest.mark.response
21def test_Fe_site_magnetization(gpw_files):
22 # Set up ground state adapter
23 calc = GPAW(gpw_files['fe_pw'], parallel=dict(domain=1))
24 gs = ResponseGroundStateAdapter(calc)
26 # Extract valid site radii range
27 rmin_a, rmax_a = get_site_radii_range(gs)
28 rmin = rmin_a[0] # Only one magnetic atom in the unit cell
29 rmax = rmax_a[0]
30 # We expect rmax to be equal to the nearest neighbour distance
31 # subtracted with the augmentation sphere radius. For a bcc lattice,
32 # nn_dist = sqrt(3) a / 2:
33 augr = gs.get_aug_radii()[0]
34 rmax_expected = np.sqrt(3) * 2.867 / 2. - augr * Bohr
35 assert abs(rmax - rmax_expected) < 1e-6
36 # Test that an error is raised outside the valid range
37 with pytest.raises(AssertionError):
38 AtomicSiteData(
39 gs, AtomicSites(indices=[0], # Too small radii
40 radii=[np.linspace(rmin * 0.8, rmin, 5)]))
41 with pytest.raises(AssertionError):
42 AtomicSiteData(
43 gs, AtomicSites(indices=[0], # Too large radii
44 radii=[np.linspace(rmax, rmax * 1.2, 5)]))
45 # Define atomic sites to span the valid range
46 rc_r = np.linspace(rmin_a[0], rmax_a[0], 100)
47 # Add the radius of the augmentation sphere explicitly
48 rc_r = np.append(rc_r, [augr * Bohr])
49 sites = AtomicSites(indices=[0], radii=[rc_r])
50 site_data = AtomicSiteData(gs, sites)
52 # Calculate site magnetization
53 magmom_ar = site_data.calculate_magnetic_moments()
54 magmom_r = magmom_ar[0]
56 # Test that a cutoff at the augmentation sphere radius reproduces
57 # the local magnetic moment of the GPAW calculation
58 magmom_at_augr = calc.get_atoms().get_magnetic_moments()[0]
59 assert abs(magmom_r[-1] - magmom_at_augr) < 4e-2
61 # Do a manual calculation of the magnetic moment using the
62 # all-electron partial waves
63 # Calculate all-electron m(r)
64 micro_setup = site_data.micro_setup_a[0]
65 m_ng = np.array([micro_setup.rgd.zeros()
66 for n in range(micro_setup.Y_nL.shape[0])])
67 for n, Y_L in enumerate(micro_setup.Y_nL):
68 n_sg = np.dot(Y_L, micro_setup.n_sLg)
69 add_spin_polarization(micro_setup.rgd, n_sg, m_ng[n, :])
70 # Integrate with varrying radii
71 m_g = integrate_lebedev(m_ng)
72 ae_magmom_r = np.array([
73 micro_setup.rgd.integrate_trapz(m_g, rcut=rcut / Bohr)
74 for rcut in rc_r])
75 # Test that values match approximately inside the augmentation sphere
76 inaug_r = rc_r <= augr * Bohr
77 assert magmom_r[inaug_r] == pytest.approx(ae_magmom_r[inaug_r], abs=3e-2)
79 # import matplotlib.pyplot as plt
80 # plt.plot(rc_r[:-1], magmom_r[:-1])
81 # plt.plot(rc_r[:-1], ae_magmom_r[:-1], zorder=0)
82 # plt.axvline(augr * Bohr, c='0.5', linestyle='--')
83 # plt.xlabel(r'$r_\mathrm{c}$ [$\mathrm{\AA}$]')
84 # plt.ylabel(r'$m$ [$\mu_\mathrm{B}$]')
85 # plt.show()
88@pytest.mark.response
89def test_Co_site_data(gpw_files):
90 # Set up ground state adapter
91 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1))
92 gs = ResponseGroundStateAdapter(calc)
94 # Extract valid site radii range
95 rmin_a, rmax_a = get_site_radii_range(gs)
96 # The valid ranges should be equal due to symmetry
97 assert abs(rmin_a[1] - rmin_a[0]) < 1e-8
98 assert abs(rmax_a[1] - rmax_a[0]) < 1e-8
99 rmin = rmin_a[0]
100 rmax = rmax_a[0]
101 # We expect rmax to be equal to the nearest neighbour distance
102 # subtracted with the augmentation sphere radius. For the hcp-lattice,
103 # nn_dist = min(a, sqrt(a^2/3 + c^2/4)):
104 augr_a = gs.get_aug_radii()
105 assert abs(augr_a[1] - augr_a[0]) < 1e-8
106 augr = augr_a[0]
107 rmax_expected = min(2.5071, np.sqrt(2.5071**2 / 3 + 4.0695**2 / 4))
108 rmax_expected -= augr * Bohr
109 assert abs(rmax - rmax_expected) < 1e-6
111 # Use radii spanning the entire valid range
112 rc_r = np.linspace(rmin, rmax, 101)
113 # Add the radius of the augmentation sphere explicitly
114 rc_r = np.append(rc_r, [augr * Bohr])
115 nr = len(rc_r)
116 # Varry the site radii together and independently
117 rc1_r = list(rc_r) + list(rc_r) + [augr * Bohr] * nr
118 rc2_r = list(rc_r) + [augr * Bohr] * nr + list(rc_r)
119 sites = AtomicSites(indices=[0, 1], radii=[rc1_r, rc2_r])
121 # Calculate site magnetization
122 magmom_ar = calculate_site_magnetization(gs, sites)
124 # Test that the magnetization inside the augmentation sphere matches
125 # the local magnetic moment of the GPAW calculation
126 magmom_at_augr_a = calc.get_atoms().get_magnetic_moments()
127 assert magmom_ar[:, -1] == pytest.approx(magmom_at_augr_a, abs=2e-2)
129 # Test consistency of varrying radii
130 assert magmom_ar[0, :nr] == pytest.approx(magmom_ar[1, :nr])
131 assert magmom_ar[0, nr:2 * nr] == pytest.approx(magmom_ar[0, :nr])
132 assert magmom_ar[0, 2 * nr:] == pytest.approx([magmom_ar[0, -1]] * nr)
133 assert magmom_ar[1, nr:2 * nr] == pytest.approx([magmom_ar[1, -1]] * nr)
134 assert magmom_ar[1, 2 * nr:] == pytest.approx(magmom_ar[1, :nr])
136 # Calculate the maximized site magnetization
137 rm_a, mm_a = maximize_site_magnetization(gs)
138 # Test radius consistency
139 assert rm_a[0] == pytest.approx(rm_a[1]) # Co site symmetry
140 assert np.average(rm_a) == pytest.approx(1.133357) # reference value
141 # Test moment consistency
142 assert mm_a[0] == pytest.approx(mm_a[1]) # Co site symmetry
143 assert np.average(mm_a) == pytest.approx(1.6362) # reference value
144 assert np.max(magmom_ar) < np.average(mm_a) < np.max(magmom_ar) * 1.01
146 # Calculate the atomic Zeeman energy
147 rc_r = rc_r[:-1]
148 sites = AtomicSites(indices=[0, 1], radii=[rc_r, rc_r])
149 EZ_ar = calculate_site_zeeman_energy(gs, sites)
150 print(EZ_ar[0, ::20])
152 # Test that the Zeeman energy comes out as expected
153 assert EZ_ar[0] == pytest.approx(EZ_ar[1])
154 assert EZ_ar[0, ::20] * 2 == pytest.approx([0.02638351, 1.41476112,
155 2.49540004, 2.79727200,
156 2.82727948, 2.83670767],
157 rel=1e-3)
159 # import matplotlib.pyplot as plt
160 # plt.subplot(1, 2, 1)
161 # plt.plot(rc_r, magmom_ar[0, :nr - 1])
162 # plt.axvline(np.average(rm_a), linestyle=':')
163 # plt.axvline(augr * Bohr, c='0.5', linestyle='--')
164 # plt.xlabel(r'$r_\mathrm{c}$ [$\mathrm{\AA}$]')
165 # plt.ylabel(r'$m$ [$\mu_\mathrm{B}$]')
166 # plt.subplot(1, 2, 2)
167 # plt.plot(rc_r, EZ_ar[0])
168 # plt.axvline(np.average(rm_a), linestyle=':')
169 # plt.axvline(augr * Bohr, c='0.5', linestyle='--')
170 # plt.xlabel(r'$r_\mathrm{c}$ [$\mathrm{\AA}$]')
171 # plt.ylabel(r'$E_\mathrm{Z}$ [eV]')
172 # plt.show()
175@pytest.mark.response
176def test_valid_site_radii_symmetry():
177 # Set up Cr2O3 crystal
178 cellpar = [4.95721, 4.95721, 13.59170, 90, 90, 120]
179 Cr_c = [0, 0, 0.34734]
180 O_c = [0.30569, 0.0, 0.25]
181 spos_ac = [Cr_c, O_c]
182 atoms = crystal('CrO',
183 spacegroup=167,
184 cellpar=cellpar,
185 basis=spos_ac,
186 primitive_cell=True,
187 pbc=True)
188 # from ase.visualize import view
189 # view(atoms)
191 # Set up calculator with a specific grid spacing and generate adapter
192 spacing = 0.1
193 calc = GPAW(mode='fd', h=spacing)
194 calc.initialize(atoms)
195 gs = DummyAdapter(calc)
197 # Generate valid site radii range
198 rmin_A, rmax_A = get_site_radii_range(gs)
199 # Test that the minimum radius corresponds loosely to the specified
200 # spacing. The correspondance would be exact for a cubic cell.
201 assert rmin_A == pytest.approx(np.ones(len(atoms)) * spacing / 2, rel=0.2)
202 # Test that all Cr and O atoms result in symmetrically equivalent maximum
203 # cutoff radii
204 CrO_dist = 1.966
205 Cr_aug_radius = 2.3 * Bohr
206 O_aug_radius = 1.3 * Bohr
207 is_Cr = np.array([c == 'Cr' for c in atoms.get_chemical_symbols()])
208 refmax_A = np.empty(len(atoms))
209 refmax_A[is_Cr] = CrO_dist - O_aug_radius
210 refmax_A[~is_Cr] = CrO_dist - Cr_aug_radius
211 assert rmax_A == pytest.approx(refmax_A, abs=0.001)
214class DummyAdapter(ResponseGroundStateAdapter):
215 def __init__(self, calc):
216 from gpaw.response.groundstate import PAWDatasetCollection
217 self.atoms = calc.atoms
218 self.gd = calc.wfs.gd
219 self.pawdatasets = PAWDatasetCollection(calc.setups)