Chimie Théorique » QMX

"""Effective medium theory potential."""

from math import sqrt, exp, log, pi

import numpy as np

import sys

from ase.data import atomic_numbers, chemical_symbols

from ase.units import Bohr

parameters = {

#          E0     s0    V0     eta2    kappa   lambda  n0

#          eV     bohr  eV     bohr^-1 bohr^-1 bohr^-1 bohr^-3

    'H':  (-2.21, 0.71, 2.132, 1.652,  2.790,  1.892,  0.00547, 'dimer'),

    'Al': (-3.28, 3.00, 1.493, 1.240,  2.000,  1.169,  0.00700, 'fcc'),

    'Cu': (-3.51, 2.67, 2.476, 1.652,  2.740,  1.906,  0.00910, 'fcc'),

    'Ag': (-2.96, 3.01, 2.132, 1.652,  2.790,  1.892,  0.00547, 'fcc'),

    'Au': (-3.80, 3.00, 2.321, 1.674,  2.873,  2.182,  0.00703, 'fcc'),

    'Ni': (-4.44, 2.60, 3.673, 1.669,  2.757,  1.948,  0.01030, 'fcc'),

    'Pd': (-3.90, 2.87, 2.773, 1.818,  3.107,  2.155,  0.00688, 'fcc'),

    'Pt': (-5.85, 2.90, 4.067, 1.812,  3.145,  2.192,  0.00802, 'fcc'),

    'C':  (-1.97, 1.18, 0.132, 3.652,  5.790,  2.892,  0.01322, 'dimer'),

    'N':  (-4.97, 1.18, 0.132, 2.652,  3.790,  2.892,  0.01222, 'dimer'),

    'O':  (-2.97, 1.25, 2.132, 3.652,  5.790,  4.892,  0.00850, 'dimer')}

beta = 1.809     #(16 * pi / 3)**(1.0 / 3) / 2**0.5

#                 but preserve historical rounding.

#eta1 = 0.5 / Bohr

class EMT:

    disabled = False  # Set to True to disable (asap does this).

    def __init__(self, fix_alloys=False):

        self.energy = None

        self.fix_alloys = fix_alloys

        if self.disabled:

            print >>sys.stderr, """

            ase.EMT has been disabled by Asap.  Most likely, you

            intended to use Asap's EMT calculator, but accidentally

            imported ase's EMT calculator after Asap's.  This could

            happen if your script contains the lines

              from asap3 import *

              from ase.calculators.emt import EMT

            Swap the two lines to solve the problem.

            (or 'from ase import *' in older versions of ASE.) Swap

            the two lines to solve the problem.

            In the UNLIKELY event that you actually wanted to use

            ase.calculators.emt.EMT although asap3 is loaded into memory, please

            reactivate it with the command

              ase.calculators.emt.EMT.disabled = False

            """

            raise RuntimeError('ase.EMT has been disabled.  ' +

                               'See message printed above.')

    def get_spin_polarized(self):

        return False

    def initialize(self, atoms):

        self.par = {}

        self.rc = 0.0

        self.numbers = atoms.get_atomic_numbers()

        maxseq = 0.0

        seen = {}

        for Z in self.numbers:

            if Z not in seen:

                seen[Z] = True

                ss = parameters[chemical_symbols[Z]][1] * Bohr

                if maxseq < ss:

                    maxseq = ss

        rc = self.rc = beta * maxseq * 0.5 * (sqrt(3) + sqrt(4))

        rr = rc * 2 * sqrt(4) / (sqrt(3) + sqrt(4))

        self.acut = np.log(9999.0) / (rr - rc)

        for Z in self.numbers:

            if Z not in self.par:

                p = parameters[chemical_symbols[Z]]

                s0 = p[1] * Bohr

                eta2 = p[3] / Bohr

                kappa = p[4] / Bohr

                x = eta2 * beta * s0

                gamma1 = 0.0

                gamma2 = 0.0

                if p[7] == 'fcc':

                    for i, n in enumerate([12, 6, 24, 12]):

                        r = s0 * beta * sqrt(i + 1)

                        x = n / (12 * (1.0 + exp(self.acut * (r - rc))))

                        gamma1 += x * exp(-eta2 * (r - beta * s0))

                        gamma2 += x * exp(-kappa / beta * (r - beta * s0))

                elif p[7] == 'dimer':

                    r = s0 * beta

                    n = 1

                    x = n / (12 * (1.0 + exp(self.acut * (r - rc))))

                    gamma1 += x * exp(-eta2 * (r - beta * s0))

                    gamma2 += x * exp(-kappa / beta * (r - beta * s0))

                else:

                    raise RuntimeError

                self.par[Z] = {'E0': p[0],

                               's0': s0,

                               'V0': p[2],

                               'eta2': eta2,

                               'kappa': kappa,

                               'lambda': p[5] / Bohr,

                               'n0': p[6] / Bohr**3,

                               'rc': rc,

                               'gamma1': gamma1,

                               'gamma2': gamma2}

                #if rc + 0.5 > self.rc:

                #    self.rc = rc + 0.5

        self.ksi = {}

        for s1, p1 in self.par.items():

            self.ksi[s1] = {}

            for s2, p2 in self.par.items():

                #self.ksi[s1][s2] = (p2['n0'] / p1['n0'] *

                #                    exp(eta1 * (p1['s0'] - p2['s0'])))

                self.ksi[s1][s2] = p2['n0'] / p1['n0']

        self.forces = np.empty((len(atoms), 3))

        self.sigma1 = np.empty(len(atoms))

        self.deds = np.empty(len(atoms))

    def update(self, atoms):

        if (self.energy is None or

            len(self.numbers) != len(atoms) or

            (self.numbers != atoms.get_atomic_numbers()).any()):

            self.initialize(atoms)

            self.calculate(atoms)

        elif ((self.positions != atoms.get_positions()).any() or

              (self.pbc != atoms.get_pbc()).any() or

              (self.cell != atoms.get_cell()).any()):

            self.calculate(atoms)

    def calculation_required(self, atoms, quantities):

        if len(quantities) == 0:

            return False

        return (self.energy is None or

                len(self.numbers) != len(atoms) or

                (self.numbers != atoms.get_atomic_numbers()).any() or

                (self.positions != atoms.get_positions()).any() or

                (self.pbc != atoms.get_pbc()).any() or

                (self.cell != atoms.get_cell()).any())

    def get_potential_energy(self, atoms):

        self.update(atoms)

        return self.energy

    def get_numeric_forces(self, atoms):

        self.update(atoms)

        p = atoms.positions

        p0 = p.copy()

        forces = np.empty_like(p)

        eps = 0.0001

        for a in range(len(p)):

            for c in range(3):

                p[a, c] += eps

                self.calculate(atoms)

                de = self.energy

                p[a, c] -= 2 * eps

                self.calculate(atoms)

                de -= self.energy

                p[a, c] += eps

                forces[a, c] = -de / (2 * eps)

        p[:] = p0

        return forces

    def get_forces(self, atoms):

        self.update(atoms)

        return self.forces.copy()

    def get_stress(self, atoms):

        raise NotImplementedError

    def calculate(self, atoms):

        self.positions = atoms.get_positions().copy()

        self.cell = atoms.get_cell().copy()

        self.pbc = atoms.get_pbc().copy()

        icell = np.linalg.inv(self.cell)

        scaled = np.dot(self.positions, icell)

        N = []

        for i in range(3):

            if self.pbc[i]:

                scaled[:, i] %= 1.0

                v = icell[:, i]

                h = 1 / sqrt(np.dot(v, v))

                N.append(int(self.rc / h) + 1)

            else:

                N.append(0)

        R = np.dot(scaled, self.cell)

        self.energy = 0.0

        self.sigma1[:] = 0.0

        self.forces[:] = 0.0

        N1, N2, N3 = N

        natoms = len(atoms)

        for i1 in range(-N1, N1 + 1):

            for i2 in range(-N2, N2 + 1):

                for i3 in range(-N3, N3 + 1):

                    C = np.dot((i1, i2, i3), self.cell)

                    Q = R + C

                    c = (i1 == 0 and i2 == 0 and i3 == 0)

                    for a1 in range(natoms):

                        Z1 = self.numbers[a1]

                        p1 = self.par[Z1]

                        ksi = self.ksi[Z1]

                        for a2 in range(natoms):

                            if c and a2 == a1:

                                continue

                            d = Q[a2] - R[a1]

                            r = sqrt(np.dot(d, d))

                            if r < self.rc + 0.5:

                                Z2 = self.numbers[a2]

                                if self.fix_alloys:

                                    p2 = self.par[Z2]

                                else:

                                    p2 = p1  # Old buggy behaviour

                                self.interact1(a1, a2, d, r, p1, p2, ksi[Z2])

        for a in range(natoms):

            Z = self.numbers[a]

            p = self.par[Z]

            try:

                ds = -log(self.sigma1[a] / 12) / (beta * p['eta2'])

            except (OverflowError, ValueError):

                self.deds[a] = 0.0

                self.energy -= p['E0']

                continue

            x = p['lambda'] * ds

            y = exp(-x)

            z = 6 * p['V0'] * exp(-p['kappa'] * ds)

            self.deds[a] = ((x * y * p['E0'] * p['lambda'] + p['kappa'] * z) /

                            (self.sigma1[a] * beta * p['eta2']))

            e = p['E0'] * ((1 + x) * y - 1) + z

            self.energy += p['E0'] * ((1 + x) * y - 1) + z

        for i1 in range(-N1, N1 + 1):

            for i2 in range(-N2, N2 + 1):

                for i3 in range(-N3, N3 + 1):

                    C = np.dot((i1, i2, i3), self.cell)

                    Q = R + C

                    c = (i1 == 0 and i2 == 0 and i3 == 0)

                    for a1 in range(natoms):

                        Z1 = self.numbers[a1]

                        p1 = self.par[Z1]

                        ksi = self.ksi[Z1]

                        for a2 in range(natoms):

                            if c and a2 == a1:

                                continue

                            d = Q[a2] - R[a1]

                            r = sqrt(np.dot(d, d))

                            if r < self.rc + 0.5:

                                Z2 = self.numbers[a2]

                                if self.fix_alloys:

                                    p2 = self.par[Z2]

                                else:

                                    p2 = p1  # Old buggy behaviour

                                self.interact2(a1, a2, d, r, p1, p2, ksi[Z2])

    def interact1(self, a1, a2, d, r, p1, p2, ksi):

        x = exp(self.acut * (r - self.rc))

        theta = 1.0 / (1.0 + x)

        y = (0.5 * p1['V0'] * exp(-p2['kappa'] * (r / beta - p2['s0'])) *

             ksi / p1['gamma2'] * theta)

        self.energy -= y

        f = y * (p2['kappa'] / beta + self.acut * theta * x) * d / r

        self.forces[a1] += f

        self.forces[a2] -= f

        self.sigma1[a1] += (exp(-p2['eta2'] * (r - beta * p2['s0'])) *

                            ksi * theta / p1['gamma1'])

    def interact2(self, a1, a2, d, r, p1, p2, ksi):

        x = exp(self.acut * (r - self.rc))

        theta = 1.0 / (1.0 + x)

        y = (exp(-p2['eta2'] * (r - beta * p2['s0'])) *

             ksi / p1['gamma1'] * theta * self.deds[a1])

        f = y * (p2['eta2'] + self.acut * theta * x) * d / r

        self.forces[a1] -= f

        self.forces[a2] += f