Coverage for gpaw/test/response/test_iron_sf_pawALDA.py: 70%

1import pytest 

2 

3import numpy as np 

4 

5from gpaw import GPAW 

6from gpaw.mpi import world 

7from gpaw.response import ResponseGroundStateAdapter, ResponseContext 

8from gpaw.response.chiks import ChiKSCalculator, SelfEnhancementCalculator 

9from gpaw.response.frequencies import ComplexFrequencyDescriptor 

10from gpaw.response.dyson import DysonSolver 

11from gpaw.response.susceptibility import (spectral_decomposition, 

12 read_eigenmode_lineshapes) 

13from gpaw.response.pair_functions import read_pair_function 

14 

15from gpaw.test import findpeak 

16from gpaw.test.gpwfile import response_band_cutoff 

17 

18 

19@pytest.mark.kspair 

20@pytest.mark.response 

21def test_response_iron_sf_pawALDA(in_tmp_dir, gpw_files, scalapack): 

22 # ---------- Inputs ---------- # 

23 

24 # Magnetic response calculation 

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

26 frq_w = np.linspace(-1.0, 3.0, 161) 

27 ecut = 100 

28 eta = 0.5 

29 rshewmin = 1e-8 

30 

31 if world.size > 1: 

32 nblocks = 2 

33 else: 

34 nblocks = 1 

35 

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

37 

38 context = ResponseContext(txt='iron_susceptibility.txt') 

39 calc = GPAW(gpw_files['fe_pw'], parallel=dict(domain=1)) 

40 nbands = response_band_cutoff['fe_pw'] 

41 gs = ResponseGroundStateAdapter(calc) 

42 

43 calc_args = (gs,) 

44 calc_kwargs = dict(context=context, 

45 nbands=nbands, 

46 ecut=ecut, 

47 gammacentered=True, 

48 bandsummation='double', 

49 nblocks=nblocks) 

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

51 xi_calc = SelfEnhancementCalculator(*calc_args, 

52 rshewmin=rshewmin, 

53 **calc_kwargs) 

54 dyson_solver = DysonSolver(context) 

55 

56 for q, q_c in enumerate(q_qc): 

57 # Calculate χ_KS^+- and Ξ^++ 

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

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

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

61 

62 # Distribute frequencies and invert dyson equation 

63 chiks = chiks.copy_with_global_frequency_distribution() 

64 xi = xi.copy_with_global_frequency_distribution() 

65 chi = dyson_solver(chiks, xi) 

66 

67 # plot_inverse_enhancement(xi) 

68 

69 # Write macroscopic component and acoustic magnon mode 

70 chi.write_macroscopic_component(f'iron_chiM_q{q}.csv') 

71 Amaj, _ = spectral_decomposition(chi) 

72 Amaj.write_eigenmode_lineshapes(f'iron_Amaj_q{q}.csv') 

73 

74 context.write_timer() 

75 

76 # plot_magnons() 

77 

78 refs_q = [ 

79 # (wpeak, Ipeak, Apeak) 

80 (0.001, 0.476, 3.104), 

81 (0.165, 0.428, 2.816)] 

82 

83 world.barrier() # wait for csv-file written above ... 

84 for q, refs in enumerate(refs_q): 

85 w_w, chiM_w, a_w = extract_data(q) 

86 wpeak1, Ipeak = findpeak(w_w, -chiM_w.imag / np.pi) 

87 wpeak2, Apeak = findpeak(w_w, a_w) 

88 

89 assert wpeak1 == pytest.approx(wpeak2, abs=0.002) # eV 

90 assert wpeak1 == pytest.approx(refs[0], abs=0.01) # eV 

91 assert Ipeak == pytest.approx(refs[1], abs=0.01) # a.u. 

92 assert Apeak == pytest.approx(refs[2], abs=0.05) # a.u. 

93 

94 

95def extract_data(q): 

96 w1_w, chiM_w = read_pair_function(f'iron_chiM_q{q}.csv') 

97 w2_w, a_wm = read_eigenmode_lineshapes(f'iron_Amaj_q{q}.csv') 

98 assert np.allclose(w1_w, w2_w) 

99 return w1_w, chiM_w, a_wm[:, 0] 

100 

101 

102def plot_inverse_enhancement(xi): 

103 import matplotlib.pyplot as plt 

104 from ase.units import Ha 

105 invenh_mywM, _ = np.linalg.eig(xi.array) 

106 invenh_wM = xi.blocks1d.all_gather(invenh_mywM) 

107 omega_w = xi.zd.omega_w * Ha 

108 

109 for M in range(invenh_wM.shape[1]): 

110 plt.subplot(1, 2, 1) 

111 plt.scatter(omega_w, invenh_wM[:, M].real) 

112 plt.subplot(1, 2, 2) 

113 plt.scatter(omega_w, invenh_wM[:, M].imag) 

114 plt.subplot(1, 2, 1) 

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

116 if world.rank == 0: 

117 plt.show() 

118 

119 

120def plot_magnons(): 

121 import matplotlib.pyplot as plt 

122 for q in range(2): 

123 w_w, chiM_w, a_w = extract_data(q) 

124 plt.subplot(1, 2, 1) 

125 plt.plot(w_w, - chiM_w.imag / np.pi, label=f'{q}') 

126 plt.subplot(1, 2, 2) 

127 plt.plot(w_w, a_w, label=f'{q}') 

128 plt.legend(title='q') 

129 if world.rank == 0: 

130 plt.show()