Coverage for gpaw/test/response/test_nicl2_sf_gssALDA.py: 28%

1import numpy as np 

2import pytest 

3 

4from gpaw.mpi import world 

5from gpaw.test import findpeak 

6 

7from gpaw.response import ResponseGroundStateAdapter 

8from gpaw.response.frequencies import ComplexFrequencyDescriptor 

9from gpaw.response.chiks import ChiKSCalculator 

10from gpaw.response.fxc_kernels import AdiabaticFXCCalculator 

11from gpaw.response.dyson import HXCKernel 

12from gpaw.response.goldstone import FMGoldstoneScaling 

13from gpaw.response.susceptibility import ChiFactory 

14from gpaw.response.pair_functions import read_pair_function 

15 

16 

17pytestmark = pytest.mark.skipif(world.size < 4, 

18 reason='too slow for world.size < 4') 

19 

20 

21@pytest.mark.kspair 

22@pytest.mark.response 

23@pytest.mark.old_gpaw_only 

24def test_nicl2_magnetic_response(in_tmp_dir, gpw_files): 

25 # ---------- Inputs ---------- # 

26 

27 q_qc = [[0., 0., 0.], 

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

29 fxc = 'ALDA' 

30 rshelmax = 0 

31 rshewmin = None 

32 bg_density = 0.004 

33 ecut = 200 

34 frq_w = np.linspace(-0.15, 0.075, 16) 

35 eta = 0.12 

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

37 nblocks = 4 

38 

39 # ---------- Script ---------- # 

40 

41 # Magnetic response calculation 

42 gs = ResponseGroundStateAdapter.from_gpw_file(gpw_files['nicl2_pw']) 

43 chiks_calc = ChiKSCalculator(gs, 

44 ecut=ecut, 

45 gammacentered=True, 

46 nblocks=nblocks) 

47 

48 # Calculate the magnetic response with and without a background density 

49 hxc_scaling = FMGoldstoneScaling() 

50 fxc_calculator = AdiabaticFXCCalculator.from_rshe_parameters( 

51 gs, chiks_calc.context, 

52 rshelmax=rshelmax, 

53 rshewmin=rshewmin) 

54 chi_factory = ChiFactory(chiks_calc, fxc_calculator) 

55 bgd_hxc_scaling = FMGoldstoneScaling() 

56 bgd_fxc_calculator = AdiabaticFXCCalculator.from_rshe_parameters( 

57 gs, chiks_calc.context, 

58 bg_density=bg_density, 

59 rshelmax=rshelmax, 

60 rshewmin=rshewmin) 

61 

62 # Check that the pseudo-density is nonnegative 

63 nt_sr, gd = gs.get_pseudo_density(gridrefinement=2) 

64 assert np.all(nt_sr >= 0.) 

65 

66 # Plot pseudo spin-density to check exponential localization into vacuum 

67 # if world.rank == 0: 

68 # import matplotlib.pyplot as plt 

69 # N_c = gd.N_c 

70 # # Plot the spin-densities above one of the Cl atoms 

71 # plt.plot(range(N_c[2]), nt_sr[0, 2 * N_c[0] // 3, N_c[1] // 3]) 

72 # plt.plot(range(N_c[2]), nt_sr[1, 2 * N_c[0] // 3, N_c[1] // 3]) 

73 # plt.yscale('log') 

74 # plt.show() 

75 

76 filestr = 'nicl2_macro_tms' 

77 bgd_filestr = 'nicl2_macro_tms_bgd' 

78 for q, q_c in enumerate(q_qc): 

79 chiks, chi = chi_factory('+-', q_c, zd, 

80 fxc=fxc, 

81 hxc_scaling=hxc_scaling, 

82 txt=filestr + '_q%d.txt' % q) 

83 chi.write_macroscopic_component(filestr + '_q%d.csv' % q) 

84 if q == 0: 

85 # Calculate kernel with background charge 

86 bgd_hxc_kernel = HXCKernel( 

87 hartree_kernel=chi_factory.get_hartree_kernel( 

88 '+-', chiks.qpd), 

89 xc_kernel=bgd_fxc_calculator(fxc, '+-', chiks.qpd)) 

90 bgd_chi = chi_factory.dyson_solver( 

91 chiks, bgd_hxc_kernel, bgd_hxc_scaling) 

92 bgd_chi.write_macroscopic_component(bgd_filestr + '_q%d.csv' % q) 

93 

94 chiks_calc.context.write_timer() 

95 world.barrier() 

96 

97 # Compare new results to test values 

98 check_magnons(filestr, hxc_scaling, 

99 test_fxcs=0.7130, 

100 test_mw0=-10.3, # meV 

101 test_mw1=-40.9, # meV 

102 test_Ipeak0=0.2306, # a.u. 

103 test_Ipeak1=0.0956, # a.u. 

104 ) 

105 check_magnons(bgd_filestr, bgd_hxc_scaling, 

106 test_fxcs=1.0826, 

107 test_mw0=-10.2, # meV 

108 test_mw1=-48.0, # meV 

109 test_Ipeak0=0.2321, # a.u. 

110 test_Ipeak1=0.0956, # a.u. 

111 ) 

112 

113 

114def check_magnons(filestr, hxc_scaling, *, 

115 test_fxcs, test_mw0, test_mw1, test_Ipeak0, test_Ipeak1): 

116 # Identify magnon peaks and extract kernel scaling 

117 w0_w, chi0_w = read_pair_function(filestr + '_q0.csv') 

118 w1_w, chi1_w = read_pair_function(filestr + '_q1.csv') 

119 

120 wpeak0, Ipeak0 = findpeak(w0_w, -chi0_w.imag / np.pi) 

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

122 mw0 = wpeak0 * 1e3 # meV 

123 mw1 = wpeak1 * 1e3 # meV 

124 

125 assert hxc_scaling.lambd is not None 

126 fxcs = hxc_scaling.lambd 

127 

128 if world.rank == 0: 

129 # import matplotlib.pyplot as plt 

130 # plt.plot(w0_w, -chi0_w.imag / np.pi) 

131 # plt.plot(w1_w, -chi1_w.imag / np.pi) 

132 # plt.show() 

133 

134 print(fxcs, mw0, mw1, Ipeak0, Ipeak1) 

135 

136 # Test fxc scaling 

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

138 

139 # Test magnon peaks 

140 assert mw0 == pytest.approx(test_mw0, abs=5.0) 

141 assert mw1 == pytest.approx(test_mw1, abs=10.0) 

142 

143 # Test peak intensities 

144 assert Ipeak0 == pytest.approx(test_Ipeak0, abs=0.01) 

145 assert Ipeak1 == pytest.approx(test_Ipeak1, abs=0.01)