index | modules | gitlab

# Source code for ase.vibrations.infrared

    # -*- coding: utf-8 -*-

    """Infrared intensities"""

    import os.path as op
    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
    from ase.utils import basestring, pickleload

    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.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
        -------------------------------------
        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.set_calculator(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
            self.ram = False

        def read(self, method='standard', direction='central'):
            self.method = method.lower()
            self.direction = direction.lower()
            assert self.method in ['standard', 'frederiksen']

            def load(fname, combined_data=None):
                if combined_data is not None:
                    try:
                        return combined_data[op.basename(fname)]
                    except KeyError:
                        return combined_data[fname]  # Old version
                return pickleload(open(fname, 'rb'))

            if direction != 'central':
                raise NotImplementedError(
                    'Only central difference is implemented at the moment.')

            if op.isfile(self.name + '.all.pckl'):
                # Open the combined pickle-file
                combined_data = load(self.name + '.all.pckl')
            else:
                combined_data = None
            # Get "static" dipole moment and forces
            name = '%s.eq.pckl' % self.name
            [forces_zero, dipole_zero] = load(name, combined_data)
            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 in self.indices:
                for i in 'xyz':
                    name = '%s.%d%s' % (self.name, a, i)
                    [fminus, dminus] = load(name + '-.pckl', combined_data)
                    [fplus, dplus] = load(name + '+.pckl', combined_data)
                    if self.nfree == 4:
                        [fminusminus, dminusminus] = load(
                            name + '--.pckl', combined_data)
                        [fplusplus, dplusplus] = load(
                            name + '++.pckl', combined_data)
                    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
            m = self.atoms.get_masses()
            H += H.copy().T
            self.H = H
            m = self.atoms.get_masses()
            self.im = np.repeat(m[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(m[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, basestring):
                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
            fd = open(out, 'w')
            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]))
            fd.close()
            # np.savetxt(out, outdata, fmt='%.3f  %15.5e  %15.5e')

    InfraRed = Infrared  # old name

Documentation v: 3.19.0

Versions
  3.20.0b1 (development)
  3.19.0 (latest stable)