Coverage for gpaw/test/ext_potential/stark_shift.py: 11%

1from math import sqrt, pi 

2 

3import pytest 

4import numpy as np 

5from ase import Atoms 

6from ase.units import Bohr, Hartree 

7 

8from gpaw.mpi import rank, size 

9from gpaw import GPAW 

10from gpaw.external import ConstantElectricField 

11from gpaw.utilities import packed_index 

12from gpaw.pair_density import PairDensity 

13 

14# Three ways to compute the polarizability of hydrogen: 

15# 1. Perturbation theory 

16# 2. Constant electric field 

17# 3. Middle of an electric dipole -- not tested here 

18 

19# Note: The analytical value for the polarizability 

20# is 4.5 a0**3 (e.g. PRA 33, 3671), while the experimental 

21# value is 4.6 a0**3 (e.g. PR 133, A629). 

22 

23 

24@pytest.mark.skip(reason='too-slow') 

25def test_stark_shift(): 

26 to_au = Hartree / Bohr**2 

27 

28 def dipole_op(c, state1, state2, k=0, s=0): 

29 # Taken from KSSingle, maybe make this accessible in 

30 # KSSingle? 

31 wfs = c.wfs 

32 gd = wfs.gd 

33 

34 kpt = None 

35 for i in wfs.kpt_u: 

36 if i.k == k and i.s == s: 

37 kpt = i 

38 

39 pd = PairDensity(c) 

40 pd.initialize(kpt, state1, state2) 

41 

42 # coarse grid contribution 

43 # <i|r|j> is the negative of the dipole moment (because of negative 

44 # e- charge) 

45 me = -gd.calculate_dipole_moment(pd.get()) 

46 

47 # augmentation contributions 

48 ma = np.zeros(me.shape) 

49 pos_av = c.get_atoms().get_positions() / Bohr 

50 for a, P_ni in kpt.P_ani.items(): 

51 Ra = pos_av[a] 

52 Pi_i = P_ni[state1] 

53 Pj_i = P_ni[state2] 

54 Delta_pL = wfs.setups[a].Delta_pL 

55 ni = len(Pi_i) 

56 

57 ma0 = 0 

58 ma1 = np.zeros(me.shape) 

59 for i in range(ni): 

60 for j in range(ni): 

61 pij = Pi_i[i] * Pj_i[j] 

62 ij = packed_index(i, j, ni) 

63 # L=0 term 

64 ma0 += Delta_pL[ij, 0] * pij 

65 # L=1 terms 

66 if wfs.setups[a].lmax >= 1: 

67 # see spherical_harmonics.py for 

68 # L=1:y L=2:z; L=3:x 

69 ma1 += np.array([ 

70 Delta_pL[ij, 3], Delta_pL[ij, 1], Delta_pL[ij, 2] 

71 ]) * pij 

72 ma += sqrt(4 * pi / 3) * ma1 + Ra * sqrt(4 * pi) * ma0 

73 gd.comm.sum(ma) 

74 

75 me += ma 

76 

77 return me * Bohr 

78 

79 # Currently only works on a single processor 

80 assert size == 1 

81 

82 maxfield = 0.01 

83 nfs = 5 # number of field 

84 nbands = 30 # number of bands 

85 h = 0.20 # grid spacing 

86 

87 debug = not False 

88 

89 if debug: 

90 txt = 'gpaw.out' 

91 else: 

92 txt = None 

93 

94 test1 = True 

95 test2 = True 

96 

97 a0 = 6.0 

98 a = Atoms('H', positions=[[a0 / 2, a0 / 2, a0 / 2]], cell=[a0, a0, a0]) 

99 

100 alpha1 = None 

101 alpha2 = None 

102 

103 # Test 1 

104 

105 if test1: 

106 c = GPAW(mode='fd', 

107 h=h, 

108 nbands=nbands + 10, 

109 spinpol=True, 

110 hund=True, 

111 xc='LDA', 

112 eigensolver='cg', 

113 convergence={ 

114 'bands': nbands, 

115 'eigenstates': 3.3e-4 

116 }, 

117 maxiter=1000, 

118 txt=txt) 

119 a.calc = c 

120 a.get_potential_energy() 

121 

122 o1 = c.get_occupation_numbers(spin=0) 

123 

124 if o1[0] > 0.0: 

125 spin = 0 

126 else: 

127 spin = 1 

128 

129 alpha = 0.0 

130 

131 ev = c.get_eigenvalues(0, spin) 

132 for i in range(1, nbands): 

133 mu_x, mu_y, mu_z = dipole_op(c, 0, i, k=0, s=spin) 

134 

135 alpha += mu_z**2 / (ev[i] - ev[0]) 

136 

137 alpha *= 2 

138 

139 if rank == 0 and debug: 

140 print('From perturbation theory:') 

141 print(' alpha = ', alpha, ' A**2/eV') 

142 print(' alpha = ', alpha * to_au, ' Bohr**3') 

143 

144 alpha1 = alpha 

145 

146 ### 

147 

148 c = GPAW( 

149 mode='fd', 

150 h=h, 

151 nbands=2, 

152 spinpol=True, 

153 hund=True, 

154 xc='LDA', 

155 # eigensolver = 'cg', 

156 convergence={ 

157 'bands': nbands, 

158 'eigenstates': 3.3e-4 

159 }, 

160 maxiter=1000, 

161 txt=txt) 

162 a.calc = c 

163 

164 # Test 2 

165 

166 if test2: 

167 e = [] 

168 e1s = [] 

169 d = [] 

170 fields = np.linspace(-maxfield, maxfield, nfs) 

171 for field in fields: 

172 if rank == 0 and debug: 

173 print(field) 

174 c = c.new(external=ConstantElectricField(field)) 

175 a.calc = c 

176 etot = a.get_potential_energy() 

177 e += [etot] 

178 ev0 = c.get_eigenvalues(0) 

179 ev1 = c.get_eigenvalues(0, 1) 

180 e1s += [min(ev0[0], ev1[0])] 

181 dip = c.get_dipole_moment() 

182 d += [dip[2]] 

183 

184 pol1, dummy = np.polyfit(fields, d, 1) 

185 pol2, dummy1, dummy2 = np.polyfit(fields, e1s, 2) 

186 

187 if rank == 0 and debug: 

188 print('From shift in 1s-state at constant electric field:') 

189 print(' alpha = ', -pol2, ' A**2/eV') 

190 print(' alpha = ', -pol2 * to_au, ' Bohr**3') 

191 

192 print('From dipole moment at constant electric field:') 

193 print(' alpha = ', pol1, ' A**2/eV') 

194 print(' alpha = ', pol1 * to_au, ' Bohr**3') 

195 

196 np.savetxt('ecf.out', np.transpose([fields, e, e1s, d])) 

197 

198 assert abs(pol1 + pol2) < 0.002 

199 alpha2 = (pol1 - pol2) / 2 

200 

201 # # This is a very, very rough test 

202 assert alpha1 == pytest.approx(alpha2, abs=0.01)