Coverage for gpaw/pes/state.py: 96%

1from math import exp, pi, sqrt 

2 

3import numpy as np 

4 

5from ase.units import Bohr, Hartree, alpha 

6from gpaw.pes import ds_prefactor 

7 

8 

9class State: 

10 """Electronic state base class.""" 

11 

12 def __init__(self, gd): 

13 self.gd = gd 

14 self.r_c = [] 

15 for c in range(3): 

16 r_G = (np.arange(gd.n_c[c], dtype=float) + 

17 gd.beg_c[c]) * gd.h_cv.diagonal()[c] 

18 self.r_c.append(r_G) 

19 

20 def get_grid(self): 

21 return self.grid_g 

22 

23 def set_grid(self, grid_g): 

24 # assert(grid_g.shape == self.gd.empty().shape) 

25 self.grid_g = grid_g 

26 

27 def get_energy(self): 

28 """The states energy in [eV]""" 

29 return self.energy * Hartree 

30 

31 def set_energy(self, energy): 

32 """Set the states energy in [eV]""" 

33 self.energy = energy / Hartree 

34 

35 

36class BoundState(State): 

37 

38 """Bound state originating from a gpaw calculation.""" 

39 

40 def __init__(self, 

41 calculator, 

42 s, n): 

43 self.gd = calculator.wfs.gd 

44 self.kpt = calculator.wfs.kpt_u[s] 

45 self.grid_g = self.kpt.psit_nG[n] 

46 self.energy = self.kpt.eps_n[n] 

47 

48 

49class H1s(State): 

50 

51 """Analytic H1s state.""" 

52 

53 def __init__(self, gd=None, center=None, Z=1): 

54 if center is None: 

55 center = 0.5 * gd.n_c * gd.h_c 

56 else: 

57 self.center = center / Bohr 

58 self.Z = Z 

59 self.energy = - Z ** 2 / 2. 

60 if gd is not None: 

61 self.grid_g = None 

62 State.__init__(self, gd) 

63 

64 def get_grid(self): 

65 if self.grid_g is None: 

66 gd = self.gd 

67 assert gd.orthogonal 

68 wf = gd.zeros(dtype=float) 

69 vr0 = np.array([0., 0., 0.]) 

70 for i in range(gd.n_c[0]): 

71 vr0[0] = (gd.beg_c[0] + i) * gd.h_cv[0, 0] 

72 for j in range(gd.n_c[1]): 

73 vr0[1] = (gd.beg_c[1] + j) * gd.h_cv[1, 1] 

74 for k in range(gd.n_c[2]): 

75 vr0[2] = (gd.beg_c[2] + k) * gd.h_cv[2, 2] 

76 vr = vr0 - self.center 

77 r = sqrt(np.dot(vr, vr)) 

78 wf[i, j, k] = exp(-self.Z * r) 

79 self.grid_g = wf * sqrt(self.Z ** 3 / pi) 

80 return self.grid_g 

81 

82 def get_me_c(self, k_c, form='L'): 

83 """Transition matrix element.""" 

84 k = sqrt(np.dot(k_c, k_c)) 

85 if form == 'L': 

86 pre = 2j 

87 elif form == 'V': 

88 pre = 1j 

89 else: 

90 raise NotImplementedError 

91 return pre * k_c * self.FT(k) 

92 

93 def get_ds(self, Ekin, form='L', units='Mb'): 

94 """Angular averaged cross section. 

95 

96 Ekin: photoelectron kinetic energy [eV]""" 

97 E = Ekin / Hartree 

98 k = sqrt(2 * E) 

99 omega = E - self.energy 

100 k_c = np.array([0., 0., k]) 

101 me_c = self.get_me_c(k_c, form) 

102 T2 = np.abs(np.dot(me_c, me_c)) / 3. 

103 pre = ds_prefactor[units] 

104 return pre * (2 * pi) ** 2 * alpha * k * 4 * pi / omega * T2 

105 

106 def FT(self, k): 

107 return sqrt(8 * self.Z ** 5) / pi / (k ** 2 + self.Z ** 2) ** 2