Vibration analysis¶
Vibrational modes¶
You can calculate the vibrational modes of an
Atoms
object in the harmonic approximation using
the Vibrations
.

class
ase.vibrations.
Vibrations
(atoms, indices=None, name=’vib’, delta=0.01, nfree=2)[source]¶ Class for calculating vibrational modes using finite difference.
The vibrational modes are calculated from a finite difference approximation of the Hessian matrix.
The summary(), get_energies() and get_frequencies() 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 firstprinciples: 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 atom and cartesian coordinate, 2 and 4 are supported. Default is 2 which will displace each atom +delta and delta for each cartesian coordinate.
Example:
>>> from ase import Atoms >>> from ase.calculators.emt import EMT >>> from ase.optimize import BFGS >>> from ase.vibrations import Vibrations >>> n2 = Atoms('N2', [(0, 0, 0), (0, 0, 1.1)], ... calculator=EMT()) >>> BFGS(n2).run(fmax=0.01) BFGS: 0 16:01:21 0.440339 3.2518 BFGS: 1 16:01:21 0.271928 0.8211 BFGS: 2 16:01:21 0.263278 0.1994 BFGS: 3 16:01:21 0.262777 0.0088 >>> vib = Vibrations(n2) >>> vib.run() Writing vib.eq.pckl Writing vib.0x.pckl Writing vib.0x+.pckl Writing vib.0y.pckl Writing vib.0y+.pckl Writing vib.0z.pckl Writing vib.0z+.pckl Writing vib.1x.pckl Writing vib.1x+.pckl Writing vib.1y.pckl Writing vib.1y+.pckl Writing vib.1z.pckl Writing vib.1z+.pckl >>> vib.summary()  # meV cm^1  0 0.0 0.0 1 0.0 0.0 2 0.0 0.0 3 2.5 20.4 4 2.5 20.4 5 152.6 1230.8  Zeropoint energy: 0.079 eV >>> vib.write_mode(1) # write last mode to trajectory file

clean
(empty_files=False)[source]¶ Remove picklefiles.
Use empty_files=True to remove only empty files.

fold
(frequencies, intensities, start=800.0, end=4000.0, npts=None, width=4.0, type=’Gaussian’, normalize=False)[source]¶ Fold frequencies and intensities within the given range and folding method (Gaussian/Lorentzian). The energy unit is cm^1. normalize=True ensures the integral over the peaks to give the intensity.

get_frequencies
(method=’standard’, direction=’central’)[source]¶ Get vibration frequencies in cm^1.

run
()[source]¶ Run the vibration calculations.
This will calculate the forces for 6 displacements per atom +/x, +/y, +/z. Only those calculations that are not already done will be started. Be aware that an interrupted calculation may produce an empty file (ending with .pckl), which must be deleted before restarting the job. Otherwise the forces will not be calculated for that displacement.
Note that the calculations for the different displacements can be done simultaneously by several independent processes. This feature relies on the existence of files and the subsequent creation of the file in case it is not found.

summary
(method=’standard’, direction=’central’, freq=None, log=<_io.StringIO object>)[source]¶ Print a summary of the vibrational frequencies.
Parameters:
 method : string
 Can be ‘standard’(default) or ‘Frederiksen’.
 direction: string
 Direction for finite differences. Can be one of ‘central’ (default), ‘forward’, ‘backward’.
 freq : numpy array
 Optional. Can be used to create a summary on a set of known frequencies.
 log : if specified, write output to a different location than
 stdout. Can be an object with a write() method or the name of a file to create.

write_dos
(out=’vibdos.dat’, start=800, end=4000, npts=None, width=10, type=’Gaussian’, method=’standard’, direction=’central’)[source]¶ Write out the vibrational density of states to file.
First column is the wavenumber in cm^1, the second column the folded vibrational density of states. Start and end points, and width of the Gaussian/Lorentzian should be given in cm^1.
name is a string that is prefixed to the names of all the files created. atoms is an Atoms object that is either at a fully relaxed ground state or at a saddle point. freeatoms is a list of atom indices for which the vibrational modes will be calculated, the rest of the atoms are considered frozen. displacements is a list of displacements, one for each free atom that are used in the finite difference method to calculate the Hessian matrix. method is 1 for backward differences, 0 for centered differences, and 1 for forward differences.
Infrared intensities¶
Infrared
is an extension of
Vibrations
, in addition to the
vibrational modes, also the infrared intensities of the modes
are calculated for an Atoms
object.

class
ase.vibrations.
Infrared
(atoms, indices=None, name=’ir’, delta=0.01, nfree=2, directions=None)[source]¶ 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 Ramanscattering activities within densityfunctional 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 theInfrared
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 firstprinciples: 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 zdirection will be considered, whereas for directions = [0, 1] only the dipole moment in the xyplane 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 prerelaxed structure of water >>> calc = Vasp(prec='Accurate', ... ediff=1E8, ... isym=0, ... idipol=4, # calculate the total dipole moment ... dipol=water.get_center_of_mass(scaled=True), ... ldipol=True) >>> water.set_calculator(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  Zeropoint 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 python
>>> 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.set_calculator(calc)
>>> ir = Infrared(bud) >>> ir.run() >>> ir.summary()

get_spectrum
(start=800, end=4000, npts=None, width=4, type=’Gaussian’, method=’standard’, direction=’central’, intensity_unit=’(D/A)2/amu’, normalize=False)[source]¶ 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.

write_spectra
(out=’irspectra.dat’, start=800, end=4000, npts=None, width=10, type=’Gaussian’, method=’standard’, direction=’central’, intensity_unit=’(D/A)2/amu’, normalize=False)[source]¶ 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.