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

1# fmt: off 

2 

3import numpy as np 

4 

5from ase.units import Bohr, Hartree 

6 

7 

8class Excitation: 

9 """Base class for a single excitation""" 

10 

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

12 """ 

13 Parameters 

14 ---------- 

15 energy: float 

16 Energy realtive to the ground state 

17 index: int 

18 Excited state index 

19 mur: list of three floats or complex numbers 

20 Length form dipole matrix element 

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

22 Velocity form dipole matrix element, default None 

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

24 Magnetic matrix element, default None 

25 """ 

26 self.energy = energy 

27 self.index = index 

28 self.mur = mur 

29 self.muv = muv 

30 self.magn = magn 

31 self.fij = 1. 

32 

33 def outstring(self): 

34 """Format yourself as a string""" 

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

36 

37 def format_me(me): 

38 string = '' 

39 if me.dtype == float: 

40 for m in me: 

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

42 else: 

43 for m in me: 

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

45 return string 

46 

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

48 if self.muv is not None: 

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

50 if self.magn is not None: 

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

52 string += '\n' 

53 

54 return string 

55 

56 @classmethod 

57 def fromstring(cls, string): 

58 """Initialize yourself from a string""" 

59 l = string.split() 

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

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

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

63 try: 

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

65 except IndexError: 

66 muv = None 

67 try: 

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

69 except IndexError: 

70 magn = None 

71 

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

73 

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

75 """Return the excitations dipole matrix element 

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

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

78 me = - self.mur 

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

80 me = - self.muv 

81 else: 

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

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

84 

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

86 """Return the oscillator strength tensor""" 

87 me = self.get_dipole_me(form) 

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

89 

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

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

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

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

94 

95 

96class ExcitationList(list): 

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

98 

99 def __init__(self): 

100 # initialise empty list 

101 super().__init__() 

102 

103 # set default energy scale to get eV 

104 self.energy_to_eV_scale = 1. 

105 

106 

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

108 tensor=False, index=0): 

109 """Evaluate the photon energy dependent polarizability 

110 from the sum over states 

111 

112 Parameters 

113 ---------- 

114 exlist: ExcitationList 

115 omega: 

116 Photon energy (eV) 

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

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

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

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

121 tensor: boolean 

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

123 index is ignored, default False 

124 

125 Returns 

126 ------- 

127 alpha: 

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

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

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

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

132 """ 

133 omega = np.asarray(omega) 

134 om2 = 1. * omega**2 

135 esc = exlist.energy_to_eV_scale 

136 

137 if tensor: 

138 if not np.isscalar(om2): 

139 om2 = om2[:, None, None] 

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

141 dtype=om2.dtype) 

142 for ex in exlist: 

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

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

145 else: 

146 alpha = np.zeros_like(om2) 

147 for ex in exlist: 

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

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

150 

151 return alpha * Bohr**2 * Hartree