Coverage for gpaw/lrtddft/excitation.py: 75%

1"""Excitation lists base classes 

2 

3""" 

4from math import sqrt 

5import numpy as np 

6 

7from ase.units import Ha 

8 

9import gpaw.mpi as mpi 

10from gpaw.io.logger import GPAWLogger 

11 

12 

13def get_filehandle(cls, filename, mode='r'): 

14 cls.filename = filename 

15 if filename.endswith('.gz'): 

16 try: 

17 import gzip 

18 return gzip.open(filename, mode + 't') 

19 except ModuleNotFoundError: 

20 pass 

21 return open(filename, mode) 

22 

23 

24class ExcitationLogger(GPAWLogger): 

25 def header(self): 

26 pass 

27 

28 

29class ExcitationList(list): 

30 """General Excitation List class. 

31 """ 

32 def __init__(self, log=None, txt='-'): 

33 # initialise empty list 

34 list.__init__(self) 

35 self.energy_to_eV_scale = Ha 

36 

37 # set output 

38 if log is not None: 

39 self.log = log 

40 else: 

41 self.log = ExcitationLogger(world=mpi.world) 

42 self.log.fd = txt 

43 

44 @property 

45 def calc(self): 

46 1 / 0 

47 

48 def get_energies(self): 

49 """Get excitation energies in Hartrees""" 

50 el = [] 

51 for ex in self: 

52 el.append(ex.get_energy()) 

53 return np.array(el) 

54 

55 def get_trk(self): 

56 """Evaluate the Thomas Reiche Kuhn sum rule""" 

57 trkm = np.zeros(3) 

58 for ex in self: 

59 me = ex.get_dipole_me() 

60 trkm += ex.get_energy() * (me.real ** 2 + me.imag ** 2) 

61 return 2. * trkm # scale to get the number of electrons XXX spinpol ? 

62 

63 def get_polarizabilities(self, lmax=7): 

64 """Calculate the Polarisabilities 

65 see Jamorski et al. J. Chem. Phys. 104 (1996) 5134""" 

66 S = np.zeros(lmax + 1) 

67 for ex in self: 

68 e = ex.get_energy() 

69 f = ex.get_oscillator_strength()[0] 

70 for l in range(lmax + 1): 

71 S[l] += e ** (-2 * l) * f 

72 return S 

73 

74 def __truediv__(self, x): 

75 return self.__mul__(1. / x) 

76 

77 __div__ = __truediv__ 

78 

79 def __rmul__(self, x): 

80 return self.__mul__(x) 

81 

82 def __mul__(self, x): 

83 """Multiply with a number""" 

84 if isinstance(x, (float, int)): 

85 result = self.__class__() 

86 result.dtype = self.dtype 

87 for kss in self: 

88 result.append(x * kss) 

89 return result 

90 else: 

91 return RuntimeError('not a number') 

92 

93 def __sub__(self, other): 

94 result = self.__class__() 

95 result.dtype = self.dtype 

96 assert len(self) == len(other) 

97 for kss, ksso in zip(self, other): 

98 result.append(kss - ksso) 

99 return result 

100 

101 def __str__(self): 

102 string = '# ' + str(type(self)) 

103 if len(self) != 0: 

104 string += ', %d excitations:' % len(self) 

105 string += '\n' 

106 for ex in self: 

107 string += '# ' + ex.__str__() + '\n' 

108 return string 

109 

110 

111class Excitation: 

112 

113 def get_energy(self): 

114 """Get the excitations energy relative to the ground state energy 

115 in Hartrees. 

116 """ 

117 return self.energy 

118 

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

120 """return the excitations dipole matrix element 

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

122 if form == 'r': 

123 # length form 

124 return self.me / sqrt(self.energy) 

125 elif form == 'v': 

126 # velocity form 

127 return - np.sqrt(self.fij) * self.muv 

128 else: 

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

130 

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

132 """Return the "oscillator strength tensor" 

133 

134 self.me is assumed to be:: 

135 

136 form='r': sqrt(f * E) * <I|r|J>, 

137 form='v': sqrt(f / E) * <I|d/(dr)|J> 

138 

139 for f = multiplicity, E = transition energy and initial and 

140 final states:: 

141 

142 |I>, |J> 

143 """ 

144 

145 if form == 'r': 

146 # length form 

147 me = self.me 

148 elif form == 'v': 

149 # velocity form 

150 me = self.muv * np.sqrt(self.fij * self.energy) 

151 else: 

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

153 

154 return 2 * np.outer(me, me.conj()) 

155 

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

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

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

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

160 

161 def get_rotatory_strength(self, form='r', units='cgs'): 

162 """Return rotatory strength""" 

163 if self.magn is None: 

164 raise RuntimeError('Magnetic moment not available.') 

165 

166 if units == 'cgs': 

167 # 10^-40 esu cm erg / G 

168 # = 3.33564095 * 10^-15 A^2 m^3 s 

169 # conversion factor after 

170 # T. B. Pedersen and A. E. Hansen, 

171 # Chem. Phys. Lett. 246 (1995) 1 

172 # pre = 471.43 

173 # From TurboMole 

174 pre = 64604.8164 

175 elif units == 'a.u.': 

176 pre = 1. 

177 else: 

178 raise RuntimeError('Unknown units >' + units + '<') 

179 

180 if form == 'r': 

181 # length form 

182 mu = self.mur 

183 elif form == 'v': 

184 # velocity form 

185 mu = self.muv 

186 else: 

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

188 

189 return -pre * np.dot(mu, self.magn) 

190 

191 def set_energy(self, E): 

192 """Set the excitations energy relative to the ground state energy""" 

193 self.energy = E