Coverage for /builds/kinetik161/ase/ase/calculators/excitation_list.py: 90.41%

1import numpy as np 

2 

3from ase.units import Bohr, Hartree 

4 

5 

6class Excitation: 

7 """Base class for a single excitation""" 

8 

9 def __init__(self, energy, index, mur, muv=None, magn=None): 

10 """ 

11 Parameters 

12 ---------- 

13 energy: float 

14 Energy realtive to the ground state 

15 index: int 

16 Excited state index 

17 mur: list of three floats or complex numbers 

18 Length form dipole matrix element 

19 muv: list of three floats or complex numbers or None 

20 Velocity form dipole matrix element, default None 

21 magn: list of three floats or complex numbers or None 

22 Magnetic matrix element, default None 

23 """ 

24 self.energy = energy 

25 self.index = index 

26 self.mur = mur 

27 self.muv = muv 

28 self.magn = magn 

29 self.fij = 1. 

30 

31 def outstring(self): 

32 """Format yourself as a string""" 

33 string = f'{self.energy:g} {self.index} ' 

34 

35 def format_me(me): 

36 string = '' 

37 if me.dtype == float: 

38 for m in me: 

39 string += f' {m:g}' 

40 else: 

41 for m in me: 

42 string += ' {0.real:g}{0.imag:+g}j'.format(m) 

43 return string 

44 

45 string += ' ' + format_me(self.mur) 

46 if self.muv is not None: 

47 string += ' ' + format_me(self.muv) 

48 if self.magn is not None: 

49 string += ' ' + format_me(self.magn) 

50 string += '\n' 

51 

52 return string 

53 

54 @classmethod 

55 def fromstring(cls, string): 

56 """Initialize yourself from a string""" 

57 l = string.split() 

58 energy = float(l.pop(0)) 

59 index = int(l.pop(0)) 

60 mur = np.array([float(l.pop(0)) for i in range(3)]) 

61 try: 

62 muv = np.array([float(l.pop(0)) for i in range(3)]) 

63 except IndexError: 

64 muv = None 

65 try: 

66 magn = np.array([float(l.pop(0)) for i in range(3)]) 

67 except IndexError: 

68 magn = None 

69 

70 return cls(energy, index, mur, muv, magn) 

71 

72 def get_dipole_me(self, form='r'): 

73 """Return the excitations dipole matrix element 

74 including the occupation factor sqrt(fij)""" 

75 if form == 'r': # length form 

76 me = - self.mur 

77 elif form == 'v': # velocity form 

78 me = - self.muv 

79 else: 

80 raise RuntimeError('Unknown form >' + form + '<') 

81 return np.sqrt(self.fij) * me 

82 

83 def get_dipole_tensor(self, form='r'): 

84 """Return the oscillator strength tensor""" 

85 me = self.get_dipole_me(form) 

86 return 2 * np.outer(me, me.conj()) * self.energy / Hartree 

87 

88 def get_oscillator_strength(self, form='r'): 

89 """Return the excitations dipole oscillator strength.""" 

90 me2_c = self.get_dipole_tensor(form).diagonal().real 

91 return np.array([np.sum(me2_c) / 3.] + me2_c.tolist()) 

92 

93 

94class ExcitationList(list): 

95 """Base class for excitions from the ground state""" 

96 

97 def __init__(self): 

98 # initialise empty list 

99 super().__init__() 

100 

101 # set default energy scale to get eV 

102 self.energy_to_eV_scale = 1. 

103 

104 

105def polarizability(exlist, omega, form='v', 

106 tensor=False, index=0): 

107 """Evaluate the photon energy dependent polarizability 

108 from the sum over states 

109 

110 Parameters 

111 ---------- 

112 exlist: ExcitationList 

113 omega: 

114 Photon energy (eV) 

115 form: {'v', 'r'} 

116 Form of the dipole matrix element, default 'v' 

117 index: {0, 1, 2, 3} 

118 0: averaged, 1,2,3:alpha_xx, alpha_yy, alpha_zz, default 0 

119 tensor: boolean 

120 if True returns alpha_ij, i,j=x,y,z 

121 index is ignored, default False 

122 

123 Returns 

124 ------- 

125 alpha: 

126 Unit (e^2 Angstrom^2 / eV). 

127 Multiply with Bohr * Ha to get (Angstrom^3) 

128 shape = (omega.shape,) if tensor == False 

129 shape = (omega.shape, 3, 3) else 

130 """ 

131 omega = np.asarray(omega) 

132 om2 = 1. * omega**2 

133 esc = exlist.energy_to_eV_scale 

134 

135 if tensor: 

136 if not np.isscalar(om2): 

137 om2 = om2[:, None, None] 

138 alpha = np.zeros(omega.shape + (3, 3), 

139 dtype=om2.dtype) 

140 for ex in exlist: 

141 alpha += ex.get_dipole_tensor(form=form) / ( 

142 (ex.energy * esc)**2 - om2) 

143 else: 

144 alpha = np.zeros_like(om2) 

145 for ex in exlist: 

146 alpha += ex.get_oscillator_strength(form=form)[index] / ( 

147 (ex.energy * esc)**2 - om2) 

148 

149 return alpha * Bohr**2 * Hartree