Chimie Théorique » QMX

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

"""Infrared intensities"""

import pickle

from math import sin, pi, sqrt, exp, log

from os import remove

from os.path import isfile

import numpy as np

import ase.units as units

from ase.io.trajectory import PickleTrajectory

from ase.parallel import rank, barrier, parprint

from ase.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. Peterson:

      "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.infrared import InfraRed

    >>> water = read('water.traj')  # read pre-relaxed structure of water molecule

    >>> 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 `eqiulibrium`: 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.infrared 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):

        assert nfree in [2, 4]

        self.atoms = atoms

        if atoms.constraints:

            print "WARNING! \n Your Atoms object is constrained. Some forces may be unintended set to zero. \n"

        self.calc = atoms.get_calculator()

        if indices is None:

            indices = range(len(atoms))

        self.indices = np.asarray(indices)

        self.nfree = nfree

        self.name = name+'-d%.3f' % delta

        self.delta = delta

        self.H = None

        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.')

        # Get "static" dipole moment and forces

        name = '%s.eq.pckl' % self.name

        [forces_zero, dipole_zero] = pickle.load(open(name))

        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] = pickle.load(open(name + '-.pckl'))

                [fplus, dplus] = pickle.load(open(name + '+.pckl'))

                if self.nfree == 4:

                    [fminusminus, dminusminus] = pickle.load(open(name + '--.pckl'))

                    [fplusplus, dplusplus] = pickle.load(open(name + '++.pckl'))

                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 = units.Debye**2*units._amu/units._me

        self.intensities = intensities*conv

    def summary(self, method='standard', direction='central'):

        hnu = self.get_energies(method, direction)

        s = 0.01 * units._e / units._c / units._hplanck

        parprint('-------------------------------------')

        parprint(' Mode    Frequency        Intensity')

        parprint('  #    meV     cm^-1   (D/Å)^2 amu^-1')

        parprint('-------------------------------------')

        for n, e in enumerate(hnu):

            if e.imag != 0:

                c = 'i'

                e = e.imag

            else:

                c = ' '

            parprint('%3d %6.1f%s  %7.1f%s  %9.4f' %

                     (n, 1000 * e, c, s * e, c, self.intensities[n]))

        parprint('-------------------------------------')

        parprint('Zero-point energy: %.3f eV' % self.get_zero_point_energy())

        parprint('Static dipole moment: %.3f D' % self.dipole_zero)

        parprint('Maximum force on atom in `eqiulibrium`: %.4f eV/Å' %

                  self.force_zero)

        parprint()

    def get_spectrum(self, start=800, end=4000, npts=None, width=4, type='Gaussian', method='standard', direction='central'):

        """Get infrared spectrum.

        The method returns wavenumbers in cm^-1 with corresonding absolute infrared intensity.

        Start and end point, and width of the Gaussian/Lorentzian should be given in cm^-1."""

        self.type = type.lower()

        assert self.type in ['gaussian', 'lorentzian']

        if not npts:

            npts = (end-start)/width*10+1

        frequencies = self.get_frequencies(method, direction).real

        intensities=self.intensities

        if type == 'lorentzian':

            intensities = intensities*width*pi/2.

        else:

            sigma = width/2./sqrt(2.*log(2.))

        #Make array with spectrum data

        spectrum = np.empty(npts,np.float)

        energies = np.empty(npts,np.float)

        ediff = (end-start)/float(npts-1)

        energies = np.arange(start, end+ediff/2, ediff)

        for i, energy in enumerate(energies):

            energies[i] = energy

            if type == 'lorentzian':

                spectrum[i] = (intensities*0.5*width/pi/((frequencies-energy)**2+0.25*width**2)).sum()

            else:

                spectrum[i] = (intensities*np.exp(-(frequencies - energy)**2/2./sigma**2)).sum()

        return [energies, spectrum]

    def write_spectra(self, out='ir-spectra.dat', start=800, end=4000, npts=None, width=10, type='Gaussian', method='standard', direction='central'):

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

        #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')

        for row in outdata:

            fd.write('%.3f  %15.5e  %15.5e \n' % (row[0], row[1], row[2]) )

        fd.close()

        #np.savetxt(out, outdata, fmt='%.3f  %15.5e  %15.5e')