1"""Infrared intensities""" 

2 

3from math import sqrt 

4from sys import stdout 

5 

6import numpy as np 

7 

8import ase.units as units 

9from ase.parallel import parprint, paropen 

10from ase.vibrations import Vibrations 

11 

12 

13class Infrared(Vibrations): 

14 """Class for calculating vibrational modes and infrared intensities 

15 using finite difference. 

16 

17 The vibrational modes are calculated from a finite difference 

18 approximation of the Dynamical matrix and the IR intensities from 

19 a finite difference approximation of the gradient of the dipole 

20 moment. The method is described in: 

21 

22 D. Porezag, M. R. Pederson: 

23 "Infrared intensities and Raman-scattering activities within 

24 density-functional theory", 

25 Phys. Rev. B 54, 7830 (1996) 

26 

27 The calculator object (calc) linked to the Atoms object (atoms) must 

28 have the attribute: 

29 

30 >>> calc.get_dipole_moment(atoms) 

31 

32 In addition to the methods included in the ``Vibrations`` class 

33 the ``Infrared`` class introduces two new methods; 

34 *get_spectrum()* and *write_spectra()*. The *summary()*, *get_energies()*, 

35 *get_frequencies()*, *get_spectrum()* and *write_spectra()* 

36 methods all take an optional *method* keyword. Use 

37 method='Frederiksen' to use the method described in: 

38 

39 T. Frederiksen, M. Paulsson, M. Brandbyge, A. P. Jauho: 

40 "Inelastic transport theory from first-principles: methodology 

41 and applications for nanoscale devices", 

42 Phys. Rev. B 75, 205413 (2007) 

43 

44 atoms: Atoms object 

45 The atoms to work on. 

46 indices: list of int 

47 List of indices of atoms to vibrate. Default behavior is 

48 to vibrate all atoms. 

49 name: str 

50 Name to use for files. 

51 delta: float 

52 Magnitude of displacements. 

53 nfree: int 

54 Number of displacements per degree of freedom, 2 or 4 are 

55 supported. Default is 2 which will displace each atom +delta 

56 and -delta in each cartesian direction. 

57 directions: list of int 

58 Cartesian coordinates to calculate the gradient 

59 of the dipole moment in. 

60 For example directions = 2 only dipole moment in the z-direction will 

61 be considered, whereas for directions = [0, 1] only the dipole 

62 moment in the xy-plane will be considered. Default behavior is to 

63 use the dipole moment in all directions. 

64 

65 Example: 

66 

67 >>> from ase.io import read 

68 >>> from ase.calculators.vasp import Vasp 

69 >>> from ase.vibrations import Infrared 

70 >>> water = read('water.traj') # read pre-relaxed structure of water 

71 >>> calc = Vasp(prec='Accurate', 

72 ... ediff=1E-8, 

73 ... isym=0, 

74 ... idipol=4, # calculate the total dipole moment 

75 ... dipol=water.get_center_of_mass(scaled=True), 

76 ... ldipol=True) 

77 >>> water.calc = calc 

78 >>> ir = Infrared(water) 

79 >>> ir.run() 

80 >>> ir.summary() 

81 ------------------------------------- 

82 Mode Frequency Intensity 

83 # meV cm^-1 (D/Å)^2 amu^-1 

84 ------------------------------------- 

85 0 16.9i 136.2i 1.6108 

86 1 10.5i 84.9i 2.1682 

87 2 5.1i 41.1i 1.7327 

88 3 0.3i 2.2i 0.0080 

89 4 2.4 19.0 0.1186 

90 5 15.3 123.5 1.4956 

91 6 195.5 1576.7 1.6437 

92 7 458.9 3701.3 0.0284 

93 8 473.0 3814.6 1.1812 

94 ------------------------------------- 

95 Zero-point energy: 0.573 eV 

96 Static dipole moment: 1.833 D 

97 Maximum force on atom in `equilibrium`: 0.0026 eV/Å 

98 

99 

100 

101 This interface now also works for calculator 'siesta', 

102 (added get_dipole_moment for siesta). 

103 

104 Example: 

105 

106 >>> #!/usr/bin/env python3 

107 

108 >>> from ase.io import read 

109 >>> from ase.calculators.siesta import Siesta 

110 >>> from ase.vibrations import Infrared 

111 

112 >>> bud = read('bud1.xyz') 

113 

114 >>> calc = Siesta(label='bud', 

115 ... meshcutoff=250 * Ry, 

116 ... basis='DZP', 

117 ... kpts=[1, 1, 1]) 

118 

119 >>> calc.set_fdf('DM.MixingWeight', 0.08) 

120 >>> calc.set_fdf('DM.NumberPulay', 3) 

121 >>> calc.set_fdf('DM.NumberKick', 20) 

122 >>> calc.set_fdf('DM.KickMixingWeight', 0.15) 

123 >>> calc.set_fdf('SolutionMethod', 'Diagon') 

124 >>> calc.set_fdf('MaxSCFIterations', 500) 

125 >>> calc.set_fdf('PAO.BasisType', 'split') 

126 >>> #50 meV = 0.003674931 * Ry 

127 >>> calc.set_fdf('PAO.EnergyShift', 0.003674931 * Ry ) 

128 >>> calc.set_fdf('LatticeConstant', 1.000000 * Ang) 

129 >>> calc.set_fdf('WriteCoorXmol', 'T') 

130 

131 >>> bud.calc = calc 

132 

133 >>> ir = Infrared(bud) 

134 >>> ir.run() 

135 >>> ir.summary() 

136 

137 """ 

138 

139 def __init__(self, atoms, indices=None, name='ir', delta=0.01, 

140 nfree=2, directions=None): 

141 Vibrations.__init__(self, atoms, indices=indices, name=name, 

142 delta=delta, nfree=nfree) 

143 if atoms.constraints: 

144 print('WARNING! \n Your Atoms object is constrained. ' 

145 'Some forces may be unintended set to zero. \n') 

146 if directions is None: 

147 self.directions = np.asarray([0, 1, 2]) 

148 else: 

149 self.directions = np.asarray(directions) 

150 self.ir = True 

151 

152 def read(self, method='standard', direction='central'): 

153 self.method = method.lower() 

154 self.direction = direction.lower() 

155 assert self.method in ['standard', 'frederiksen'] 

156 

157 if direction != 'central': 

158 raise NotImplementedError( 

159 'Only central difference is implemented at the moment.') 

160 

161 disp = self._eq_disp() 

162 forces_zero = disp.forces() 

163 dipole_zero = disp.dipole() 

164 self.dipole_zero = (sum(dipole_zero**2)**0.5) / units.Debye 

165 self.force_zero = max([sum((forces_zero[j])**2)**0.5 

166 for j in self.indices]) 

167 

168 ndof = 3 * len(self.indices) 

169 H = np.empty((ndof, ndof)) 

170 dpdx = np.empty((ndof, 3)) 

171 r = 0 

172 

173 for a, i in self._iter_ai(): 

174 disp_minus = self._disp(a, i, -1) 

175 disp_plus = self._disp(a, i, 1) 

176 

177 fminus = disp_minus.forces() 

178 dminus = disp_minus.dipole() 

179 

180 fplus = disp_plus.forces() 

181 dplus = disp_plus.dipole() 

182 

183 if self.nfree == 4: 

184 disp_mm = self._disp(a, i, -2) 

185 disp_pp = self._disp(a, i, 2) 

186 fminusminus = disp_mm.forces() 

187 dminusminus = disp_mm.dipole() 

188 

189 fplusplus = disp_pp.forces() 

190 dplusplus = disp_pp.dipole() 

191 if self.method == 'frederiksen': 

192 fminus[a] += -fminus.sum(0) 

193 fplus[a] += -fplus.sum(0) 

194 if self.nfree == 4: 

195 fminusminus[a] += -fminus.sum(0) 

196 fplusplus[a] += -fplus.sum(0) 

197 if self.nfree == 2: 

198 H[r] = (fminus - fplus)[self.indices].ravel() / 2.0 

199 dpdx[r] = (dminus - dplus) 

200 if self.nfree == 4: 

201 H[r] = (-fminusminus + 8 * fminus - 8 * fplus + 

202 fplusplus)[self.indices].ravel() / 12.0 

203 dpdx[r] = (-dplusplus + 8 * dplus - 8 * dminus + 

204 dminusminus) / 6.0 

205 H[r] /= 2 * self.delta 

206 dpdx[r] /= 2 * self.delta 

207 for n in range(3): 

208 if n not in self.directions: 

209 dpdx[r][n] = 0 

210 dpdx[r][n] = 0 

211 r += 1 

212 # Calculate eigenfrequencies and eigenvectors 

213 masses = self.atoms.get_masses() 

214 H += H.copy().T 

215 self.H = H 

216 

217 self.im = np.repeat(masses[self.indices]**-0.5, 3) 

218 omega2, modes = np.linalg.eigh(self.im[:, None] * H * self.im) 

219 self.modes = modes.T.copy() 

220 

221 # Calculate intensities 

222 dpdq = np.array([dpdx[j] / sqrt(masses[self.indices[j // 3]] * 

223 units._amu / units._me) 

224 for j in range(ndof)]) 

225 dpdQ = np.dot(dpdq.T, modes) 

226 dpdQ = dpdQ.T 

227 intensities = np.array([sum(dpdQ[j]**2) for j in range(ndof)]) 

228 # Conversion factor: 

229 s = units._hbar * 1e10 / sqrt(units._e * units._amu) 

230 self.hnu = s * omega2.astype(complex)**0.5 

231 # Conversion factor from atomic units to (D/Angstrom)^2/amu. 

232 conv = (1.0 / units.Debye)**2 * units._amu / units._me 

233 self.intensities = intensities * conv 

234 

235 def intensity_prefactor(self, intensity_unit): 

236 if intensity_unit == '(D/A)2/amu': 

237 return 1.0, '(D/Å)^2 amu^-1' 

238 elif intensity_unit == 'km/mol': 

239 # conversion factor from Porezag PRB 54 (1996) 7830 

240 return 42.255, 'km/mol' 

241 else: 

242 raise RuntimeError('Intensity unit >' + intensity_unit + 

243 '< unknown.') 

244 

245 def summary(self, method='standard', direction='central', 

246 intensity_unit='(D/A)2/amu', log=stdout): 

247 hnu = self.get_energies(method, direction) 

248 s = 0.01 * units._e / units._c / units._hplanck 

249 iu, iu_string = self.intensity_prefactor(intensity_unit) 

250 if intensity_unit == '(D/A)2/amu': 

251 iu_format = '%9.4f' 

252 elif intensity_unit == 'km/mol': 

253 iu_string = ' ' + iu_string 

254 iu_format = ' %7.1f' 

255 if isinstance(log, str): 

256 log = paropen(log, 'a') 

257 

258 parprint('-------------------------------------', file=log) 

259 parprint(' Mode Frequency Intensity', file=log) 

260 parprint(' # meV cm^-1 ' + iu_string, file=log) 

261 parprint('-------------------------------------', file=log) 

262 for n, e in enumerate(hnu): 

263 if e.imag != 0: 

264 c = 'i' 

265 e = e.imag 

266 else: 

267 c = ' ' 

268 e = e.real 

269 parprint(('%3d %6.1f%s %7.1f%s ' + iu_format) % 

270 (n, 1000 * e, c, s * e, c, iu * self.intensities[n]), 

271 file=log) 

272 parprint('-------------------------------------', file=log) 

273 parprint('Zero-point energy: %.3f eV' % self.get_zero_point_energy(), 

274 file=log) 

275 parprint('Static dipole moment: %.3f D' % self.dipole_zero, file=log) 

276 parprint('Maximum force on atom in `equilibrium`: %.4f eV/Å' % 

277 self.force_zero, file=log) 

278 parprint(file=log) 

279 

280 def get_spectrum(self, start=800, end=4000, npts=None, width=4, 

281 type='Gaussian', method='standard', direction='central', 

282 intensity_unit='(D/A)2/amu', normalize=False): 

283 """Get infrared spectrum. 

284 

285 The method returns wavenumbers in cm^-1 with corresponding 

286 absolute infrared intensity. 

287 Start and end point, and width of the Gaussian/Lorentzian should 

288 be given in cm^-1. 

289 normalize=True ensures the integral over the peaks to give the 

290 intensity. 

291 """ 

292 frequencies = self.get_frequencies(method, direction).real 

293 intensities = self.intensities 

294 return self.fold(frequencies, intensities, 

295 start, end, npts, width, type, normalize) 

296 

297 def write_spectra(self, out='ir-spectra.dat', start=800, end=4000, 

298 npts=None, width=10, 

299 type='Gaussian', method='standard', direction='central', 

300 intensity_unit='(D/A)2/amu', normalize=False): 

301 """Write out infrared spectrum to file. 

302 

303 First column is the wavenumber in cm^-1, the second column the 

304 absolute infrared intensities, and 

305 the third column the absorbance scaled so that data runs 

306 from 1 to 0. Start and end 

307 point, and width of the Gaussian/Lorentzian should be given 

308 in cm^-1.""" 

309 energies, spectrum = self.get_spectrum(start, end, npts, width, 

310 type, method, direction, 

311 normalize) 

312 

313 # Write out spectrum in file. First column is absolute intensities. 

314 # Second column is absorbance scaled so that data runs from 1 to 0 

315 spectrum2 = 1. - spectrum / spectrum.max() 

316 outdata = np.empty([len(energies), 3]) 

317 outdata.T[0] = energies 

318 outdata.T[1] = spectrum 

319 outdata.T[2] = spectrum2 

320 with open(out, 'w') as fd: 

321 fd.write('# %s folded, width=%g cm^-1\n' % (type.title(), width)) 

322 iu, iu_string = self.intensity_prefactor(intensity_unit) 

323 if normalize: 

324 iu_string = 'cm ' + iu_string 

325 fd.write('# [cm^-1] %14s\n' % ('[' + iu_string + ']')) 

326 for row in outdata: 

327 fd.write('%.3f %15.5e %15.5e \n' % 

328 (row[0], iu * row[1], row[2]))