"""Infrared intensities"""
from math import sqrt
from sys import stdout
import numpy as np
import ase.units as units
from ase.parallel import parprint, paropen
from ase.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))
r = 0
for a, i in 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
r += 1
# 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. - 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('# %s folded, width=%g cm^-1\n' % (type.title(), width))
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]))