Coverage for /builds/ase/ase/ase/vibrations/infrared.py: 50.39%
127 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""Infrared intensities"""
5from math import sqrt
6from sys import stdout
8import numpy as np
10import ase.units as units
11from ase.parallel import paropen, parprint
12from ase.vibrations.vibrations import Vibrations
15class Infrared(Vibrations):
16 """Class for calculating vibrational modes and infrared intensities
17 using finite difference.
19 The vibrational modes are calculated from a finite difference
20 approximation of the Dynamical matrix and the IR intensities from
21 a finite difference approximation of the gradient of the dipole
22 moment. The method is described in:
24 D. Porezag, M. R. Pederson:
25 "Infrared intensities and Raman-scattering activities within
26 density-functional theory",
27 Phys. Rev. B 54, 7830 (1996)
29 The calculator object (calc) linked to the Atoms object (atoms) must
30 have the attribute:
32 >>> calc.get_dipole_moment(atoms)
34 In addition to the methods included in the ``Vibrations`` class
35 the ``Infrared`` class introduces two new methods;
36 *get_spectrum()* and *write_spectra()*. The *summary()*, *get_energies()*,
37 *get_frequencies()*, *get_spectrum()* and *write_spectra()*
38 methods all take an optional *method* keyword. Use
39 method='Frederiksen' to use the method described in:
41 T. Frederiksen, M. Paulsson, M. Brandbyge, A. P. Jauho:
42 "Inelastic transport theory from first-principles: methodology
43 and applications for nanoscale devices",
44 Phys. Rev. B 75, 205413 (2007)
46 atoms: Atoms object
47 The atoms to work on.
48 indices: list of int
49 List of indices of atoms to vibrate. Default behavior is
50 to vibrate all atoms.
51 name: str
52 Name to use for files.
53 delta: float
54 Magnitude of displacements.
55 nfree: int
56 Number of displacements per degree of freedom, 2 or 4 are
57 supported. Default is 2 which will displace each atom +delta
58 and -delta in each cartesian direction.
59 directions: list of int
60 Cartesian coordinates to calculate the gradient
61 of the dipole moment in.
62 For example directions = 2 only dipole moment in the z-direction will
63 be considered, whereas for directions = [0, 1] only the dipole
64 moment in the xy-plane will be considered. Default behavior is to
65 use the dipole moment in all directions.
67 Example:
69 >>> from ase.io import read
70 >>> from ase.calculators.vasp import Vasp
71 >>> from ase.vibrations import Infrared
72 >>> water = read('water.traj') # read pre-relaxed structure of water
73 >>> calc = Vasp(prec='Accurate',
74 ... ediff=1E-8,
75 ... isym=0,
76 ... idipol=4, # calculate the total dipole moment
77 ... dipol=water.get_center_of_mass(scaled=True),
78 ... ldipol=True)
79 >>> water.calc = calc
80 >>> ir = Infrared(water)
81 >>> ir.run()
82 >>> ir.summary()
83 -------------------------------------
84 Mode Frequency Intensity
85 # meV cm^-1 (D/Å)^2 amu^-1
86 -------------------------------------
87 0 16.9i 136.2i 1.6108
88 1 10.5i 84.9i 2.1682
89 2 5.1i 41.1i 1.7327
90 3 0.3i 2.2i 0.0080
91 4 2.4 19.0 0.1186
92 5 15.3 123.5 1.4956
93 6 195.5 1576.7 1.6437
94 7 458.9 3701.3 0.0284
95 8 473.0 3814.6 1.1812
96 -------------------------------------
97 Zero-point energy: 0.573 eV
98 Static dipole moment: 1.833 D
99 Maximum force on atom in `equilibrium`: 0.0026 eV/Å
103 This interface now also works for calculator 'siesta',
104 (added get_dipole_moment for siesta).
106 Example:
108 >>> #!/usr/bin/env python3
110 >>> from ase.io import read
111 >>> from ase.calculators.siesta import Siesta
112 >>> from ase.vibrations import Infrared
114 >>> bud = read('bud1.xyz')
116 >>> calc = Siesta(label='bud',
117 ... meshcutoff=250 * Ry,
118 ... basis='DZP',
119 ... kpts=[1, 1, 1])
121 >>> calc.set_fdf('DM.MixingWeight', 0.08)
122 >>> calc.set_fdf('DM.NumberPulay', 3)
123 >>> calc.set_fdf('DM.NumberKick', 20)
124 >>> calc.set_fdf('DM.KickMixingWeight', 0.15)
125 >>> calc.set_fdf('SolutionMethod', 'Diagon')
126 >>> calc.set_fdf('MaxSCFIterations', 500)
127 >>> calc.set_fdf('PAO.BasisType', 'split')
128 >>> #50 meV = 0.003674931 * Ry
129 >>> calc.set_fdf('PAO.EnergyShift', 0.003674931 * Ry )
130 >>> calc.set_fdf('LatticeConstant', 1.000000 * Ang)
131 >>> calc.set_fdf('WriteCoorXmol', 'T')
133 >>> bud.calc = calc
135 >>> ir = Infrared(bud)
136 >>> ir.run()
137 >>> ir.summary()
139 """
141 def __init__(self, atoms, indices=None, name='ir', delta=0.01,
142 nfree=2, directions=None):
143 Vibrations.__init__(self, atoms, indices=indices, name=name,
144 delta=delta, nfree=nfree)
145 if atoms.constraints:
146 print('WARNING! \n Your Atoms object is constrained. '
147 'Some forces may be unintended set to zero. \n')
148 if directions is None:
149 self.directions = np.asarray([0, 1, 2])
150 else:
151 self.directions = np.asarray(directions)
152 self.ir = True
154 def read(self, method='standard', direction='central'):
155 self.method = method.lower()
156 self.direction = direction.lower()
157 assert self.method in ['standard', 'frederiksen']
159 if direction != 'central':
160 raise NotImplementedError(
161 'Only central difference is implemented at the moment.')
163 disp = self._eq_disp()
164 forces_zero = disp.forces()
165 dipole_zero = disp.dipole()
166 self.dipole_zero = (sum(dipole_zero**2)**0.5) / units.Debye
167 self.force_zero = max(
168 sum((forces_zero[j])**2)**0.5 for j in self.indices)
170 ndof = 3 * len(self.indices)
171 H = np.empty((ndof, ndof))
172 dpdx = np.empty((ndof, 3))
173 for r, (a, i) in enumerate(self._iter_ai()):
174 disp_minus = self._disp(a, i, -1)
175 disp_plus = self._disp(a, i, 1)
177 fminus = disp_minus.forces()
178 dminus = disp_minus.dipole()
180 fplus = disp_plus.forces()
181 dplus = disp_plus.dipole()
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()
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 # Calculate eigenfrequencies and eigenvectors
212 masses = self.atoms.get_masses()
213 H += H.copy().T
214 self.H = H
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()
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
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.')
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')
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)
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.
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)
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.
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)
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(f'# {type.title()} folded, width={width:g} cm^-1\n')
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]))