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

1import pytest 

2 

3import numpy as np 

4 

5from gpaw import GPAW 

6from gpaw.response import ResponseGroundStateAdapter, ResponseContext 

7from gpaw.response.frequencies import ComplexFrequencyDescriptor 

8from gpaw.response.chiks import ChiKSCalculator, SelfEnhancementCalculator 

9from gpaw.response.dyson import DysonSolver 

10from gpaw.response.goldstone import ( 

11 NewFMGoldstoneScaling, 

12 RefinedFMGoldstoneScaling, 

13) 

14from gpaw.response.susceptibility import (spectral_decomposition, 

15 read_eigenmode_lineshapes) 

16 

17from gpaw.test import findpeak 

18from gpaw.test.gpwfile import response_band_cutoff 

19 

20 

21@pytest.mark.kspair 

22@pytest.mark.response 

23def test_response_cobalt_sf_gsspawALDA(in_tmp_dir, gpw_files): 

24 # ---------- Inputs ---------- # 

25 

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

27 frq_w = np.linspace(-0.5, 2.0, 101) 

28 eta = 0.2 

29 

30 rshelmax = 0 

31 ecut = 150 

32 pos_eigs = 5 

33 nmodes = 2 # majority modes 

34 nblocks = 'max' 

35 

36 # ---------- Script ---------- # 

37 

38 # Read ground state data 

39 context = ResponseContext(txt='cobalt_susceptibility.txt') 

40 calc = GPAW(gpw_files['co_pw'], parallel=dict(domain=1)) 

41 nbands = response_band_cutoff['co_pw'] 

42 gs = ResponseGroundStateAdapter(calc) 

43 

44 # Set up response calculators 

45 calc_args = (gs,) 

46 calc_kwargs = dict(context=context, 

47 nbands=nbands, 

48 ecut=ecut, 

49 gammacentered=True, 

50 bandsummation='pairwise', 

51 nblocks=nblocks) 

52 chiks_calc = ChiKSCalculator(*calc_args, **calc_kwargs) 

53 xi_calc = SelfEnhancementCalculator(*calc_args, 

54 rshelmax=rshelmax, 

55 **calc_kwargs) 

56 base_scaling = NewFMGoldstoneScaling.from_xi_calculator(xi_calc) 

57 refi_scaling = RefinedFMGoldstoneScaling.from_xi_calculator(xi_calc) 

58 dyson_solver = DysonSolver(context) 

59 

60 for q, q_c in enumerate(q_qc): 

61 # Calculate χ_KS^+-(q,z) and Ξ^++(q,z) 

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

63 chiks = chiks_calc.calculate('+-', q_c, zd) 

64 xi = xi_calc.calculate('+-', q_c, zd) 

65 

66 # Distribute frequencies and invert dyson equation 

67 chiks = chiks.copy_with_global_frequency_distribution() 

68 xi = xi.copy_with_global_frequency_distribution() 

69 chi = dyson_solver(chiks, xi) 

70 # Test ability to apply Goldstone scaling when inverting the dyson eq. 

71 scaled_chi = dyson_solver(chiks, xi, hxc_scaling=base_scaling) 

72 refined_chi = dyson_solver(chiks, xi, hxc_scaling=refi_scaling) 

73 # Simulate a "restart" of the GoldstoneScaling 

74 base_scaling = NewFMGoldstoneScaling(lambd=base_scaling.lambd) 

75 refi_scaling = RefinedFMGoldstoneScaling(lambd=refi_scaling.lambd) 

76 

77 # Calculate majority spectral function 

78 Amaj, _ = spectral_decomposition(chi, pos_eigs=pos_eigs) 

79 Amaj.write_eigenmode_lineshapes( 

80 f'cobalt_Amaj_q{q}.csv', nmodes=nmodes) 

81 sAmaj, _ = spectral_decomposition(scaled_chi, pos_eigs=pos_eigs) 

82 sAmaj.write_eigenmode_lineshapes( 

83 f'cobalt_sAmaj_q{q}.csv', nmodes=nmodes) 

84 rAmaj, _ = spectral_decomposition(refined_chi, pos_eigs=pos_eigs) 

85 rAmaj.write_eigenmode_lineshapes( 

86 f'cobalt_rAmaj_q{q}.csv', nmodes=nmodes) 

87 

88 # Store Re ξ^++(q=0,ω), to test the self-enhancement after scaling 

89 if q == 0: 

90 _, xi_mW = get_mode_projections( 

91 chiks, xi, sAmaj, lambd=base_scaling.lambd, nmodes=nmodes) 

92 xi0_w = xi_mW[0].real 

93 _, xi_mW = get_mode_projections( 

94 chiks, xi, sAmaj, lambd=refi_scaling.lambd, nmodes=nmodes) 

95 rxi0_w = xi_mW[0].real 

96 

97 # plot_enhancement(chiks, xi, Amaj, sAmaj, 

98 # lambd=base_scaling.lambd, nmodes=nmodes) 

99 # plot_enhancement(chiks, xi, Amaj, rAmaj, 

100 # lambd=refi_scaling.lambd, nmodes=nmodes) 

101 

102 context.write_timer() 

103 

104 # Compare scaling coefficient to reference 

105 assert base_scaling.lambd == pytest.approx(1.0541, abs=0.001) 

106 assert refi_scaling.lambd == pytest.approx(1.0526, abs=0.001) 

107 # Test that Re ξ^++(q=0,ω) ≾ 1 at ω=0 

108 w0 = np.argmin(np.abs(frq_w)) 

109 assert xi0_w[w0] == pytest.approx(0.987, abs=0.01) 

110 assert rxi0_w[w0] == pytest.approx(0.986, abs=0.01) 

111 

112 # Compare magnon peaks to reference data 

113 refs_mqa = [ 

114 # Acoustic 

115 [ 

116 # q_Γ 

117 [ 

118 # (wpeak, Apeak) 

119 (0.085, 7.895), # unscaled 

120 (-0.002, 7.980), # scaled 

121 ], 

122 # q_M / 2 

123 [ 

124 # (wpeak, Apeak) 

125 (0.320, 5.828), # unscaled 

126 (0.245, 6.215), # scaled 

127 ], 

128 ], 

129 # Optical 

130 [ 

131 # q_Γ 

132 [ 

133 # (wpeak, Apeak) 

134 (0.904, 3.493), # unscaled 

135 (0.860, 3.395), # scaled 

136 ], 

137 # q_M / 2 

138 [ 

139 # (wpeak, Apeak) 

140 (0.857, 2.988), # unscaled 

141 (0.721, 3.163), # scaled 

142 ], 

143 ], 

144 ] 

145 for a, Astr in zip([0, 1, 1], ['Amaj', 'sAmaj', 'rAmaj']): 

146 for q in range(len(q_qc)): 

147 w_w, a_wm = read_eigenmode_lineshapes(f'cobalt_{Astr}_q{q}.csv') 

148 for m in range(nmodes): 

149 wpeak, Apeak = findpeak(w_w, a_wm[:, m]) 

150 refw, refA = refs_mqa[m][q][a] 

151 print(m, q, a, wpeak, Apeak) 

152 assert wpeak == pytest.approx(refw, abs=0.01) # eV 

153 assert Apeak == pytest.approx(refA, abs=0.05) # a.u. 

154 if q == 0 and m == 0 and Astr == 'rAmaj': 

155 # Check that the gap error is completely removed when using 

156 # the refined scaling 

157 assert wpeak == pytest.approx(0.0, abs=1e-4) # eV 

158 

159 

160def get_mode_projections(chiks, xi, Amaj, *, lambd, nmodes): 

161 """Project χ_KS^+-(q,z) and Ξ^++(q,z) onto the magnon mode vectors.""" 

162 wm = Amaj.get_eigenmode_frequency(nmodes=nmodes) 

163 v_Gm = Amaj.get_eigenvectors_at_frequency(wm, nmodes=nmodes) 

164 chiks_wm = np.zeros((chiks.blocks1d.nlocal, nmodes), dtype=complex) 

165 xi_wm = np.zeros((xi.blocks1d.nlocal, nmodes), dtype=complex) 

166 for m, v_G in enumerate(v_Gm.T): 

167 chiks_wm[:, m] = np.conj(v_G) @ chiks.array @ v_G # chiks_wGG 

168 xi_wm[:, m] = np.conj(v_G) @ xi.array @ v_G # xi_wGG 

169 chiks_mW = chiks.blocks1d.all_gather(chiks_wm * lambd).T 

170 xi_mW = xi.blocks1d.all_gather(xi_wm * lambd).T 

171 return chiks_mW, xi_mW 

172 

173 

174def plot_enhancement(chiks, xi, Amaj0, sAmaj, *, lambd, nmodes): 

175 import matplotlib.pyplot as plt 

176 from gpaw.mpi import world 

177 from ase.units import Ha 

178 

179 for Amaj, _lambd in zip([Amaj0, sAmaj], [1., lambd]): 

180 a_mW = Amaj.get_eigenmode_lineshapes(nmodes=nmodes).T 

181 chiks_mW, xi_mW = get_mode_projections( 

182 chiks, xi, Amaj, lambd=_lambd, nmodes=nmodes) 

183 for m in range(nmodes): 

184 plt.subplot(1, nmodes, m + 1) 

185 plt.plot(chiks.zd.omega_w * Ha, -chiks_mW[m].imag / np.pi) 

186 plt.plot(xi.zd.omega_w * Ha, xi_mW[m].real) 

187 plt.axhline(1., c='0.5') 

188 plt.plot(Amaj.omega_w, a_mW[m]) 

189 

190 if world.rank == 0: 

191 plt.show() 

192 world.barrier()