1import numpy as np 

2 

3from ase.units import Hartree, Bohr 

4 

5 

6class Excitation: 

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

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

9 """ 

10 Parameters 

11 ---------- 

12 energy: float 

13 Energy realtive to the ground state 

14 index: int 

15 Excited state index 

16 mur: list of three floats or complex numbers 

17 Length form dipole matrix element 

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

19 Velocity form dipole matrix element, default None 

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

21 Magnetic matrix element, default None 

22 """ 

23 self.energy = energy 

24 self.index = index 

25 self.mur = mur 

26 self.muv = muv 

27 self.magn = magn 

28 self.fij = 1. 

29 

30 def outstring(self): 

31 """Format yourself as a string""" 

32 string = '{0:g} {1} '.format(self.energy, self.index) 

33 

34 def format_me(me): 

35 string = '' 

36 if me.dtype == float: 

37 for m in me: 

38 string += ' {0:g}'.format(m) 

39 else: 

40 for m in me: 

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

42 return string 

43 

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

45 if self.muv is not None: 

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

47 if self.magn is not None: 

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

49 string += '\n' 

50 

51 return string 

52 

53 @classmethod 

54 def fromstring(cls, string): 

55 """Initialize yourself from a string""" 

56 l = string.split() 

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

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

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

60 try: 

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

62 except IndexError: 

63 muv = None 

64 try: 

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

66 except IndexError: 

67 magn = None 

68 

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

70 

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

72 """Return the excitations dipole matrix element 

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

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

75 me = - self.mur 

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

77 me = - self.muv 

78 else: 

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

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

81 

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

83 """Return the oscillator strength tensor""" 

84 me = self.get_dipole_me(form) 

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

86 

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

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

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

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

91 

92 

93class ExcitationList(list): 

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

95 def __init__(self): 

96 # initialise empty list 

97 super().__init__() 

98 

99 # set default energy scale to get eV 

100 self.energy_to_eV_scale = 1. 

101 

102 

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

104 tensor=False, index=0): 

105 """Evaluate the photon energy dependent polarizability 

106 from the sum over states 

107 

108 Parameters 

109 ---------- 

110 exlist: ExcitationList 

111 omega: 

112 Photon energy (eV) 

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

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

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

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

117 tensor: boolean 

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

119 index is ignored, default False 

120 

121 Returns 

122 ------- 

123 alpha: 

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

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

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

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

128 """ 

129 omega = np.asarray(omega) 

130 om2 = 1. * omega**2 

131 esc = exlist.energy_to_eV_scale 

132 

133 if tensor: 

134 if not np.isscalar(om2): 

135 om2 = om2[:, None, None] 

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

137 dtype=om2.dtype) 

138 for ex in exlist: 

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

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

141 else: 

142 alpha = np.zeros_like(om2) 

143 for ex in exlist: 

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

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

146 

147 return alpha * Bohr**2 * Hartree