Skip to main content
`Ctrl`+`K`

#

# Source code for ase.vibrations.infrared

    # fmt: off

    """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

    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)

        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)

        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(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]))