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 def __init__(self, atoms, indices=None, name='ir', delta=0.01, 

139 nfree=2, directions=None): 

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

141 delta=delta, nfree=nfree) 

142 if atoms.constraints: 

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

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

145 if directions is None: 

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

147 else: 

148 self.directions = np.asarray(directions) 

149 self.ir = True 

150 

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

152 self.method = method.lower() 

153 self.direction = direction.lower() 

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

155 

156 if direction != 'central': 

157 raise NotImplementedError( 

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

159 

160 disp = self._eq_disp() 

161 forces_zero = disp.forces() 

162 dipole_zero = disp.dipole() 

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

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

165 for j in self.indices]) 

166 

167 ndof = 3 * len(self.indices) 

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

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

170 r = 0 

171 

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

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

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

175 

176 fminus = disp_minus.forces() 

177 dminus = disp_minus.dipole() 

178 

179 fplus = disp_plus.forces() 

180 dplus = disp_plus.dipole() 

181 

182 if self.nfree == 4: 

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

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

185 fminusminus = disp_mm.forces() 

186 dminusminus = disp_mm.dipole() 

187 

188 fplusplus = disp_pp.forces() 

189 dplusplus = disp_pp.dipole() 

190 if self.method == 'frederiksen': 

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

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

193 if self.nfree == 4: 

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

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

196 if self.nfree == 2: 

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

198 dpdx[r] = (dminus - dplus) 

199 if self.nfree == 4: 

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

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

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

203 dminusminus) / 6.0 

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

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

206 for n in range(3): 

207 if n not in self.directions: 

208 dpdx[r][n] = 0 

209 dpdx[r][n] = 0 

210 r += 1 

211 # Calculate eigenfrequencies and eigenvectors 

212 masses = self.atoms.get_masses() 

213 H += H.copy().T 

214 self.H = H 

215 

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

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

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

219 

220 # Calculate intensities 

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

222 units._amu / units._me) 

223 for j in range(ndof)]) 

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

225 dpdQ = dpdQ.T 

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

227 # Conversion factor: 

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

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

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

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

232 self.intensities = intensities * conv 

233 

234 def intensity_prefactor(self, intensity_unit): 

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

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

237 elif intensity_unit == 'km/mol': 

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

239 return 42.255, 'km/mol' 

240 else: 

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

242 '< unknown.') 

243 

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

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

246 hnu = self.get_energies(method, direction) 

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

248 iu, iu_string = self.intensity_prefactor(intensity_unit) 

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

250 iu_format = '%9.4f' 

251 elif intensity_unit == 'km/mol': 

252 iu_string = ' ' + iu_string 

253 iu_format = ' %7.1f' 

254 if isinstance(log, str): 

255 log = paropen(log, 'a') 

256 

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

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

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

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

261 for n, e in enumerate(hnu): 

262 if e.imag != 0: 

263 c = 'i' 

264 e = e.imag 

265 else: 

266 c = ' ' 

267 e = e.real 

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

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

270 file=log) 

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

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

273 file=log) 

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

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

276 self.force_zero, file=log) 

277 parprint(file=log) 

278 

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

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

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

282 """Get infrared spectrum. 

283 

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

285 absolute infrared intensity. 

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

287 be given in cm^-1. 

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

289 intensity. 

290 """ 

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

292 intensities = self.intensities 

293 return self.fold(frequencies, intensities, 

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

295 

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

297 npts=None, width=10, 

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

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

300 """Write out infrared spectrum to file. 

301 

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

303 absolute infrared intensities, and 

304 the third column the absorbance scaled so that data runs 

305 from 1 to 0. Start and end 

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

307 in cm^-1.""" 

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

309 type, method, direction, 

310 normalize) 

311 

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

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

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

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

316 outdata.T[0] = energies 

317 outdata.T[1] = spectrum 

318 outdata.T[2] = spectrum2 

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

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

321 iu, iu_string = self.intensity_prefactor(intensity_unit) 

322 if normalize: 

323 iu_string = 'cm ' + iu_string 

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

325 for row in outdata: 

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

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