Coverage for gpaw/lrtddft/spectrum.py: 57%

1import sys 

2import numpy as np 

3 

4from ase.units import _hbar, _c, _e, _me, Hartree, Bohr 

5from gpaw import __version__ as version 

6from gpaw.utilities.folder import Folder 

7 

8 

9def get_folded_spectrum( 

10 exlist=None, 

11 emin=None, 

12 emax=None, 

13 de=None, 

14 energyunit='eV', 

15 folding='Gauss', 

16 width=0.08, # Gauss/Lorentz width 

17 form='r'): 

18 """Return folded spectrum.""" 

19 

20 x = [] 

21 y = [] 

22 for ex in exlist: 

23 x.append(ex.get_energy() * Hartree) 

24 y.append(ex.get_oscillator_strength(form)) 

25 

26 if energyunit == 'nm': 

27 # transform to experimentally used wavelength [nm] 

28 x = 1.e+9 * 2 * np.pi * _hbar * _c / _e / np.array(x) 

29 elif energyunit != 'eV': 

30 raise RuntimeError('currently only eV and nm are supported') 

31 

32 if folding is None: 

33 return x, y 

34 else: 

35 return Folder(width, folding).fold(x, y, de, emin, emax) 

36 

37 

38def spectrum(exlist=None, 

39 filename=None, 

40 emin=None, 

41 emax=None, 

42 de=None, 

43 energyunit='eV', 

44 folding='Gauss', 

45 width=0.08, # Gauss/Lorentz width 

46 comment=None, 

47 form='r'): 

48 """Write out a folded spectrum. 

49 

50 Parameters 

51 ---------- 

52 exlist: ExcitationList 

53 filename: 

54 File name for the output file, STDOUT if not given 

55 emin: 

56 min. energy, set to cover all energies if not given 

57 emax: 

58 max. energy, set to cover all energies if not given 

59 de: 

60 energy spacing 

61 energyunit: {'eV', 'nm'} 

62 Energy unit, default 'eV' 

63 folding: {'Gauss', 'Lorentz'} 

64 width: 

65 folding width in terms of the chosen energyunit 

66 """ 

67 

68 # output 

69 out = sys.stdout 

70 if filename is not None: 

71 with open(filename, 'w') as out: 

72 if comment: 

73 print('#', comment, file=out) 

74 

75 print('# Photoabsorption spectrum from linear response TD-DFT', 

76 file=out) 

77 print('# GPAW version:', version, file=out) 

78 if folding is not None: # fold the spectrum 

79 print(f'# {folding} folded, width={width:g} [{energyunit}]', 

80 file=out) 

81 if form == 'r': 

82 out.write('# length form') 

83 else: 

84 assert form == 'v' 

85 out.write('# velocity form') 

86 print('# om [%s] osz osz x osz y osz z' 

87 % energyunit, file=out) 

88 

89 energies, values = get_folded_spectrum(exlist, emin, emax, de, 

90 energyunit, folding, 

91 width, form) 

92 

93 for e, val in zip(energies, values): 

94 print('%10.5f %12.7e %12.7e %11.7e %11.7e' % 

95 (e, val[0], val[1], val[2], val[3]), file=out) 

96 

97 

98def get_adsorbance_pre_factor(atoms): 

99 """Return the absorbance pre-factor for solids. Unit m^-1. 

100 

101 Use this factor to multiply the folded oscillatior strength 

102 obtained from spectrum() 

103 

104 Robert C. Hilborn, Am. J. Phys. 50, 982 (1982) 

105 """ 

106 return np.pi * _e**2 / 2. / _me / _c 

107 

108 

109def rotatory_spectrum(exlist=None, 

110 filename=None, 

111 emin=None, 

112 emax=None, 

113 de=None, 

114 energyunit='eV', 

115 folding='Gauss', 

116 width=0.08, # Gauss/Lorentz width 

117 comment=None): 

118 """Write out a folded rotatory spectrum. 

119 

120 See spectrum() for explanation of the parameters. 

121 """ 

122 

123 # output 

124 out = sys.stdout 

125 if filename is not None: 

126 with open(filename, 'w') as out: 

127 if comment: 

128 print('#', comment, file=out) 

129 

130 print('# Rotatory spectrum from linear response TD-DFT', file=out) 

131 print('# GPAW version:', version, file=out) 

132 if folding is not None: # fold the spectrum 

133 print(f'# {folding} folded, width={width:g} [{energyunit}]', 

134 file=out) 

135 print(f'# om [{energyunit}] R [cgs]', file=out) 

136 

137 x = [] 

138 y = [] 

139 for ex in exlist: 

140 x.append(ex.get_energy() * Hartree) 

141 y.append(ex.get_rotatory_strength()) 

142 

143 if energyunit == 'nm': 

144 # transform to experimentally used wavelength [nm] 

145 x = 1.e+9 * 2 * np.pi * _hbar * _c / _e / np.array(x) 

146 y = np.array(y) 

147 elif energyunit != 'eV': 

148 raise RuntimeError('currently only eV and nm are supported') 

149 

150 if folding is None: 

151 energies, values = x, y 

152 else: 

153 energies, values = Folder(width, folding).fold(x, y, de, 

154 emin, emax) 

155 for e, val in zip(energies, values): 

156 print(f'{e:10.5f} {val:12.7e}', file=out) 

157 

158 

159class Writer(Folder): 

160 """Object to write data to a file""" 

161 def __init__(self, folding=None, width=0.08): 

162 """Evaluate the polarizability from sum over states. 

163 

164 Parameters 

165 ---------- 

166 folding : string or None(default) 

167 Name of the folding function, see gpaw/folder.py for options. 

168 None means do not fold at all. 

169 width : float 

170 The width to be transferred to the folding function. 

171 """ 

172 self.folding = folding 

173 if folding is not None: 

174 Folder.__init__(self, width, folding) 

175 

176 def write(self, filename=None, 

177 emin=None, emax=None, de=None, 

178 comment=None): 

179 

180 out = sys.stdout 

181 if filename is not None: 

182 with open(filename, 'w') as out: 

183 print('#', self.title, file=out) 

184 print('# GPAW version:', version, file=out) 

185 if comment: 

186 print('#', comment, file=out) 

187 if self.folding is not None: 

188 print('# {} folded, width={:g} [eV]'.format(self.folding, 

189 self.width), file=out) 

190 energies, values = self.fold(self.energies, self.values, 

191 de, emin, emax) 

192 else: 

193 energies, values = self.energies, self.values 

194 

195 print('#', self.fields, file=out) 

196 

197 for e, val in zip(energies, values): 

198 string = '%10.5f' % e 

199 for vf in val: 

200 string += ' %12.7e' % vf 

201 print(string, file=out) 

202 

203 

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

205 tensor=False, index=0): 

206 """Evaluate the polarizability from sum over states. 

207 

208 Parameters 

209 ---------- 

210 exlist: ExcitationList 

211 omega: 

212 Photon energy (eV) 

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

214 Form of the dipole matrix element 

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

216 0: averaged, 1,2,3:alpha_xx, alpha_yy, alpha_zz 

217 tensor: 

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

219 index is ignored 

220 

221 Returns 

222 ------- 

223 alpha: 

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

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

226 """ 

227 if tensor: 

228 alpha = np.zeros(np.array(omega).shape + (3, 3), dtype=complex) 

229 for ex in exlist: 

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

231 (ex.energy * Hartree)**2 - omega**2) 

232 else: 

233 alpha = np.zeros_like(omega) 

234 for ex in exlist: 

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

236 (ex.energy * Hartree)**2 - omega**2) 

237 

238 return alpha * Bohr**2 * Hartree