Chimie Théorique » QMX

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

"""Vibrational modes."""

import pickle

from math import sin, pi, sqrt

from os import remove

from os.path import isfile

import sys

import numpy as np

import ase.units as units

from ase.io.trajectory import PickleTrajectory

from ase.parallel import rank, barrier, paropen

class Vibrations:

    """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 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 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.optimizers 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  19:16:06        0.042171       2.9357

    BFGS:   1  19:16:07        0.104197       3.9270

    BFGS:   2  19:16:07        0.000963       0.4142

    BFGS:   3  19:16:07        0.000027       0.0698

    BFGS:   4  19:16:07        0.000000       0.0010

    >>> vib = Vibrations(n2)

    >>> vib.run()

    >>> vib.summary()

    ---------------------

      #    meV     cm^-1

    ---------------------

      0    0.0i      0.0i

      1    0.0i      0.0i

      2    0.0i      0.0i

      3    1.6      13.1 

      4    1.6      13.1 

      5  232.7    1877.2 

    ---------------------

    Zero-point energy: 0.118 eV

    Thermodynamic properties at 298.00 K

    Enthalpy: 0.050 eV

    Entropy : 0.648 meV/K

    T*S     : 0.193 eV

    E->G    : -0.025 eV

    >>> vib.write_mode(-1)  # write last mode to trajectory file

    """

    def __init__(self, atoms, indices=None, name='vib', delta=0.01, nfree=2):

        assert nfree in [2, 4]

        self.atoms = atoms

        if indices is None:

            indices = range(len(atoms))

        self.indices = np.asarray(indices)

        self.name = name

        self.delta = delta

        self.nfree = nfree

        self.H = None

        self.ir = None

    def run(self):

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

        filename = self.name + '.eq.pckl'

        if not isfile(filename):

            barrier()

            forces = self.atoms.get_forces()

            if self.ir:

                dipole = self.calc.get_dipole_moment(self.atoms)

            if rank == 0:

                fd = open(filename, 'w')

                if self.ir:

                    pickle.dump([forces, dipole], fd)

                    sys.stdout.write(

                        'Writing %s, dipole moment = (%.6f %.6f %.6f)\n' % 

                        (filename, dipole[0], dipole[1], dipole[2]))

                else:

                    pickle.dump(forces, fd)

                    sys.stdout.write('Writing %s\n' % filename)

                fd.close()

            sys.stdout.flush()

        p = self.atoms.positions.copy()

        for a in self.indices:

            for i in range(3):

                for sign in [-1, 1]:

                    for ndis in range(1, self.nfree//2+1):

                        filename = ('%s.%d%s%s.pckl' %

                                    (self.name, a, 'xyz'[i], ndis*' +-'[sign]))

                        if isfile(filename):

                            continue

                        barrier()

                        self.atoms.positions[a, i] = (p[a, i] +

                                                      ndis * sign * self.delta)

                        forces = self.atoms.get_forces()

                        if self.ir:

                            dipole = self.calc.get_dipole_moment(self.atoms)

                        if rank == 0:

                            fd = open(filename, 'w')

                            if self.ir:

                                pickle.dump([forces, dipole], fd)

                                sys.stdout.write(

                                    'Writing %s, ' % filename +

                                    'dipole moment = (%.6f %.6f %.6f)\n' % 

                                    (dipole[0], dipole[1], dipole[2]))

                            else:

                                pickle.dump(forces, fd)

                                sys.stdout.write('Writing %s\n' % filename)

                            fd.close()

                        sys.stdout.flush()

                        self.atoms.positions[a, i] = p[a, i]

        self.atoms.set_positions(p)

    def clean(self):

        if isfile(self.name + '.eq.pckl'):

            remove(self.name + '.eq.pckl')

        for a in self.indices:

            for i in 'xyz':

                for sign in '-+':

                    for ndis in range(1, self.nfree/2+1):

                        name = '%s.%d%s%s.pckl' % (self.name, a, i, ndis*sign)

                        if isfile(name):

                            remove(name)

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

        self.method = method.lower()

        self.direction = direction.lower()

        assert self.method in ['standard', 'frederiksen']

        assert self.direction in ['central', 'forward', 'backward']

        n = 3 * len(self.indices)

        H = np.empty((n, n))

        r = 0

        if direction != 'central':

            feq = pickle.load(open(self.name + '.eq.pckl'))

        for a in self.indices:

            for i in 'xyz':

                name = '%s.%d%s' % (self.name, a, i)

                fminus = pickle.load(open(name + '-.pckl'))

                fplus = pickle.load(open(name + '+.pckl'))

                if self.method == 'frederiksen':

                    fminus[a] -= fminus.sum(0)

                    fplus[a] -= fplus.sum(0)

                if self.nfree == 4:

                    fminusminus = pickle.load(open(name + '--.pckl'))

                    fplusplus = pickle.load(open(name + '++.pckl'))

                    if self.method == 'frederiksen':

                        fminusminus[a] -= fminusminus.sum(0)

                        fplusplus[a] -= fplusplus.sum(0)

                if self.direction == 'central':

                    if self.nfree == 2:

                        H[r] = .5 * (fminus - fplus)[self.indices].ravel()

                    else:

                        H[r] = H[r] = (-fminusminus +

                                       8 * fminus -

                                       8 * fplus +

                                       fplusplus)[self.indices].ravel() / 12.0

                elif self.direction == 'forward':

                    H[r] = (feq - fplus)[self.indices].ravel()

                else: # self.direction == 'backward':

                    H[r] = (fminus - feq)[self.indices].ravel()

                H[r] /= 2 * self.delta

                r += 1

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

        # Conversion factor:

        s = units._hbar * 1e10 / sqrt(units._e * units._amu)

        self.hnu = s * omega2.astype(complex)**0.5

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

        """Get vibration energies in eV."""

        if (self.H is None or method.lower() != self.method or

            direction.lower() != self.direction):

            self.read(method, direction)

        return self.hnu

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

        """Get vibration frequencies in cm^-1."""

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

        return s * self.get_energies(method, direction)

    def summary(self, method='standard', direction='central', T=298., 

                threshold=10, freq=None, log=sys.stdout):

        """Print a summary of the frequencies and derived thermodynamic

        properties. The equations for the calculation of the enthalpy and 

        entropy diverge for zero frequencies and a threshold value is used

        to ignore extremely low frequencies (default = 10 cm^-1).

        Parameters:

        method : string

            Can be 'standard'(default) or 'Frederiksen'.

        direction: string

            Direction for finite differences. Can be one of 'central'

            (default), 'forward', 'backward'.

        T : float

            Temperature in K at which thermodynamic properties are calculated.

        threshold : float

            Threshold value for low frequencies (default 10 cm^-1).

        freq : numpy array

            Optional. Can be used to create a summary on a set of known

            frequencies.

        destination : string or object with a .write() method

            Where to print the summary info. Default is to print to

            standard output. If another string is given, it creates that

            text file and writes the output to it. If a file object (or any

            object with a .write(string) function) is given, it writes to that

            file.

        Notes:

        The enthalpy and entropy calculations are very sensitive to low

        frequencies and the threshold value must be used with caution."""

        if isinstance(log, str):

            log = paropen(log, 'a')

        write = log.write

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

        if freq != None:

            hnu = freq / s

        else:

            hnu = self.get_energies(method, direction)

        write('---------------------\n')

        write('  #    meV     cm^-1\n')

        write('---------------------\n')

        for n, e in enumerate(hnu):

            if e.imag != 0:

                c = 'i'

                e = e.imag

            else:

                c = ' '

            write('%3d %6.1f%s  %7.1f%s\n' % (n, 1000 * e, c, s * e, c))

        write('---------------------\n')

        write('Zero-point energy: %.3f eV\n' %

              self.get_zero_point_energy(freq=freq))

        write('Thermodynamic properties at %.2f K\n' % T)

        write('Enthalpy: %.3f eV\n' % self.get_enthalpy(method=method,

                                                        direction=direction,

                                                        T=T,

                                                        threshold=threshold,

                                                        freq=freq))

        write('Entropy : %.3f meV/K\n' % 

              (1E3 * self.get_entropy(method=method,

                                      direction=direction,

                                      T=T,

                                      threshold=threshold,

                                      freq=freq)))

        write('T*S     : %.3f eV\n' %

              (T * self.get_entropy(method=method,

                                    direction=direction,

                                    T=T,

                                    threshold=threshold,

                                    freq=freq)))

        write('E->G    : %.3f eV\n' %

              (self.get_zero_point_energy(freq=freq) +

               self.get_enthalpy(method=method,

                                 direction=direction,

                                 T=T,

                                 threshold=threshold,

                                 freq=freq) -

               T * self.get_entropy(method=method,

                                    direction=direction,

                                    T=T,

                                    threshold=threshold,

                                    freq=freq)))

    def get_zero_point_energy(self, freq=None):

        if freq is None:

            return 0.5 * self.hnu.real.sum()

        else:

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

            return 0.5 * freq.real.sum() / s

    def get_enthalpy(self, method='standard', direction='central', T=298.0,

                     threshold=10, freq=None):

        H = 0.0

        if freq is None:

            freq = self.get_frequencies(method=method, direction=direction)

        for f in freq:

            if f.imag == 0 and f.real >= threshold:

                # The formula diverges for f->0

                x = (f.real * 100 * units._hplanck * units._c) / units._k

                H += units._k / units._e * x / (np.exp(x / T) - 1)

        return H

    def get_entropy(self, method='standard', direction='central', T=298.0,

                    threshold=10, freq=None):

        S = 0.0

        if freq is None:

            freq=self.get_frequencies(method=method, direction=direction)

        for f in freq:

            if f.imag == 0 and f.real >= threshold:

                # The formula diverges for f->0

                x = (f.real * 100 * units._hplanck * units._c) / units._k

                S += (units._k / units._e * (((x / T) / (np.exp(x / T) - 1)) -

                                             np.log(1 - np.exp(-x / T))))

        return S

    def get_mode(self, n):

        mode = np.zeros((len(self.atoms), 3))

        mode[self.indices] = (self.modes[n] * self.im).reshape((-1, 3))

        return mode

    def write_mode(self, n, kT=units.kB * 300, nimages=30):

        """Write mode to trajectory file."""

        mode = self.get_mode(n) * sqrt(kT / abs(self.hnu[n]))

        p = self.atoms.positions.copy()

        n %= 3 * len(self.indices)

        traj = PickleTrajectory('%s.%d.traj' % (self.name, n), 'w')

        calc = self.atoms.get_calculator()

        self.atoms.set_calculator()

        for x in np.linspace(0, 2 * pi, nimages, endpoint=False):

            self.atoms.set_positions(p + sin(x) * mode)

            traj.write(self.atoms)

        self.atoms.set_positions(p)

        self.atoms.set_calculator(calc)

        traj.close()

    def write_jmol(self):

        """Writes file for viewing of the modes with jmol."""

        fd = open(self.name + '.xyz', 'w')

        symbols = self.atoms.get_chemical_symbols()

        f = self.get_frequencies()

        for n in range(3 * len(self.atoms)):

            fd.write('%6d\n' % len(self.atoms))

            if f[n].imag != 0:

                c = 'i'

                f[n] = f[n].imag

            else:

                c = ' ' 

            fd.write('Mode #%d, f = %.1f%s cm^-1' % (n, f[n], c))

            if self.ir:

                fd.write(', I = %.4f (D/Å)^2 amu^-1.\n' % self.intensities[n])

            else:

                fd.write('.\n')

            mode = self.get_mode(n)

            for i, pos in enumerate(self.atoms.positions):

                fd.write('%2s %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f \n' % 

                         (symbols[i], pos[0], pos[1], pos[2],

                          mode[i,0], mode[i,1], mode[i,2]))

        fd.close()