Coverage for gpaw/pes/ds_beta.py: 94%

1from math import sin, cos, pi, sqrt 

2 

3import numpy as np 

4 

5from ase.units import Bohr, Hartree, alpha 

6from gpaw.fd_operators import Gradient 

7from gpaw.utilities.gl_quadrature import GaussLegendre 

8from gpaw.pes import ds_prefactor 

9from gpaw.pes.state import H1s 

10from gpaw.pes.continuum import PlaneWave 

11 

12# debug 

13# from gpaw.mpi import rank 

14 

15 

16class CrossSectionBeta: 

17 

18 def __init__(self, 

19 initial=None, 

20 final=None, 

21 r0=[0, 0, 0], # center of mass vector 

22 form='L', 

23 ngauss=8): 

24 

25 self.initial = initial 

26 self.r0 = np.array(r0) / Bohr 

27 self.form = form 

28 if final is None: 

29 self.final = PlaneWave(initial.gd) 

30 else: 

31 self.final = final 

32 

33 self.Ekin = None 

34 

35 # Gauss-Legendre weights and abscissas 

36 self.gl = {} 

37 self.gl['x'] = GaussLegendre(-1., 1., ngauss) 

38 self.gl['phi'] = GaussLegendre(0, 2 * pi, ngauss) 

39 self.gl['psi'] = self.gl['phi'] 

40 self.angle = {} 

41 

42 # sin and cos of the magic angle (54.7 deg) 

43 self.costhm = 1. / sqrt(3) 

44 self.sinthm = sqrt(2. / 3.) 

45 

46 def calculate(self, Ekin): 

47 """Calculate the necessary overlaps.""" 

48 

49 Ekin = Ekin / Hartree 

50 if self.Ekin == Ekin: 

51 return 

52 self.Ekin = Ekin 

53 

54 # photoelectron momentum 

55 self.k = sqrt(2 * self.Ekin) 

56 

57 for angle in ['x', 'phi', 'psi']: 

58 self.angle[angle] = self.gl[angle].get_x()[0] 

59 self.T20, self.T2m = self.gi_x() 

60 

61 # we need the average 

62 self.T20 /= 8 * pi ** 2 

63 self.T2m /= 8 * pi ** 2 

64 

65 def get_omega(self): 

66 """Return the necessary photon energy.""" 

67 return self.Ekin - self.initial.get_energy() / Hartree 

68 

69 def get_beta(self, Ekin=None): 

70 """Return the asymmetry parameter. 

71 

72 E: photoelectron kinetic energy [eV] 

73 """ 

74 if Ekin is not None: 

75 self.calculate(Ekin) 

76 return self.T20 / self.T2m - 1. 

77 

78 def get_ds(self, Ekin=None, units='Mb'): 

79 """Return the total cross section. 

80 

81 Ekin: photoelectron kinetic energy [eV] 

82 units: 

83 'Mb', 1 Mb = 1.e-22 m**2 

84 'Ang', 1 A**2 = 1.e-20 m**2 

85 'a.u.', 1 a_0**2 = 2.8e-21 m**2 

86 as output units 

87 """ 

88 if Ekin is not None: 

89 self.calculate(Ekin) 

90 try: 

91 pre = ds_prefactor[units] 

92 except KeyError: 

93 print('Unknown units: >' + units + '<') 

94 raise 

95 

96# me_c = self.initial.get_me_c(np.array([0., 0., self.k]), self.form) 

97# T2mana = np.abs(np.dot(me_c,me_c)) / 3. 

98# print "T2m:", T2mana, self.T2m 

99 

100 omega = self.get_omega() 

101 

102 # integration over momentum agles 

103 pre *= self.k * 4 * pi 

104 

105# print omega, self.initial.get_ds(self.k, omega, self.form), \ 

106# (self.k * 4 * pi * (2 * pi)**2 / 137.0359895 * self.T2m / omega) 

107 

108 return pre * ((2 * pi) ** 2 * alpha * self.T2m / omega) 

109 

110 def gauss_integrate(self, angle, function): 

111 T20 = 0. 

112 T2m = 0. 

113 gl = self.gl[angle] 

114 for x, w in zip(gl.get_x(), gl.get_w()): 

115 self.angle[angle] = x 

116 t20, t2m = function() 

117 T20 += w * t20 

118 T2m += w * t2m 

119# print "<gauss_integrate> angle=", angle, T2m 

120 return T20, T2m 

121 

122 def gi_x(self): 

123 """Gauss integrate over x=cos(theta)""" 

124 return self.gauss_integrate('x', self.gi_phi) 

125 

126 def gi_phi(self): 

127 """Gauss integrate over phi""" 

128 return self.gauss_integrate('phi', self.gi_psi) 

129 

130 def gi_psi(self): 

131 """Gauss integrate over psi""" 

132 return self.gauss_integrate('psi', self.integrand) 

133 

134 def integrand(self): 

135 

136 # polarisation in the direction of vk 

137 costh = self.angle['x'] 

138 sinth = sqrt(1. - costh ** 2) 

139 sinphi = sin(self.angle['phi']) 

140 cosphi = cos(self.angle['phi']) 

141 eps0 = np.array([sinth * cosphi, 

142 sinth * sinphi, 

143 costh]) 

144 vk = self.k * eps0 

145 

146 # polarisation at the magic angle 

147 costhm = self.costhm 

148 sinthm = self.sinthm 

149 sinpsi = sin(self.angle['psi']) 

150 cospsi = cos(self.angle['psi']) 

151 epsm = np.array([sinthm * (cosphi * sinpsi * costh + 

152 sinphi * cospsi) + 

153 costhm * cosphi * sinth, 

154 sinthm * (sinphi * sinpsi * costh - 

155 cosphi * cospsi) + 

156 costhm * sinphi * sinth, 

157 costhm * costh - sinthm * sinth * sinpsi]) 

158 

159 # initial and final state on the grid 

160 initial_G = self.initial.get_grid() 

161 final_G = self.final.get_grid(vk, self.r0) 

162 ini_analyt = H1s(self.initial.gd, self.r0) 

163 

164 gd = self.initial.gd 

165 if self.form == 'L': 

166 if_G = initial_G * final_G 

167 omega = self.get_omega() 

168 if 0: 

169 me_c = [] 

170 for c in range(3): 

171 xyz_G = ((np.arange(gd.n_c[c]) + gd.beg_c[c]) * gd.h_c[c] 

172 - self.r0[c]) 

173 shape = [1, 1, 1] 

174 shape[c] = -1 

175 xyz_G.shape = tuple(shape) 

176 np.resize(xyz_G, gd.n_c) 

177 me_c.append(gd.integrate(if_G * xyz_G)) 

178 me_c = np.array(me_c) * omega 

179 else: 

180 me_c = gd.calculate_dipole_moment(if_G) 

181 me_c += self.r0 * gd.integrate(if_G) 

182 me_c *= -omega 

183 elif self.form == 'V': 

184 dtype = final_G.dtype 

185 phase_cd = np.ones((3, 2), dtype) 

186 if not hasattr(gd, 'ddr'): 

187 gd.ddr = [Gradient(gd, c, dtype=dtype).apply for c in range(3)] 

188 dfinal_G = gd.empty(dtype=dtype) 

189 me_c = np.empty(3, dtype=dtype) 

190 for c in range(3): 

191 gd.ddr[c](final_G, dfinal_G, phase_cd) 

192 me_c[c] = gd.integrate(initial_G * dfinal_G) 

193 else: 

194 raise NotImplementedError 

195 

196 if 0: 

197 omega = self.get_omega() 

198 me_analyt = ini_analyt.get_me_c(vk, self.form)[0].imag 

199 me = me_c[0].imag 

200 

201 def ds(me): 

202 return self.k / omega * me ** 2 

203 

204 print(omega, ds(me_analyt), ds(me), me_analyt, me) 

205# print 'analyt', self.initial.get_me_c(vk, self.form) 

206# print 'num', me_c 

207# print 'analyt/num', self.initial.get_me_c(vk, self.form) / me_c 

208 

209 # return the squared matrix elements 

210 T2 = [] 

211 for eps in [eps0, epsm]: 

212 me = np.dot(eps, me_c) 

213# print "eps, T2:", eps, (me * me.conj()).real 

214 T2.append((me * me.conj()).real) 

215# print vk, T2 

216 return T2[0], T2[1]