Coverage for gpaw/test/response/test_iron_sf_gssALDA.py: 85%

1""" 

2Calculate the magnetic response in iron using ALDA. 

3 

4Fast test, where the kernel is scaled to fulfill the Goldstone theorem. 

5""" 

6 

7# Workflow modules 

8from pathlib import Path 

9import pytest 

10import numpy as np 

11 

12# Script modules 

13from gpaw.test import findpeak 

14from gpaw.mpi import world 

15 

16from gpaw.response import ResponseGroundStateAdapter 

17from gpaw.response.chiks import ChiKSCalculator 

18from gpaw.response.susceptibility import (ChiFactory, spectral_decomposition, 

19 read_eigenmode_lineshapes) 

20from gpaw.response.localft import LocalGridFTCalculator, LocalPAWFTCalculator 

21from gpaw.response.fxc_kernels import FXCKernel, AdiabaticFXCCalculator 

22from gpaw.response.dyson import HXCKernel 

23from gpaw.response.goldstone import FMGoldstoneScaling 

24from gpaw.response.pair_functions import read_pair_function 

25 

26 

27def set_up_fxc_calculators(gs, context): 

28 fxckwargs_and_identifiers = [] 

29 

30 # Set up grid calculator (without file buffer) 

31 localft_calc = LocalGridFTCalculator(gs, context) 

32 fxckwargs_grid = {'fxc_calculator': AdiabaticFXCCalculator(localft_calc), 

33 'hxc_scaling': FMGoldstoneScaling()} 

34 fxckwargs_and_identifiers.append((fxckwargs_grid, 'grid')) 

35 

36 # Set up paw calculator (with file buffer) 

37 localft_calc = LocalPAWFTCalculator(gs, context, rshelmax=0) 

38 fxckwargs_paw = {'fxc_calculator': AdiabaticFXCCalculator(localft_calc), 

39 'fxc_file': Path('paw_ALDA_fxc.npz'), 

40 'hxc_scaling': FMGoldstoneScaling()} 

41 fxckwargs_and_identifiers.append((fxckwargs_paw, 'paw')) 

42 

43 return fxckwargs_and_identifiers 

44 

45 

46def get_test_values(identifier): 

47 test_mw1 = 0. # meV 

48 test_Ipeak1 = 7.48 # a.u. 

49 test_Ipeak3 = 19.83 # a.u. 

50 if identifier == 'grid': 

51 test_fxcs = 1.059 

52 test_mw2 = 363. # meV 

53 test_Ipeak2 = 3.47 # a.u. 

54 test_Ipeak4 = 9.20 # a.u. 

55 elif identifier == 'paw': 

56 test_fxcs = 1.131 

57 test_mw2 = 352. # meV 

58 test_Ipeak2 = 3.35 # a.u. 

59 test_Ipeak4 = 9.14 # a.u. 

60 else: 

61 raise ValueError(f'Invalid identifier {identifier}') 

62 

63 return (test_fxcs, test_mw1, test_mw2, 

64 test_Ipeak1, test_Ipeak2, test_Ipeak3, test_Ipeak4) 

65 

66 

67@pytest.mark.kspair 

68@pytest.mark.response 

69def test_response_iron_sf_gssALDA(in_tmp_dir, gpw_files): 

70 # ---------- Inputs ---------- # 

71 

72 nbands = 6 

73 q_qc = [[0.0, 0.0, 0.0], [0.0, 0.0, 1. / 4.]] # Two q-points along G-N 

74 frq_qw = [np.linspace(-0.080, 0.120, 26), np.linspace(0.250, 0.450, 26)] 

75 fxc = 'ALDA' 

76 ecut = 300 

77 eta = 0.1 

78 if world.size > 1: 

79 nblocks = 2 

80 else: 

81 nblocks = 1 

82 

83 # ---------- Script ---------- # 

84 

85 gs = ResponseGroundStateAdapter.from_gpw_file(gpw_files['fe_pw']) 

86 chiks_calc = ChiKSCalculator(gs, 

87 nbands=nbands, 

88 ecut=ecut, 

89 gammacentered=True, 

90 nblocks=nblocks) 

91 fxckwargs_and_identifiers = set_up_fxc_calculators( 

92 gs, chiks_calc.context) 

93 

94 chi_factory = ChiFactory( 

95 chiks_calc, 

96 # Use the first fxc_calculator for the ChiFactory 

97 fxc_calculator=fxckwargs_and_identifiers[0][0]['fxc_calculator']) 

98 

99 for q in range(2): 

100 complex_frequencies = frq_qw[q] + 1.j * eta 

101 

102 # Calculate chi using the various fxc calculators 

103 for f, (fxckwargs, identifier) in enumerate(fxckwargs_and_identifiers): 

104 if f == 0: 

105 # The first kernel is used directly through the ChiFactory 

106 chiks, chi = chi_factory('+-', q_qc[q], complex_frequencies, 

107 fxc=fxc, 

108 hxc_scaling=fxckwargs['hxc_scaling']) 

109 else: # We wish to test storing the remaining kernels in files 

110 assert 'fxc_file' in fxckwargs 

111 fxc_file = fxckwargs['fxc_file'] 

112 if q == 0: # Calculate and save kernel 

113 assert not fxc_file.is_file() 

114 fxc_calculator = fxckwargs['fxc_calculator'] 

115 fxc_kernel = fxc_calculator(fxc, '+-', chiks.qpd) 

116 fxc_kernel.save(fxc_file) 

117 else: # Reuse kernel from previous calculation 

118 assert fxc_file.is_file() 

119 fxc_kernel = FXCKernel.from_file(fxc_file) 

120 

121 # Calculate many-body susceptibility 

122 chi = chi_factory.dyson_solver( 

123 chiks, HXCKernel( 

124 hartree_kernel=chi_factory.get_hartree_kernel( 

125 '+-', chiks.qpd), 

126 xc_kernel=fxc_kernel), 

127 hxc_scaling=fxckwargs['hxc_scaling']) 

128 

129 chi.write_macroscopic_component(identifier + '_iron_dsus' 

130 + '_%d.csv' % (q + 1)) 

131 Amaj, _ = spectral_decomposition(chi, pos_eigs=10, neg_eigs=0) 

132 Amaj.write_eigenmode_lineshapes(identifier + '_Amaj_modes' 

133 + '_%d.csv' % (q + 1), nmodes=1) 

134 

135 chi_factory.context.write_timer() 

136 

137 world.barrier() 

138 

139 # plot_comparison('grid', 'paw') 

140 

141 # Compare results to test values 

142 for fxckwargs, identifier in fxckwargs_and_identifiers: 

143 fxcs = fxckwargs['hxc_scaling'].lambd 

144 (_, _, mw1, Ipeak1, _, _, mw2, Ipeak2, 

145 _, _, mw3, Ipeak3, _, _, mw4, Ipeak4) = extract_data(identifier) 

146 

147 print(fxcs, mw1, mw2, Ipeak1, Ipeak2, mw3, Ipeak3, mw4, Ipeak4) 

148 

149 (test_fxcs, test_mw1, test_mw2, 

150 test_Ipeak1, test_Ipeak2, 

151 test_Ipeak3, test_Ipeak4) = get_test_values(identifier) 

152 

153 # fxc scaling: 

154 assert fxcs == pytest.approx(test_fxcs, abs=0.005) 

155 

156 # Magnon peak: 

157 assert mw1 == pytest.approx(test_mw1, abs=20.) 

158 assert mw2 == pytest.approx(test_mw2, abs=50.) 

159 assert mw3 == pytest.approx(mw1, abs=10.) 

160 assert mw4 == pytest.approx(mw2, abs=10.) 

161 

162 # Scattering function intensity: 

163 assert Ipeak1 == pytest.approx(test_Ipeak1, abs=0.5) 

164 assert Ipeak2 == pytest.approx(test_Ipeak2, abs=0.5) 

165 assert Ipeak3 == pytest.approx(test_Ipeak3, abs=0.5) 

166 assert Ipeak4 == pytest.approx(test_Ipeak4, abs=1.1) 

167 

168 

169def extract_data(identifier): 

170 # Read data 

171 w1_w, chi1_w = read_pair_function(identifier + '_iron_dsus_1.csv') 

172 w2_w, chi2_w = read_pair_function(identifier + '_iron_dsus_2.csv') 

173 w3_w, s3_wm = read_eigenmode_lineshapes(identifier + '_Amaj_modes_1.csv') 

174 w4_w, s4_wm = read_eigenmode_lineshapes(identifier + '_Amaj_modes_2.csv') 

175 s3_w = s3_wm[:, 0] 

176 s4_w = s4_wm[:, 0] 

177 

178 # Spectral function 

179 S1_w = -chi1_w.imag 

180 S2_w = -chi2_w.imag 

181 

182 # Identify peaks 

183 wpeak1, Ipeak1 = findpeak(w1_w, S1_w) 

184 wpeak2, Ipeak2 = findpeak(w2_w, S2_w) 

185 wpeak3, Ipeak3 = findpeak(w3_w, s3_w) 

186 wpeak4, Ipeak4 = findpeak(w4_w, s4_w) 

187 

188 # Peak positions in meV 

189 mw1 = wpeak1 * 1000 

190 mw2 = wpeak2 * 1000 

191 mw3 = wpeak3 * 1000 

192 mw4 = wpeak4 * 1000 

193 

194 return (w1_w, S1_w, mw1, Ipeak1, w2_w, S2_w, mw2, Ipeak2, 

195 w3_w, s3_w, mw3, Ipeak3, w4_w, s4_w, mw4, Ipeak4) 

196 

197 

198def plot_comparison(identifier1, identifier2): 

199 (w11_w, S11_w, _, _, w12_w, S12_w, _, _, 

200 w13_w, s13_w, _, _, w14_w, s14_w, _, _) = extract_data(identifier1) 

201 (w21_w, S21_w, _, _, w22_w, S22_w, _, _, 

202 w23_w, s23_w, _, _, w24_w, s24_w, _, _) = extract_data(identifier2) 

203 

204 import matplotlib.pyplot as plt 

205 plt.subplot(2, 2, 1) 

206 plt.plot(w11_w, S11_w) 

207 plt.plot(w21_w, S21_w) 

208 plt.subplot(2, 2, 2) 

209 plt.plot(w12_w, S12_w) 

210 plt.plot(w22_w, S22_w) 

211 plt.subplot(2, 2, 3) 

212 plt.plot(w13_w, s13_w) 

213 plt.plot(w23_w, s23_w) 

214 plt.subplot(2, 2, 4) 

215 plt.plot(w14_w, s14_w) 

216 plt.plot(w24_w, s24_w) 

217 plt.show()