Coverage for gpaw/test/response/test_mft.py: 92%
204 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1"""Calculate the Heisenberg exchange constants in Fe and Co using the MFT.
2Test with unrealisticly loose parameters to catch if the numerics change.
3"""
5# General modules
6import pytest
7import numpy as np
9# Script modules
10from gpaw import GPAW
12from gpaw.response import ResponseGroundStateAdapter, ResponseContext
13from gpaw.response.chiks import ChiKSCalculator
14from gpaw.response.localft import LocalFTCalculator, LocalPAWFTCalculator
15from gpaw.response.site_data import (AtomicSites,
16 calculate_site_magnetization,
17 calculate_site_zeeman_energy)
18from gpaw.response.mft import (IsotropicExchangeCalculator,
19 HeisenbergExchangeCalculator,
20 calculate_single_particle_site_magnetization,
21 calculate_pair_site_magnetization,
22 calculate_single_particle_site_zeeman_energy,
23 calculate_pair_site_zeeman_energy,
24 calculate_exchange_parameters)
25from gpaw.response.site_kernels import (SphericalSiteKernels,
26 CylindricalSiteKernels,
27 ParallelepipedicSiteKernels)
28from gpaw.response.heisenberg import (calculate_single_site_magnon_energies,
29 calculate_fm_magnon_energies)
30from gpaw.test.gpwfile import response_band_cutoff
31from gpaw.test.response.test_chiks import generate_qrel_q, get_q_c
34@pytest.mark.response
35@pytest.mark.kspair
36def test_Fe_bcc(in_tmp_dir, gpw_files):
37 # ---------- Inputs ---------- #
39 # MFT calculation
40 ecut = 50
41 # Do the high symmetry points of the bcc lattice
42 q_qc = np.array([[0, 0, 0], # Gamma
43 [0.5, -0.5, 0.5], # H
44 [0.0, 0.0, 0.5], # N
45 ])
46 # Define site kernels to test
47 # Test a single site of spherical and cylindrical geometries
48 rc_pa = np.array([[1.0], [1.5], [2.0]])
49 hc_pa = np.array([[1.0], [1.5], [2.0]])
50 ez_pav = np.array([[[1., 0., 0.]], [[0., 1., 0.]], [[0., 0., 1.]]])
52 # ---------- Script ---------- #
54 # Extract the ground state fixture
55 calc = GPAW(gpw_files['fe_pw'], parallel=dict(domain=1))
56 nbands = response_band_cutoff['fe_pw']
57 atoms = calc.atoms
59 # Set up site kernels with a single site
60 positions = atoms.get_positions()
61 sitekernels = SphericalSiteKernels(positions, rc_pa)
62 sitekernels.append(CylindricalSiteKernels(positions, ez_pav,
63 rc_pa, hc_pa))
64 # Set up a kernel to fill out the entire unit cell
65 sitekernels.append(ParallelepipedicSiteKernels(positions,
66 [[atoms.get_cell()]]))
68 # Initialize the Heisenberg exchange calculator
69 gs = ResponseGroundStateAdapter(calc)
70 context = ResponseContext()
71 chiks_calc = ChiKSCalculator(gs, context,
72 ecut=ecut, nbands=nbands, gammacentered=True)
73 localft_calc = LocalFTCalculator.from_rshe_parameters(gs, context)
74 isoexch_calc = IsotropicExchangeCalculator(chiks_calc, localft_calc)
76 # Allocate array for the exchange constants
77 nq = len(q_qc)
78 nsites = sitekernels.nsites
79 npartitions = sitekernels.npartitions
80 J_qabp = np.empty((nq, nsites, nsites, npartitions), dtype=complex)
82 # Calcualate the exchange constant for each q-point
83 for q, q_c in enumerate(q_qc):
84 J_qabp[q] = isoexch_calc(q_c, sitekernels)
86 # Since we only have a single site, reduce the array
87 J_qp = J_qabp[:, 0, 0, :]
89 # Calculate the magnon energies
90 mm = 2.21
91 mm_ap = mm * np.ones((1, npartitions)) # Magnetic moments
92 mw_qp = calculate_fm_magnon_energies(J_qabp, q_qc, mm_ap)[:, 0, :]
94 # Compare results to test values
95 test_J_pq = np.array(
96 [[2.1907596825086455, 1.172424411323134, 1.6060583789867644],
97 [2.612428039019977, 1.2193926800088601, 1.7635196888465006],
98 [6.782367391186284, 0.2993922109834177, 1.9346016211386057],
99 [1.5764800860123762, 0.8365204592352894, 1.1648584638500161],
100 [2.4230224513213234, 1.2179759558303274, 1.6691805687218078],
101 [5.35668502504496, 0.3801778545994659, 1.6948968244858478],
102 [2.523580017606111, 1.21779750159267, 1.7637120466695273]])
103 test_mw_pq = np.array(
104 [[0.0, 0.9215703811633589, 0.5291414511510236],
105 [0.0, 1.2606654832679791, 0.7682428508357253],
106 [0.0, 5.866945864436984, 4.38711834393455],
107 [0.0, 0.6696467210652369, 0.3725082553505521],
108 [0.0, 1.0905398149239784, 0.682209848506349],
109 [0.0, 4.503626398593207, 3.313835475619106],
110 [0.0, 1.181703634401304, 0.6876633221145555]])
112 # Exchange constants
113 assert J_qp.imag == pytest.approx(0.0)
114 assert J_qp.T.real == pytest.approx(test_J_pq, rel=2e-3)
116 # Magnon energies
117 assert mw_qp.T == pytest.approx(test_mw_pq, rel=2e-3)
120@pytest.mark.response
121@pytest.mark.kspair
122def test_Co_hcp(in_tmp_dir, gpw_files):
123 # ---------- Inputs ---------- #
125 # MFT calculation
126 ecut = 100
127 # Do high symmetry points of the hcp lattice
128 q_qc = np.array([[0, 0, 0], # Gamma
129 [0.5, 0., 0.], # M
130 [0., 0., 0.5] # A
131 ])
133 # Use spherical site kernels in a radius range which should yield
134 # stable results
135 rc_pa = np.array([[1.0, 1.0], [1.1, 1.1], [1.2, 1.2]])
137 # Unfortunately, the usage of symmetry leads to such extensive repetition
138 # of random noise, that one cannot trust individual values of J very well.
139 # This is improved when increasing the number of k-points, but the problem
140 # never completely vanishes
141 J_atol = 1.e-2
142 J_rtol = 5.e-2
143 # However, derived physical values have an increased error cancellation due
144 # to their collective nature.
145 mw_rtol = 25e-3 # relative tolerance of absolute results
146 mw_ctol = 5.e-2 # relative tolerance on kernel and eta self-consistency
148 # ---------- Script ---------- #
150 # Extract the ground state fixture
151 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1))
152 nbands = response_band_cutoff['co_pw']
153 atoms = calc.get_atoms()
155 # Set up spherical site kernels
156 positions = atoms.get_positions()
157 sitekernels = SphericalSiteKernels(positions, rc_pa)
159 # Set up a site kernel to fill out the entire unit cell
160 cell_cv = atoms.get_cell()
161 cc_v = np.sum(cell_cv, axis=0) / 2. # Unit cell center
162 ucsitekernels = ParallelepipedicSiteKernels([cc_v], [[cell_cv]])
164 # Initialize the exchange calculator with and without symmetry
165 gs = ResponseGroundStateAdapter(calc)
166 context = ResponseContext()
167 chiks_calc0 = ChiKSCalculator(gs, context, qsymmetry=False,
168 ecut=ecut, nbands=nbands, gammacentered=True)
169 localft_calc = LocalPAWFTCalculator(gs, context)
170 isoexch_calc0 = IsotropicExchangeCalculator(chiks_calc0, localft_calc)
171 chiks_calc1 = ChiKSCalculator(gs, context,
172 ecut=ecut, nbands=nbands, gammacentered=True)
173 isoexch_calc1 = IsotropicExchangeCalculator(chiks_calc1, localft_calc)
175 # Allocate array for the spherical site exchange constants
176 nq = len(q_qc)
177 nsites = sitekernels.nsites
178 npartitions = sitekernels.npartitions
179 J_qabp = np.empty((nq, nsites, nsites, npartitions), dtype=complex)
181 # Allocate array for the unit cell site exchange constants
182 Juc_qs = np.empty((nq, 2), dtype=complex)
184 # Calcualate the exchange constants for each q-point
185 for q, q_c in enumerate(q_qc):
186 J_qabp[q] = isoexch_calc0(q_c, sitekernels)
187 chiksr_buffer = isoexch_calc0._chiksr
188 Juc_qs[q, 0] = isoexch_calc0(q_c, ucsitekernels)[0, 0, 0]
189 assert isoexch_calc0._chiksr is chiksr_buffer, \
190 'Two subsequent IsotropicExchangeCalculator calls with the same '\
191 'q_c, should reuse, not update, the chiks buffer'
193 Juc_qs[q, 1] = isoexch_calc1(q_c, ucsitekernels)[0, 0, 0]
195 # Calculate the magnon energy
196 mom = atoms.get_magnetic_moment()
197 mm_ap = mom / 2.0 * np.ones((nsites, npartitions))
198 mw_qnp = calculate_fm_magnon_energies(J_qabp, q_qc, mm_ap)
199 mw_qnp = np.sort(mw_qnp, axis=1) # Make sure the eigenvalues are sorted
200 mwuc_qs = calculate_single_site_magnon_energies(Juc_qs, q_qc, mom)
202 # Compare results to test values
203 print(J_qabp[..., 1], mw_qnp[..., 1], mwuc_qs[:, 0])
204 test_J_qab = np.array([[[1.23106207 - 0.j, 0.25816335 - 0.j],
205 [0.25816335 + 0.j, 1.23106207 + 0.j]],
206 [[0.88823839 + 0.j, 0.07345416 - 0.04947835j],
207 [0.07345416 + 0.04947835j, 0.88823839 + 0.j]],
208 [[1.09349955 - 0.j, 0.00000010 - 0.01176761j],
209 [0.00000010 + 0.01176761j, 1.09349955 - 0.j]]])
210 test_mw_qn = np.array([[0., 0.64793939],
211 [0.64304039, 0.86531921],
212 [0.48182997, 0.51136436]])
213 test_mwuc_q = np.array([0., 0.69678659, 0.44825874])
215 # Exchange constants
216 assert J_qabp[..., 1] == pytest.approx(test_J_qab, abs=J_atol, rel=J_rtol)
218 # Magnon energies
219 assert np.all(np.abs(mw_qnp[0, 0, :]) < 1.e-8) # Goldstone theorem
220 assert np.allclose(mwuc_qs[0, :], 0.) # Goldstone
221 assert mw_qnp[1:, 0, 1] == pytest.approx(test_mw_qn[1:, 0], rel=mw_rtol)
222 assert mw_qnp[:, 1, 1] == pytest.approx(test_mw_qn[:, 1], rel=mw_rtol)
223 assert mwuc_qs[1:, 0] == pytest.approx(test_mwuc_q[1:], rel=mw_rtol)
225 # Check self-consistency of results
226 # We should be in a radius range, where the magnon energies don't change
227 assert np.allclose(mw_qnp[1:, 0, ::2],
228 test_mw_qn[1:, 0, np.newaxis], rtol=mw_ctol)
229 assert np.allclose(mw_qnp[:, 1, ::2],
230 test_mw_qn[:, 1, np.newaxis], rtol=mw_ctol)
231 # Check that symmetry toggle do not change the magnon energies
232 assert np.allclose(mwuc_qs[1:, 0], mwuc_qs[1:, 1], rtol=mw_ctol)
235@pytest.mark.response
236@pytest.mark.kspair
237@pytest.mark.parametrize('qrel', generate_qrel_q())
238def test_Co_site_magnetization_sum_rule(in_tmp_dir, gpw_files, qrel):
239 # Set up ground state adapter and basic parameters
240 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1))
241 gs = ResponseGroundStateAdapter(calc)
242 sites = get_co_sites(gs)
243 context = 'Co_sum_rule.txt'
244 nbands = response_band_cutoff['co_pw']
246 # Get wave vector to test
247 q_c = get_q_c('co_pw', qrel)
249 # Calculate site magnetization
250 magmom_ar = calculate_site_magnetization(gs, sites)
252 # ----- Single-particle site magnetization ----- #
254 # Test that the single-particle site magnetization matches a conventional
255 # calculation based on the density
256 sp_magmom_ar = calculate_single_particle_site_magnetization(
257 gs, sites, context=context)
258 assert sp_magmom_ar == pytest.approx(magmom_ar, rel=5e-3)
260 # ----- Two-particle site magnetization ----- #
262 magmom_abr = calculate_pair_site_magnetization(
263 gs, sites, context=context, q_c=q_c, nbands=nbands)
265 # Test that the site pair magnetization is a positive-valued diagonal
266 # real array
267 tp_magmom_ra = magmom_abr.diagonal()
268 assert np.all(tp_magmom_ra.real > 0)
269 assert np.all(np.abs(tp_magmom_ra.imag) / tp_magmom_ra.real < 1e-6)
270 assert np.all(np.abs(np.diagonal(np.fliplr( # off-diagonal elements
271 magmom_abr))) / tp_magmom_ra.real < 5e-2)
273 # Test that the magnetic moments on the two Co atoms are identical
274 tp_magmom_ar = magmom_abr.diagonal().T.real
275 assert tp_magmom_ar[0] == pytest.approx(tp_magmom_ar[1], rel=1e-4)
277 # Test that the result more or less matches a conventional calculation at
278 # close-packing
279 assert np.average(tp_magmom_ar, axis=0)[-1] == pytest.approx(
280 np.average(magmom_ar, axis=0)[-1], rel=5e-2)
282 # Test values against reference
283 print(np.average(tp_magmom_ar, axis=0)[::2])
284 assert np.average(tp_magmom_ar, axis=0)[::2] == pytest.approx(
285 np.array([3.91823444e-04, 1.45641911e-01, 6.85939109e-01,
286 1.18813171e+00, 1.49761591e+00, 1.58954270e+00]), rel=5e-2)
288 # import matplotlib.pyplot as plt
289 # from ase.units import Bohr
290 # rc_r = sites.rc_ap[0] * Bohr
291 # plt.plot(rc_r, magmom_ar[0], '-o', mec='k')
292 # plt.plot(rc_r, sp_magmom_ar[0], '-o', mec='k', zorder=0)
293 # plt.plot(rc_r, tp_magmom_ar[0], '-o', mec='k', zorder=1)
294 # plt.xlabel(r'$r_\mathrm{c}$ [$\mathrm{\AA}$]')
295 # plt.ylabel(r'$m$ [$\mu_\mathrm{B}$]')
296 # plt.title(str(q_c))
297 # plt.show()
300@pytest.mark.response
301@pytest.mark.kspair
302@pytest.mark.parametrize('qrel', generate_qrel_q())
303def test_Co_site_zeeman_energy_sum_rule(in_tmp_dir, gpw_files, qrel):
304 # Set up ground state adapter and atomic site data
305 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1))
306 gs = ResponseGroundStateAdapter(calc)
307 sites = get_co_sites(gs)
308 context = ResponseContext('Co_sum_rule.txt')
309 nbands = response_band_cutoff['co_pw']
311 # Get wave vector to test
312 q_c = get_q_c('co_pw', qrel)
314 # Calculate the site Zeeman energy
315 EZ_ar = calculate_site_zeeman_energy(gs, sites)
317 # ----- Single-particle site Zeeman energy ----- #
319 # Test that the results match a conventional calculation
320 sp_EZ_ar = calculate_single_particle_site_zeeman_energy(
321 gs, sites, context=context)
322 assert sp_EZ_ar == pytest.approx(EZ_ar, rel=5e-3)
324 # ----- Two-particle site Zeeman energy ----- #
326 EZ_abr = calculate_pair_site_zeeman_energy(
327 gs, sites, context=context, q_c=q_c, nbands=nbands)
329 # Test that the pair site Zeeman energy is a positive-valued diagonal
330 # real array
331 tp_EZ_ra = EZ_abr.diagonal()
332 assert np.all(tp_EZ_ra.real > 0)
333 assert np.all(np.abs(tp_EZ_ra.imag) / tp_EZ_ra.real < 1e-4)
334 assert np.all(np.abs(np.diagonal(np.fliplr( # off-diagonal elements
335 EZ_abr))) / tp_EZ_ra.real < 5e-2)
337 # Test that the Zeeman energy on the two Co atoms is identical
338 tp_EZ_ar = EZ_abr.diagonal().T.real
339 assert tp_EZ_ar[0] == pytest.approx(tp_EZ_ar[1], rel=1e-4)
341 # Test values against reference
342 print(np.average(tp_EZ_ar, axis=0)[::2])
343 assert np.average(tp_EZ_ar, axis=0)[::2] * 2. == pytest.approx(
344 np.array([3.68344584e-04, 3.13780575e-01, 1.35409600e+00,
345 2.14237563e+00, 2.52032513e+00, 2.61406726e+00]), rel=5e-2)
347 # Test ability to distribute band and spin transitions
348 if context.comm.size > 1:
349 assert calculate_pair_site_zeeman_energy(
350 gs, sites, context=context,
351 q_c=q_c, nbands=nbands, nblocks='max') == pytest.approx(EZ_abr)
353 # import matplotlib.pyplot as plt
354 # from ase.units import Bohr
355 # rc_r = sites.rc_ap[0] * Bohr
356 # plt.plot(rc_r, EZ_ar[0], '-o', mec='k')
357 # plt.plot(rc_r, sp_EZ_ar[0], '-o', mec='k', zorder=0)
358 # plt.plot(rc_r, tp_EZ_ar[0], '-o', mec='k', zorder=1)
359 # plt.xlabel(r'$r_\mathrm{c}$ [$\mathrm{\AA}$]')
360 # plt.ylabel(r'$E_\mathrm{Z}$ [eV]')
361 # plt.title(str(q_c))
362 # plt.show()
365def get_Co_exchange_reference(qrel):
366 if qrel == 0.:
367 return np.array([0.186, 0.731, 1.048, 1.14, 1.145])
368 elif qrel == 0.25:
369 return np.array([0.166, 0.656, 0.943, 1.029, 1.035])
370 elif qrel == 0.5:
371 return np.array([0.151, 0.594, 0.856, 0.936, 0.943])
372 raise ValueError(qrel)
375@pytest.mark.response
376@pytest.mark.kspair
377@pytest.mark.parametrize('qrel', generate_qrel_q())
378def test_Co_exchange(in_tmp_dir, gpw_files, qrel):
379 # Set up ground state adapter and atomic site data
380 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1))
381 gs = ResponseGroundStateAdapter(calc)
382 context = ResponseContext('Co_exchange.txt')
383 sites = get_co_sites(gs)
384 nbands = response_band_cutoff['co_pw']
386 # Get wave vector to test
387 q_c = get_q_c('co_pw', qrel)
389 # Calculate the exchange parameters
390 J_abr = calculate_exchange_parameters(
391 gs, sites, q_c, context=context, nbands=nbands, nblocks=1)
393 # Test that J is hermitian in (a,b)
394 assert np.transpose(J_abr, (1, 0, 2)).conj() == pytest.approx(J_abr)
395 # Test the Co site symmetry
396 J_ar = J_abr.diagonal().T
397 assert J_ar[0] == pytest.approx(J_ar[1], rel=1e-4)
398 # Test against reference values
399 Jref_r = get_Co_exchange_reference(qrel)
400 assert J_ar[0, 2::2] == pytest.approx(Jref_r, abs=2e-3)
402 # Test ability to distribute band and spin transitions
403 if context.comm.size > 1:
404 assert calculate_exchange_parameters(
405 gs, sites, q_c, context=context, nbands=nbands,
406 nblocks='max') == pytest.approx(J_abr)
408 # Calculate the magnon energy at the Gamma-point
409 if qrel == 0.:
410 mom = calc.get_atoms().get_magnetic_moment()
411 mm_ar = mom / 2. * np.ones(J_ar.shape)
412 mw_nr = calculate_fm_magnon_energies(
413 np.array([J_abr]), np.array([q_c]), mm_ar)[0]
414 mw_nr = np.sort(mw_nr, axis=0) # Make sure the eigenvalues are sorted
416 assert mw_nr[0] == pytest.approx(0.) # Goldstone theorem
417 assert mw_nr[1, 2::2] == pytest.approx(np.array(
418 [0.092, 0.358, 0.504, 0.539, 0.539]), abs=1e-3)
420 # import matplotlib.pyplot as plt
421 # from ase.units import Bohr
422 # rc_r = sites.rc_ap[0] * Bohr
423 # for n, mw_r in enumerate(mw_nr):
424 # plt.plot(rc_r, mw_r, '-o', mec='k', label=str(n))
425 # plt.legend()
426 # plt.xlabel(r'$r_\mathrm{c}$ [$\mathrm{\AA}$]')
427 # plt.ylabel(r'$\hbar\omega$ [eV]')
428 # plt.title(str(q_c))
429 # plt.show()
432@pytest.mark.response
433@pytest.mark.kspair
434@pytest.mark.parallel
435def test_heisenberg_distribution_over_transitions(in_tmp_dir, gpw_files):
436 # Set up ground state adapter and atomic site data
437 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1))
438 gs = ResponseGroundStateAdapter(calc)
439 sites = get_co_sites(gs)
441 # To properly test the distributions over transitions, we need a value of
442 # nbands, which produces a number of band and spin transitions, which isn't
443 # divisible by the number of blocks.
444 nbands = response_band_cutoff['co_pw'] - 3
445 context = ResponseContext('distributed.txt')
446 calc = HeisenbergExchangeCalculator(
447 gs, sites, context=context, nbands=nbands, nblocks='max')
448 assert context.comm.size % 2 == 0
449 assert len(calc.transitions) % 2 > 0
450 J_abr = calc(q_c=[0, 0, 0]).array
452 # Test that the result doesn't depend on nblocks
453 context.new_txt_and_timer('undistributed.txt')
454 ref_calc = HeisenbergExchangeCalculator(
455 gs, sites, context=context, nbands=nbands, nblocks=1)
456 assert J_abr == pytest.approx(ref_calc(q_c=[0, 0, 0]).array)
457 assert np.all(np.abs(J_abr) > 1e-10) # check that some actual work is done
460# ---------- Test functionality ---------- #
463def get_co_sites(gs):
464 from gpaw.response.site_data import get_site_radii_range
465 # Set up site radii
466 rmin_a, _ = get_site_radii_range(gs)
467 # Make sure that the two sites do not overlap
468 nn_dist = min(2.5071, np.sqrt(2.5071**2 / 3 + 4.0695**2 / 4))
469 rc_r = np.linspace(rmin_a[0], nn_dist / 2, 11)
470 return AtomicSites(indices=[0, 1], radii=[rc_r, rc_r])