Coverage for gpaw/test/response/test_afm_hchain_sf_gssALDA.py: 98%

1"""Simple magnons calculation in a 1D hydrogen chain.""" 

2 

3# General modules 

4import pytest 

5import numpy as np 

6 

7# Script modules 

8from ase import Atoms 

9from ase.dft.kpoints import monkhorst_pack 

10 

11from gpaw import PW, GPAW 

12from gpaw.mpi import world 

13from gpaw.test import findpeak 

14 

15from gpaw.response import ResponseGroundStateAdapter 

16from gpaw.response.frequencies import ComplexFrequencyDescriptor 

17from gpaw.response.chiks import ChiKSCalculator 

18from gpaw.response.susceptibility import ChiFactory 

19from gpaw.response.fxc_kernels import AdiabaticFXCCalculator 

20from gpaw.response.goldstone import AFMGoldstoneScaling 

21from gpaw.response.pair_functions import read_pair_function 

22 

23 

24@pytest.mark.old_gpaw_only # interpolate=3 for PW-mode not implemented! 

25@pytest.mark.kspair 

26@pytest.mark.response 

27@pytest.mark.parametrize('from_file', [False, True]) 

28def test_response_afm_hchain_gssALDA(in_tmp_dir, from_file): 

29 # ---------- Inputs ---------- # 

30 

31 # Part 1: Ground state calculation 

32 # Define atomic structure 

33 a = 2.5 

34 vfactor = 0.8 

35 mm = 1. 

36 # Ground state calculator options 

37 xc = 'LDA' 

38 kpts = 12 

39 nbands = 2 * (1 + 0) # 1s + 0 empty shell bands 

40 ebands = 2 * 1 # Include also 2s bands for numerical consistency 

41 pw = 250 

42 conv = {'bands': nbands} 

43 

44 # # Part 2: Magnetic response calculation 

45 q_qc = [[0., 0., 0.], 

46 [1. / 6., 0., 0.], 

47 [1. / 3., 0., 0.]] 

48 fxc = 'ALDA' 

49 hxc_scaling = AFMGoldstoneScaling() 

50 rshelmax = -1 

51 rshewmin = 1e-8 

52 ecut = 120 

53 frq_w = np.linspace(-0.6, 0.6, 41) 

54 eta = 0.24 

55 zd = ComplexFrequencyDescriptor.from_array(frq_w + 1.j * eta) 

56 if world.size % 4 == 0: 

57 nblocks = 4 

58 elif world.size % 2 == 0: 

59 nblocks = 2 

60 else: 

61 nblocks = 1 

62 

63 # ---------- Script ---------- # 

64 

65 # Part 1: Ground state calculation 

66 

67 Hatom = Atoms('H', 

68 cell=[a, 0, 0], 

69 # Use pbc to allow for real-space density interpolation 

70 pbc=[1, 1, 1]) 

71 Hatom.center(vacuum=vfactor * a, axis=(1, 2)) 

72 Hchain = Hatom.repeat((2, 1, 1)) 

73 Hchain.set_initial_magnetic_moments([mm, -mm]) 

74 

75 calc = GPAW(xc=xc, 

76 mode=PW(pw, 

77 # Interpolate the density in real space 

78 interpolation=3), 

79 kpts=monkhorst_pack((kpts, 1, 1)), 

80 nbands=nbands + ebands, 

81 convergence=conv, 

82 symmetry={'point_group': True}, 

83 parallel={'domain': 1}) 

84 

85 Hchain.calc = calc 

86 Hchain.get_potential_energy() 

87 

88 # Part 2: Magnetic response calculation 

89 if from_file: 

90 calc.write('gs.gpw', mode='all') 

91 gs = ResponseGroundStateAdapter.from_gpw_file('gs.gpw') 

92 else: 

93 gs = ResponseGroundStateAdapter(calc) 

94 chiks_calc = ChiKSCalculator(gs, 

95 nbands=nbands, 

96 ecut=ecut, 

97 gammacentered=True, 

98 nblocks=nblocks) 

99 fxc_calculator = AdiabaticFXCCalculator.from_rshe_parameters( 

100 gs, chiks_calc.context, 

101 rshelmax=rshelmax, 

102 rshewmin=rshewmin) 

103 chi_factory = ChiFactory(chiks_calc, fxc_calculator) 

104 

105 for q, q_c in enumerate(q_qc): 

106 filename = 'h-chain_macro_tms_q%d.csv' % q 

107 txt = 'h-chain_macro_tms_q%d.txt' % q 

108 _, chi = chi_factory('+-', q_c, zd, 

109 fxc=fxc, 

110 hxc_scaling=hxc_scaling, 

111 txt=txt) 

112 chi.write_macroscopic_component(filename) 

113 

114 chi_factory.context.write_timer() 

115 world.barrier() 

116 

117 # Part 3: Identify magnon peak in finite q scattering function 

118 w0_w, chi0_w = read_pair_function('h-chain_macro_tms_q0.csv') 

119 w1_w, chi1_w = read_pair_function('h-chain_macro_tms_q1.csv') 

120 w2_w, chi2_w = read_pair_function('h-chain_macro_tms_q2.csv') 

121 

122 wpeak1, Ipeak1 = findpeak(w1_w, -chi1_w.imag / np.pi) 

123 wpeak2, Ipeak2 = findpeak(w2_w, -chi2_w.imag / np.pi) 

124 mw1 = calculate_afm_mw(wpeak1, eta) * 1000 # to meV 

125 mw2 = calculate_afm_mw(wpeak2, eta) * 1000 # to meV 

126 

127 # Part 4: compare new results to test values 

128 test_fxcs = 1.04744 

129 test_mw1 = 285. # meV 

130 test_mw2 = 494. # meV # Remark that mw2 < 2 * mw1 (linear dispersion) 

131 test_Ipeak1 = 0.0131 # a.u. 

132 test_Ipeak2 = 0.0290 

133 

134 # Test fxc_scaling: 

135 fxcs = hxc_scaling.lambd 

136 assert abs(fxcs - test_fxcs) < 0.005 

137 

138 # Magnon peak at q=1/3 q_X: 

139 assert abs(mw1 - test_mw1) < 10. 

140 

141 # Magnon peak at q=2/3 q_X: 

142 assert abs(mw2 - test_mw2) < 10. 

143 

144 # Scattering function intensity: 

145 assert abs(Ipeak1 - test_Ipeak1) < 0.005 

146 assert abs(Ipeak2 - test_Ipeak2) < 0.005 

147 

148 # Check that spectrum at q=0 vanishes 

149 chitol = 1e-3 * np.abs(chi1_w.imag).max() 

150 assert np.abs(chi0_w.imag).max() < chitol 

151 

152 # Check that the spectrum is antisymmetric around q=0 

153 assert np.allclose(w0_w[19::-1] + w0_w[21:], 0.) 

154 assert np.allclose(chi0_w.imag[19::-1] + chi0_w.imag[21:], 0., 

155 atol=0.01 * chitol) 

156 assert np.allclose(chi1_w.imag[19::-1] + chi1_w.imag[21:], 0., 

157 atol=chitol) 

158 assert np.allclose(chi2_w.imag[19::-1] + chi2_w.imag[21:], 0., 

159 atol=chitol) 

160 

161 

162# ---------- Script functionality ---------- # 

163 

164 

165def calculate_afm_mw(peak_frequency, eta): 

166 """Assume lorentzian lineshape to compute the magnon frequency 

167 from the afm magnon lineshape peak position.""" 

168 return np.sqrt(2 * np.sqrt(peak_frequency**4 + eta**2 * peak_frequency**2) 

169 - peak_frequency**2 - eta**2)