Coverage for gpaw/test/response/test_cobalt_sf_gssALDA.py: 100%
131 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1# General modules
2import pytest
3import numpy as np
5# Script modules
6from gpaw import GPAW
7from gpaw.test import findpeak
8from gpaw.response import ResponseGroundStateAdapter, ResponseContext
9from gpaw.response.chiks import ChiKSCalculator
10from gpaw.response.susceptibility import (ChiFactory, spectral_decomposition,
11 EigendecomposedSpectrum,
12 read_full_spectral_weight,
13 read_eigenmode_lineshapes)
14from gpaw.response.fxc_kernels import AdiabaticFXCCalculator
15from gpaw.response.goldstone import FMGoldstoneScaling
16from gpaw.response.pair_functions import read_susceptibility_array
17from gpaw.test.gpwfile import response_band_cutoff
20@pytest.mark.kspair
21@pytest.mark.response
22def test_response_cobalt_sf_gssALDA(in_tmp_dir, gpw_files):
23 # ---------- Inputs ---------- #
25 fxc = 'ALDA'
26 q_qc = [[0.0, 0.0, 0.0], [1. / 4., 0.0, 0.0]] # Two q-points along G-M
27 # We make sure to have exactly 49 frequency points, so that one rank will
28 # have no block distributed frequencies when world.size == 8
29 frq_w = np.linspace(-0.6, 1.2, 49)
30 eta = 0.2
32 rshelmax = 0
33 hxc_scaling = FMGoldstoneScaling()
34 ecut = 250
35 reduced_ecut = 100 # ecut for eigenmode analysis
36 pos_eigs = 2 # majority modes
37 neg_eigs = 0 # minority modes
38 nblocks = 'max'
40 # ---------- Script ---------- #
42 # Initialize objects to calculat Chi
43 context = ResponseContext()
44 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1))
45 nbands = response_band_cutoff['co_pw']
46 gs = ResponseGroundStateAdapter(calc)
47 chiks_calc = ChiKSCalculator(gs, context,
48 nbands=nbands,
49 ecut=ecut,
50 gammacentered=True,
51 nblocks=nblocks)
52 fxc_calculator = AdiabaticFXCCalculator.from_rshe_parameters(
53 gs, context, rshelmax=rshelmax)
54 chi_factory = ChiFactory(chiks_calc, fxc_calculator)
56 for q, q_c in enumerate(q_qc):
57 complex_frequencies = frq_w + 1.j * eta
58 chiks, chi = chi_factory('+-', q_c, complex_frequencies,
59 fxc=fxc, hxc_scaling=hxc_scaling)
61 # Check that the dissipative part of the susceptibility can be
62 # reconstructed from the EigendecomposedSpectrum
63 chid = chi.copy_dissipative_part()
64 spectrum = EigendecomposedSpectrum.from_chi(chi)
65 assert np.allclose(-np.pi * spectrum.A_wGG, chid.array)
67 # Perform a spectral decomposition of the magnetic excitations in the
68 # Kohn-Sham and many-body susceptibilities, writing the full spectral
69 # weight of both the majority and minority spectra
70 Aksmaj, Aksmin = spectral_decomposition(chiks)
71 Aksmaj.write_full_spectral_weight(f'Aksmaj_q{q}.dat')
72 Aksmin.write_full_spectral_weight(f'Aksmin_q{q}.dat')
73 Amaj, Amin = spectral_decomposition(chi,
74 pos_eigs=pos_eigs,
75 neg_eigs=neg_eigs)
76 Amaj.write_full_spectral_weight(f'Amaj_q{q}.dat')
77 Amin.write_full_spectral_weight(f'Amin_q{q}.dat')
78 # Write also the majority magnon modes of the many-body system
79 Amaj.write_eigenmode_lineshapes(f'Amaj_modes_q{q}.dat', nmodes=2)
81 # Perform an analogous spectral decomposition of the many-body
82 # susceptibility in a reduced plane-wave basis and write the decomposed
83 # spectrum to a file
84 chi = chi.copy_with_reduced_ecut(reduced_ecut)
85 chi.write_diagonal(f'chiwG_q{q}.npz')
86 Amaj, Amin = spectral_decomposition(chi,
87 pos_eigs=pos_eigs,
88 neg_eigs=neg_eigs)
89 Amaj.write(f'Amaj_q{q}.npz')
90 Amin.write(f'Amin_q{q}.npz')
91 assert f'{fxc},+-' in chi_factory.fxc_kernel_cache
93 context.write_timer()
94 context.comm.barrier()
96 # Read data
97 wksmaj0_w, Aksmaj0_w = read_full_spectral_weight('Aksmaj_q0.dat')
98 wksmin0_w, Aksmin0_w = read_full_spectral_weight('Aksmin_q0.dat')
99 wksmaj1_w, Aksmaj1_w = read_full_spectral_weight('Aksmaj_q1.dat')
100 wksmin1_w, Aksmin1_w = read_full_spectral_weight('Aksmin_q1.dat')
101 wmaj0_w, Amaj0_w = read_full_spectral_weight('Amaj_q0.dat')
102 wmin0_w, Amin0_w = read_full_spectral_weight('Amin_q0.dat')
103 wmaj1_w, Amaj1_w = read_full_spectral_weight('Amaj_q1.dat')
104 wmin1_w, Amin1_w = read_full_spectral_weight('Amin_q1.dat')
105 wa0_w, a0_wm = read_eigenmode_lineshapes('Amaj_modes_q0.dat')
106 wa1_w, a1_wm = read_eigenmode_lineshapes('Amaj_modes_q1.dat')
107 w0_w, _, chi0_wG = read_susceptibility_array('chiwG_q0.npz')
108 w1_w, _, chi1_wG = read_susceptibility_array('chiwG_q1.npz')
109 Amaj0 = EigendecomposedSpectrum.from_file('Amaj_q0.npz')
110 Amaj1 = EigendecomposedSpectrum.from_file('Amaj_q1.npz')
112 # Find acoustic magnon mode
113 wpeak00, Ipeak00 = findpeak(w0_w, -chi0_wG[:, 0].imag)
114 wpeak01, Ipeak01 = findpeak(w1_w, -chi1_wG[:, 0].imag)
115 # Find optical magnon mode
116 wpeak10, Ipeak10 = findpeak(w0_w, -chi0_wG[:, 1].imag)
117 wpeak11, Ipeak11 = findpeak(w1_w, -chi1_wG[:, 1].imag)
119 # Extract the eigenmodes from the eigendecomposed spectrum in the reduced
120 # plane-wave basis, restricting the frequencies to be non-negative
121 wmask = frq_w >= 0.
122 wm0 = Amaj0.get_eigenmode_frequency(nmodes=pos_eigs, wmask=wmask)
123 ar0_wm = Amaj0.get_eigenmodes_at_frequency(wm0, nmodes=pos_eigs)
124 wm1 = Amaj1.get_eigenmode_frequency(nmodes=pos_eigs, wmask=wmask)
125 ar1_wm = Amaj1.get_eigenmodes_at_frequency(wm1, nmodes=pos_eigs)
127 # Find peaks in eigenmodes (in full and reduced basis)
128 mpeak00, Speak00 = findpeak(wa0_w, a0_wm[:, 0])
129 mpeak01, Speak01 = findpeak(wa1_w, a1_wm[:, 0])
130 mpeak10, Speak10 = findpeak(wa0_w, a0_wm[:, 1])
131 mpeak11, Speak11 = findpeak(wa1_w, a1_wm[:, 1])
132 mrpeak00, Srpeak00 = findpeak(Amaj0.omega_w, ar0_wm[:, 0])
133 mrpeak01, Srpeak01 = findpeak(Amaj1.omega_w, ar1_wm[:, 0])
134 mrpeak10, Srpeak10 = findpeak(Amaj0.omega_w, ar0_wm[:, 1])
135 mrpeak11, Srpeak11 = findpeak(Amaj1.omega_w, ar1_wm[:, 1])
137 # Calculate the majority spectral enhancement at the acoustic magnon maxima
138 w0 = np.argmin(np.abs(Amaj0.omega_w - wpeak00))
139 w1 = np.argmin(np.abs(Amaj1.omega_w - wpeak01))
140 enh0 = Amaj0_w[w0] / Aksmaj0_w[w0]
141 enh1 = Amaj1_w[w1] / Aksmaj1_w[w1]
143 # Calculate the minority spectral enhancement at 600 meV (corresponding to
144 # -600 meV on original frequency grid)
145 min_enh0 = Amin0_w[0] / Aksmin0_w[0]
146 min_enh1 = Amin1_w[0] / Aksmin1_w[0]
148 if context.comm.rank == 0:
149 # import matplotlib.pyplot as plt
150 # # Plot the magnon lineshapes
151 # # q=0
152 # plt.subplot(2, 3, 1)
153 # plt.plot(w0_w, -chi0_wG[:, 0].imag)
154 # plt.axvline(wpeak00, c='0.5', linewidth=0.8)
155 # plt.plot(w0_w, -chi0_wG[:, 1].imag)
156 # plt.axvline(wpeak10, c='0.5', linewidth=0.8)
157 # plt.plot(Amaj0.omega_w, Amaj0.s_we[:, 0])
158 # plt.plot(Amaj0.omega_w, Amaj0.s_we[:, 1])
159 # plt.plot(wa0_w, a0_wm[:, 0])
160 # plt.plot(wa0_w, a0_wm[:, 1])
161 # plt.plot(Amaj0.omega_w, ar0_wm[:, 0])
162 # plt.plot(Amaj0.omega_w, ar0_wm[:, 1])
163 # # q=1
164 # plt.subplot(2, 3, 4)
165 # plt.plot(w1_w, -chi1_wG[:, 0].imag)
166 # plt.axvline(wpeak01, c='0.5', linewidth=0.8)
167 # plt.plot(w1_w, -chi1_wG[:, 1].imag)
168 # plt.axvline(wpeak11, c='0.5', linewidth=0.8)
169 # plt.plot(Amaj1.omega_w, Amaj1.s_we[:, 0])
170 # plt.plot(Amaj1.omega_w, Amaj1.s_we[:, 1])
171 # plt.plot(wa1_w, a1_wm[:, 0])
172 # plt.plot(wa1_w, a1_wm[:, 1])
173 # plt.plot(Amaj1.omega_w, ar1_wm[:, 0])
174 # plt.plot(Amaj1.omega_w, ar1_wm[:, 1])
175 # # Plot full spectral weight of majority excitations
176 # # q=0
177 # plt.subplot(2, 3, 2)
178 # plt.plot(wmaj0_w, Amaj0_w)
179 # plt.plot(wksmaj0_w, Aksmaj0_w)
180 # plt.plot(Amaj0.omega_w, Amaj0.A_w)
181 # plt.axvline(wpeak00, c='0.5', linewidth=0.8)
182 # # q=1
183 # plt.subplot(2, 3, 5)
184 # plt.plot(wmaj1_w, Amaj1_w)
185 # plt.plot(wksmaj1_w, Aksmaj1_w)
186 # plt.plot(Amaj1.omega_w, Amaj1.A_w)
187 # plt.axvline(wpeak01, c='0.5', linewidth=0.8)
188 # # Plot full spectral weight of minority excitations
189 # # q=0
190 # plt.subplot(2, 3, 3)
191 # plt.plot(wmin0_w, Amin0_w)
192 # plt.plot(wksmin0_w, Aksmin0_w)
193 # # q=1
194 # plt.subplot(2, 3, 6)
195 # plt.plot(wmin1_w, Amin1_w)
196 # plt.plot(wksmin1_w, Aksmin1_w)
197 # plt.show()
199 # Print values
200 print(Amaj0.omega_w[wm0], Amaj1.omega_w[wm1])
201 print(hxc_scaling.lambd)
202 print(wpeak00, wpeak01, wpeak10, wpeak11)
203 print(Ipeak00, Ipeak01, Ipeak10, Ipeak11)
204 print(mpeak00, mpeak01, mpeak10, mpeak11)
205 print(Speak00, Speak01, Speak10, Speak11)
206 print(mrpeak00, mrpeak01, mrpeak10, mrpeak11)
207 print(Srpeak00, Srpeak01, Srpeak10, Srpeak11)
208 print(enh0, enh1, min_enh0, min_enh1)
210 # Test that the mode extraction frequency is indeed non-negative
211 assert Amaj0.omega_w[wm0] >= 0.
212 assert Amaj1.omega_w[wm1] >= 0.
214 # Test kernel scaling
215 assert hxc_scaling.lambd == pytest.approx(0.9675, abs=0.005)
217 # Test magnon frequencies
218 assert wpeak00 == pytest.approx(-0.0281, abs=0.005)
219 assert wpeak01 == pytest.approx(0.218, abs=0.01)
220 assert wpeak10 == pytest.approx(0.508, abs=0.01)
221 assert wpeak11 == pytest.approx(0.637, abs=0.01)
223 # Test magnon amplitudes
224 assert Ipeak00 == pytest.approx(2.897, abs=0.01)
225 assert Ipeak01 == pytest.approx(2.245, abs=0.01)
226 assert Ipeak10 == pytest.approx(1.090, abs=0.01)
227 assert Ipeak11 == pytest.approx(1.023, abs=0.01)
229 # Test magnon frequency consistency
230 assert mpeak00 == pytest.approx(wpeak00, abs=0.005)
231 assert mpeak01 == pytest.approx(wpeak01, abs=0.01)
232 assert mpeak10 == pytest.approx(wpeak10, abs=0.01)
233 assert mpeak11 == pytest.approx(wpeak11, abs=0.01)
234 assert mrpeak00 == pytest.approx(wpeak00, abs=0.005)
235 assert mrpeak01 == pytest.approx(wpeak01, abs=0.01)
236 assert mrpeak10 == pytest.approx(wpeak10, abs=0.01)
237 assert mrpeak11 == pytest.approx(wpeak11, abs=0.01)
239 # Test magnon mode eigenvalues at extrema
240 assert Speak00 == pytest.approx(8.409, abs=0.02)
241 assert Speak01 == pytest.approx(6.734, abs=0.02)
242 assert Speak10 == pytest.approx(3.800, abs=0.02)
243 assert Speak11 == pytest.approx(3.683, abs=0.02)
244 assert Srpeak00 == pytest.approx(6.402, abs=0.02)
245 assert Srpeak01 == pytest.approx(5.087, abs=0.02)
246 assert Srpeak10 == pytest.approx(2.837, abs=0.02)
247 assert Srpeak11 == pytest.approx(2.692, abs=0.02)
249 # Test enhancement factors
250 assert enh0 == pytest.approx(36.77, abs=0.1)
251 assert enh1 == pytest.approx(24.10, abs=0.1)
252 assert min_enh0 == pytest.approx(1.141, abs=0.01)
253 assert min_enh1 == pytest.approx(1.162, abs=0.01)