"""Infrared intensities"""
from math import sqrt
from sys import stdout
import numpy as np
import ase.units as units
from ase.parallel import paropen, parprint
from ase.vibrations.vibrations import Vibrations
[docs]
class Infrared(Vibrations):
"""Class for calculating vibrational modes and infrared intensities
using finite difference.
The vibrational modes are calculated from a finite difference
approximation of the Dynamical matrix and the IR intensities from
a finite difference approximation of the gradient of the dipole
moment. The method is described in:
D. Porezag, M. R. Pederson:
"Infrared intensities and Raman-scattering activities within
density-functional theory",
Phys. Rev. B 54, 7830 (1996)
The calculator object (calc) linked to the Atoms object (atoms) must
have the attribute:
>>> calc.get_dipole_moment(atoms)
In addition to the methods included in the ``Vibrations`` class
the ``Infrared`` class introduces two new methods;
*get_spectrum()* and *write_spectra()*. The *summary()*, *get_energies()*,
*get_frequencies()*, *get_spectrum()* and *write_spectra()*
methods all take an optional *method* keyword. Use
method='Frederiksen' to use the method described in:
T. Frederiksen, M. Paulsson, M. Brandbyge, A. P. Jauho:
"Inelastic transport theory from first-principles: methodology
and applications for nanoscale devices",
Phys. Rev. B 75, 205413 (2007)
atoms: Atoms object
The atoms to work on.
indices: list of int
List of indices of atoms to vibrate. Default behavior is
to vibrate all atoms.
name: str
Name to use for files.
delta: float
Magnitude of displacements.
nfree: int
Number of displacements per degree of freedom, 2 or 4 are
supported. Default is 2 which will displace each atom +delta
and -delta in each cartesian direction.
directions: list of int
Cartesian coordinates to calculate the gradient
of the dipole moment in.
For example directions = 2 only dipole moment in the z-direction will
be considered, whereas for directions = [0, 1] only the dipole
moment in the xy-plane will be considered. Default behavior is to
use the dipole moment in all directions.
Example:
>>> from ase.io import read
>>> from ase.calculators.vasp import Vasp
>>> from ase.vibrations import Infrared
>>> water = read('water.traj') # read pre-relaxed structure of water
>>> calc = Vasp(prec='Accurate',
... ediff=1E-8,
... isym=0,
... idipol=4, # calculate the total dipole moment
... dipol=water.get_center_of_mass(scaled=True),
... ldipol=True)
>>> water.calc = calc
>>> ir = Infrared(water)
>>> ir.run()
>>> ir.summary()
-------------------------------------
Mode Frequency Intensity
# meV cm^-1 (D/Å)^2 amu^-1
-------------------------------------
0 16.9i 136.2i 1.6108
1 10.5i 84.9i 2.1682
2 5.1i 41.1i 1.7327
3 0.3i 2.2i 0.0080
4 2.4 19.0 0.1186
5 15.3 123.5 1.4956
6 195.5 1576.7 1.6437
7 458.9 3701.3 0.0284
8 473.0 3814.6 1.1812
-------------------------------------
Zero-point energy: 0.573 eV
Static dipole moment: 1.833 D
Maximum force on atom in `equilibrium`: 0.0026 eV/Å
This interface now also works for calculator 'siesta',
(added get_dipole_moment for siesta).
Example:
>>> #!/usr/bin/env python3
>>> from ase.io import read
>>> from ase.calculators.siesta import Siesta
>>> from ase.vibrations import Infrared
>>> bud = read('bud1.xyz')
>>> calc = Siesta(label='bud',
... meshcutoff=250 * Ry,
... basis='DZP',
... kpts=[1, 1, 1])
>>> calc.set_fdf('DM.MixingWeight', 0.08)
>>> calc.set_fdf('DM.NumberPulay', 3)
>>> calc.set_fdf('DM.NumberKick', 20)
>>> calc.set_fdf('DM.KickMixingWeight', 0.15)
>>> calc.set_fdf('SolutionMethod', 'Diagon')
>>> calc.set_fdf('MaxSCFIterations', 500)
>>> calc.set_fdf('PAO.BasisType', 'split')
>>> #50 meV = 0.003674931 * Ry
>>> calc.set_fdf('PAO.EnergyShift', 0.003674931 * Ry )
>>> calc.set_fdf('LatticeConstant', 1.000000 * Ang)
>>> calc.set_fdf('WriteCoorXmol', 'T')
>>> bud.calc = calc
>>> ir = Infrared(bud)
>>> ir.run()
>>> ir.summary()
"""
def __init__(
self,
atoms,
indices=None,
name='ir',
delta=0.01,
nfree=2,
directions=None,
):
Vibrations.__init__(
self, atoms, indices=indices, name=name, delta=delta, nfree=nfree
)
if atoms.constraints:
print(
'WARNING! \n Your Atoms object is constrained. '
'Some forces may be unintended set to zero. \n'
)
if directions is None:
self.directions = np.asarray([0, 1, 2])
else:
self.directions = np.asarray(directions)
self.ir = True
def read(self, method='standard', direction='central'):
self.method = method.lower()
self.direction = direction.lower()
assert self.method in ['standard', 'frederiksen']
if direction != 'central':
raise NotImplementedError(
'Only central difference is implemented at the moment.'
)
disp = self._eq_disp()
forces_zero = disp.forces()
dipole_zero = disp.dipole()
self.dipole_zero = (sum(dipole_zero**2) ** 0.5) / units.Debye
self.force_zero = max(
sum((forces_zero[j]) ** 2) ** 0.5 for j in self.indices
)
ndof = 3 * len(self.indices)
H = np.empty((ndof, ndof))
dpdx = np.empty((ndof, 3))
for r, (a, i) in enumerate(self._iter_ai()):
disp_minus = self._disp(a, i, -1)
disp_plus = self._disp(a, i, 1)
fminus = disp_minus.forces()
dminus = disp_minus.dipole()
fplus = disp_plus.forces()
dplus = disp_plus.dipole()
if self.nfree == 4:
disp_mm = self._disp(a, i, -2)
disp_pp = self._disp(a, i, 2)
fminusminus = disp_mm.forces()
dminusminus = disp_mm.dipole()
fplusplus = disp_pp.forces()
dplusplus = disp_pp.dipole()
if self.method == 'frederiksen':
fminus[a] += -fminus.sum(0)
fplus[a] += -fplus.sum(0)
if self.nfree == 4:
fminusminus[a] += -fminus.sum(0)
fplusplus[a] += -fplus.sum(0)
if self.nfree == 2:
H[r] = (fminus - fplus)[self.indices].ravel() / 2.0
dpdx[r] = dminus - dplus
if self.nfree == 4:
H[r] = (-fminusminus + 8 * fminus - 8 * fplus + fplusplus)[
self.indices
].ravel() / 12.0
dpdx[r] = (
-dplusplus + 8 * dplus - 8 * dminus + dminusminus
) / 6.0
H[r] /= 2 * self.delta
dpdx[r] /= 2 * self.delta
for n in range(3):
if n not in self.directions:
dpdx[r][n] = 0
dpdx[r][n] = 0
# Calculate eigenfrequencies and eigenvectors
masses = self.atoms.get_masses()
H += H.copy().T
self.H = H
self.im = np.repeat(masses[self.indices] ** -0.5, 3)
omega2, modes = np.linalg.eigh(self.im[:, None] * H * self.im)
self.modes = modes.T.copy()
# Calculate intensities
dpdq = np.array(
[
dpdx[j]
/ sqrt(masses[self.indices[j // 3]] * units._amu / units._me)
for j in range(ndof)
]
)
dpdQ = np.dot(dpdq.T, modes)
dpdQ = dpdQ.T
intensities = np.array([sum(dpdQ[j] ** 2) for j in range(ndof)])
# Conversion factor:
s = units._hbar * 1e10 / sqrt(units._e * units._amu)
self.hnu = s * omega2.astype(complex) ** 0.5
# Conversion factor from atomic units to (D/Angstrom)^2/amu.
conv = (1.0 / units.Debye) ** 2 * units._amu / units._me
self.intensities = intensities * conv
def intensity_prefactor(self, intensity_unit):
if intensity_unit == '(D/A)2/amu':
return 1.0, '(D/Å)^2 amu^-1'
elif intensity_unit == 'km/mol':
# conversion factor from Porezag PRB 54 (1996) 7830
return 42.255, 'km/mol'
else:
raise RuntimeError(
'Intensity unit >' + intensity_unit + '< unknown.'
)
def summary(
self,
method='standard',
direction='central',
intensity_unit='(D/A)2/amu',
log=stdout,
):
hnu = self.get_energies(method, direction)
s = 0.01 * units._e / units._c / units._hplanck
iu, iu_string = self.intensity_prefactor(intensity_unit)
if intensity_unit == '(D/A)2/amu':
iu_format = '%9.4f'
elif intensity_unit == 'km/mol':
iu_string = ' ' + iu_string
iu_format = ' %7.1f'
if isinstance(log, str):
log = paropen(log, 'a')
parprint('-------------------------------------', file=log)
parprint(' Mode Frequency Intensity', file=log)
parprint(' # meV cm^-1 ' + iu_string, file=log)
parprint('-------------------------------------', file=log)
for n, e in enumerate(hnu):
if e.imag != 0:
c = 'i'
e = e.imag
else:
c = ' '
e = e.real
parprint(
('%3d %6.1f%s %7.1f%s ' + iu_format)
% (n, 1000 * e, c, s * e, c, iu * self.intensities[n]),
file=log,
)
parprint('-------------------------------------', file=log)
parprint(
'Zero-point energy: %.3f eV' % self.get_zero_point_energy(),
file=log,
)
parprint('Static dipole moment: %.3f D' % self.dipole_zero, file=log)
parprint(
'Maximum force on atom in `equilibrium`: %.4f eV/Å'
% self.force_zero,
file=log,
)
parprint(file=log)
[docs]
def get_spectrum(
self,
start=800,
end=4000,
npts=None,
width=4,
type='Gaussian',
method='standard',
direction='central',
intensity_unit='(D/A)2/amu',
normalize=False,
):
"""Get infrared spectrum.
The method returns wavenumbers in cm^-1 with corresponding
absolute infrared intensity.
Start and end point, and width of the Gaussian/Lorentzian should
be given in cm^-1.
normalize=True ensures the integral over the peaks to give the
intensity.
"""
frequencies = self.get_frequencies(method, direction).real
intensities = self.intensities
return self.fold(
frequencies, intensities, start, end, npts, width, type, normalize
)
[docs]
def write_spectra(
self,
out='ir-spectra.dat',
start=800,
end=4000,
npts=None,
width=10,
type='Gaussian',
method='standard',
direction='central',
intensity_unit='(D/A)2/amu',
normalize=False,
):
"""Write out infrared spectrum to file.
First column is the wavenumber in cm^-1, the second column the
absolute infrared intensities, and
the third column the absorbance scaled so that data runs
from 1 to 0. Start and end
point, and width of the Gaussian/Lorentzian should be given
in cm^-1."""
energies, spectrum = self.get_spectrum(
start, end, npts, width, type, method, direction, normalize
)
# Write out spectrum in file. First column is absolute intensities.
# Second column is absorbance scaled so that data runs from 1 to 0
spectrum2 = 1.0 - spectrum / spectrum.max()
outdata = np.empty([len(energies), 3])
outdata.T[0] = energies
outdata.T[1] = spectrum
outdata.T[2] = spectrum2
with open(out, 'w') as fd:
fd.write(f'# {type.title()} folded, width={width:g} cm^-1\n')
iu, iu_string = self.intensity_prefactor(intensity_unit)
if normalize:
iu_string = 'cm ' + iu_string
fd.write('# [cm^-1] %14s\n' % ('[' + iu_string + ']'))
for row in outdata:
fd.write(
'%.3f %15.5e %15.5e \n' % (row[0], iu * row[1], row[2])
)
