Coverage for /builds/kinetik161/ase/ase/vibrations/infrared.py: 51.16%
129 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Infrared intensities"""
3from math import sqrt
4from sys import stdout
6import numpy as np
8import ase.units as units
9from ase.parallel import paropen, parprint
10from ase.vibrations.vibrations import Vibrations
13class Infrared(Vibrations):
14 """Class for calculating vibrational modes and infrared intensities
15 using finite difference.
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:
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)
27 The calculator object (calc) linked to the Atoms object (atoms) must
28 have the attribute:
30 >>> calc.get_dipole_moment(atoms)
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:
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)
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.
65 Example:
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/Å
101 This interface now also works for calculator 'siesta',
102 (added get_dipole_moment for siesta).
104 Example:
106 >>> #!/usr/bin/env python3
108 >>> from ase.io import read
109 >>> from ase.calculators.siesta import Siesta
110 >>> from ase.vibrations import Infrared
112 >>> bud = read('bud1.xyz')
114 >>> calc = Siesta(label='bud',
115 ... meshcutoff=250 * Ry,
116 ... basis='DZP',
117 ... kpts=[1, 1, 1])
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')
131 >>> bud.calc = calc
133 >>> ir = Infrared(bud)
134 >>> ir.run()
135 >>> ir.summary()
137 """
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
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']
157 if direction != 'central':
158 raise NotImplementedError(
159 'Only central difference is implemented at the moment.')
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])
168 ndof = 3 * len(self.indices)
169 H = np.empty((ndof, ndof))
170 dpdx = np.empty((ndof, 3))
171 r = 0
173 for a, i in 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 r += 1
212 # Calculate eigenfrequencies and eigenvectors
213 masses = self.atoms.get_masses()
214 H += H.copy().T
215 self.H = H
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()
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
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.')
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')
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)
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.
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)
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.
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)
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(f'# {type.title()} folded, width={width:g} cm^-1\n')
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]))