Coverage for gpaw/test/corehole/test_sh2.py: 100%

1import pytest 

2from ase.build import molecule 

3 

4import gpaw.mpi as mpi 

5from gpaw import GPAW 

6from gpaw.test import gen 

7from gpaw.xas import XAS 

8 

9 

10def folding_is_normalized(xas: XAS, dks, rel: float = 1e-5) -> bool: 

11 _, ys_cn = xas.get_oscillator_strength(dks=dks) 

12 

13 ys_summed_c = ys_cn.sum(axis=1) 

14 xf, yf_cn = xas.get_spectra(fwhm=0.5, dks=dks) 

15 dxf = xf[1:] - xf[:-1] 

16 assert dxf == pytest.approx(dxf[0]) 

17 yf_summed_c = yf_cn.sum(axis=1) * dxf[0] 

18 

19 return yf_summed_c == pytest.approx(ys_summed_c, rel=rel) 

20 

21 

22@pytest.fixture 

23def s1s1ch_name(): 

24 setupname = 'S1s1ch' 

25 gen('S', name=setupname, corehole=(1, 0, 1), gpernode=30, write_xml=True) 

26 return setupname 

27 

28 

29@pytest.fixture 

30def s2p1ch_name(): 

31 setupname = 'S2p1ch' 

32 gen('S', name=setupname, corehole=(2, 1, 1), gpernode=30, write_xml=True) 

33 return setupname 

34 

35 

36def sh2(setupname): 

37 atoms = molecule('SH2') 

38 atoms.center(3) 

39 

40 nbands = 6 

41 atoms.calc = GPAW(mode='fd', h=0.3, nbands=nbands, 

42 setups={'S': setupname}, txt=None) 

43 atoms.get_potential_energy() 

44 

45 return atoms 

46 

47 

48@pytest.fixture 

49def sh2_s1s1ch(s1s1ch_name): 

50 return sh2(s1s1ch_name) 

51 

52 

53@pytest.fixture 

54def sh2_s2p1ch(s2p1ch_name): 

55 return sh2(s2p1ch_name) 

56 

57 

58def test_sulphur_2p_spin_io(in_tmp_dir, add_cwd_to_setup_paths, s2p1ch_name): 

59 """Make sure this calculation does not fail 

60 because of get_spin_contamination""" 

61 atoms = molecule('SH2') 

62 atoms.center(3) 

63 

64 atoms.set_initial_magnetic_moments([1, 0, 0]) 

65 atoms.calc = GPAW(mode='fd', h=0.3, spinpol=True, 

66 setups={'S': s2p1ch_name}, txt=None, 

67 convergence={ 

68 'energy': 0.1, 'density': 0.1, 'eigenstates': 0.1}) 

69 atoms.get_potential_energy() 

70 

71 

72def test_sulphur_1s_xas_tp(in_tmp_dir, add_cwd_to_setup_paths, sh2_s1s1ch): 

73 nbands = 6 

74 nocc = 4 # for SH2 

75 

76 dks = 20 

77 xas = XAS(sh2_s1s1ch.calc) 

78 x, y_cn = xas.get_oscillator_strength(dks=dks) 

79 assert y_cn.shape == (3, nbands - nocc) 

80 assert x[0] == dks 

81 assert xas.nocc == nocc 

82 

83 assert folding_is_normalized(xas, dks) 

84 

85 

86def test_sulphur_1s_xas_XCH(in_tmp_dir, add_cwd_to_setup_paths, sh2_s1s1ch): 

87 atoms = sh2_s1s1ch 

88 

89 nbands = 6 

90 nocc = 4 # for SH2 

91 dks = 20 

92 

93 calc = sh2_s1s1ch.calc.new(charge=-1) 

94 atoms.calc = calc 

95 

96 atoms[0].magmom = 1 

97 atoms.get_potential_energy() 

98 

99 xas = XAS(atoms.calc, relative_index_lumo=-1) 

100 x, y_cn = xas.get_oscillator_strength(dks=dks) 

101 assert xas.nocc == nocc 

102 assert y_cn.shape == (3, nbands - nocc) 

103 assert x[0] == dks 

104 assert folding_is_normalized(xas, dks) 

105 

106 

107def test_sulphur_2p_xas(in_tmp_dir, add_cwd_to_setup_paths, sh2_s2p1ch): 

108 dks = 20 

109 

110 xas = XAS(sh2_s2p1ch.calc) 

111 assert folding_is_normalized(xas, dks) 

112 

113 

114def test_proj(in_tmp_dir, add_cwd_to_setup_paths, sh2_s1s1ch): 

115 atoms = sh2_s1s1ch 

116 

117 dks = 20 

118 xas0 = XAS(atoms.calc) 

119 mefname = 'me.dat.npz' 

120 proj = [[1, 0, 0]] 

121 xas0.write(mefname) 

122 x0, y0_cn = xas0.get_oscillator_strength( 

123 dks=dks, proj=proj, proj_xyz=False) 

124 x1, y1_cn = xas0.get_oscillator_strength( 

125 dks=dks, proj_xyz=True) 

126 

127 xas1 = XAS().restart(mefname) 

128 x0_1, y0_1_cn = xas1.get_oscillator_strength( 

129 dks=dks, proj=proj, proj_xyz=False) 

130 x1_1, y1_1_cn = xas1.get_oscillator_strength( 

131 dks=dks, proj_xyz=True) 

132 

133 assert y1_cn.shape[0] == y0_cn.shape[0] + 2 

134 assert y1_cn.shape[1] == y0_cn.shape[1] 

135 assert x1 == pytest.approx(x0) 

136 assert y1_cn[0] == pytest.approx(y0_cn[0]) 

137 

138 assert x0 == pytest.approx(x0_1) 

139 assert x1 == pytest.approx(x1_1) 

140 assert y0_cn == pytest.approx(y0_1_cn) 

141 assert y1_cn == pytest.approx(y1_1_cn) 

142 

143 

144def test_parallel(in_tmp_dir, add_cwd_to_setup_paths, s2p1ch_name): 

145 atoms = molecule('SH2') 

146 atoms.center(3) 

147 

148 # serial calculation 

149 fserial = f'serial_xas_rank{mpi.world.rank}.npz' 

150 comm = mpi.world.new_communicator([mpi.world.rank]) 

151 print('serial, rank, size:', mpi.world.rank, comm.size) 

152 atoms.calc = GPAW(mode='fd', h=0.3, setups={'S': s2p1ch_name}, 

153 txt=None, communicator=comm) 

154 

155 print('serial, atoms.calc.world.size:', atoms.calc.world.size) 

156 atoms.get_potential_energy() 

157 

158 import time 

159 t0 = time.time() 

160 xas_s = XAS(atoms.calc) 

161 xas_s.write(fserial) 

162 t1 = time.time() 

163 print(t1 - t0) 

164 

165 # parallel calculation 

166 fparallel = 'parallel_xas.npz' 

167 atoms.calc = GPAW(mode='fd', h=0.3, setups={'S': s2p1ch_name}, txt=None) 

168 print('parallel, atoms.calc.world.size:', atoms.calc.world.size) 

169 atoms.get_potential_energy() 

170 

171 t0 = time.time() 

172 xas_p = XAS(atoms.calc) 

173 xas_p.write(fparallel) 

174 t1 = time.time() 

175 print(t1 - t0) 

176 

177 dks = 20 

178 xs, ys = xas_s.get_oscillator_strength(dks=dks) 

179 xp, yp = xas_p.get_oscillator_strength(dks=dks) 

180 

181 assert xs == pytest.approx(xp) 

182 assert ys == pytest.approx(yp) 

183 

184 assert xs == pytest.approx(xp) 

185 assert ys == pytest.approx(yp) 

186 

187 

188def test_io(in_tmp_dir, add_cwd_to_setup_paths, sh2_s1s1ch): 

189 """Test that a direct calculation gives the same results as a calculation 

190 restarted from matrix element file.""" 

191 dks = 20 

192 medata = 'xasme.dat' 

193 

194 xas1 = XAS(sh2_s1s1ch.calc) 

195 xas1.write(medata) 

196 

197 # define the XAS object by reading 

198 xas2 = XAS().restart(medata) 

199 

200 x1, y1 = xas1.get_oscillator_strength(dks=dks) 

201 x2, y2 = xas2.get_oscillator_strength(dks=dks) 

202 assert x1 == pytest.approx(x2) 

203 assert y1 == pytest.approx(y2)