root / ase / calculators / emt.py @ 4
Historique | Voir | Annoter | Télécharger (11,07 ko)
| 1 |
"""Effective medium theory potential."""
|
|---|---|
| 2 |
|
| 3 |
from math import sqrt, exp, log, pi |
| 4 |
|
| 5 |
import numpy as np |
| 6 |
import sys |
| 7 |
|
| 8 |
from ase.data import atomic_numbers, chemical_symbols |
| 9 |
from ase.units import Bohr |
| 10 |
|
| 11 |
|
| 12 |
parameters = {
|
| 13 |
# E0 s0 V0 eta2 kappa lambda n0
|
| 14 |
# eV bohr eV bohr^-1 bohr^-1 bohr^-1 bohr^-3
|
| 15 |
'H': (-2.21, 0.71, 2.132, 1.652, 2.790, 1.892, 0.00547, 'dimer'), |
| 16 |
'Al': (-3.28, 3.00, 1.493, 1.240, 2.000, 1.169, 0.00700, 'fcc'), |
| 17 |
'Cu': (-3.51, 2.67, 2.476, 1.652, 2.740, 1.906, 0.00910, 'fcc'), |
| 18 |
'Ag': (-2.96, 3.01, 2.132, 1.652, 2.790, 1.892, 0.00547, 'fcc'), |
| 19 |
'Au': (-3.80, 3.00, 2.321, 1.674, 2.873, 2.182, 0.00703, 'fcc'), |
| 20 |
'Ni': (-4.44, 2.60, 3.673, 1.669, 2.757, 1.948, 0.01030, 'fcc'), |
| 21 |
'Pd': (-3.90, 2.87, 2.773, 1.818, 3.107, 2.155, 0.00688, 'fcc'), |
| 22 |
'Pt': (-5.85, 2.90, 4.067, 1.812, 3.145, 2.192, 0.00802, 'fcc'), |
| 23 |
'C': (-1.97, 1.18, 0.132, 3.652, 5.790, 2.892, 0.01322, 'dimer'), |
| 24 |
'N': (-4.97, 1.18, 0.132, 2.652, 3.790, 2.892, 0.01222, 'dimer'), |
| 25 |
'O': (-2.97, 1.25, 2.132, 3.652, 5.790, 4.892, 0.00850, 'dimer')} |
| 26 |
|
| 27 |
beta = 1.809 #(16 * pi / 3)**(1.0 / 3) / 2**0.5 |
| 28 |
# but preserve historical rounding.
|
| 29 |
|
| 30 |
#eta1 = 0.5 / Bohr
|
| 31 |
|
| 32 |
class EMT: |
| 33 |
disabled = False # Set to True to disable (asap does this). |
| 34 |
|
| 35 |
def __init__(self, fix_alloys=False): |
| 36 |
self.energy = None |
| 37 |
self.fix_alloys = fix_alloys
|
| 38 |
if self.disabled: |
| 39 |
print >>sys.stderr, """ |
| 40 |
ase.EMT has been disabled by Asap. Most likely, you
|
| 41 |
intended to use Asap's EMT calculator, but accidentally
|
| 42 |
imported ase's EMT calculator after Asap's. This could
|
| 43 |
happen if your script contains the lines
|
| 44 |
|
| 45 |
from asap3 import *
|
| 46 |
from ase.calculators.emt import EMT
|
| 47 |
Swap the two lines to solve the problem.
|
| 48 |
|
| 49 |
(or 'from ase import *' in older versions of ASE.) Swap
|
| 50 |
the two lines to solve the problem.
|
| 51 |
|
| 52 |
In the UNLIKELY event that you actually wanted to use
|
| 53 |
ase.calculators.emt.EMT although asap3 is loaded into memory, please
|
| 54 |
reactivate it with the command
|
| 55 |
ase.calculators.emt.EMT.disabled = False
|
| 56 |
"""
|
| 57 |
raise RuntimeError('ase.EMT has been disabled. ' + |
| 58 |
'See message printed above.')
|
| 59 |
|
| 60 |
def get_spin_polarized(self): |
| 61 |
return False |
| 62 |
|
| 63 |
def initialize(self, atoms): |
| 64 |
self.par = {}
|
| 65 |
self.rc = 0.0 |
| 66 |
self.numbers = atoms.get_atomic_numbers()
|
| 67 |
maxseq = 0.0
|
| 68 |
seen = {}
|
| 69 |
for Z in self.numbers: |
| 70 |
if Z not in seen: |
| 71 |
seen[Z] = True
|
| 72 |
ss = parameters[chemical_symbols[Z]][1] * Bohr
|
| 73 |
if maxseq < ss:
|
| 74 |
maxseq = ss |
| 75 |
rc = self.rc = beta * maxseq * 0.5 * (sqrt(3) + sqrt(4)) |
| 76 |
rr = rc * 2 * sqrt(4) / (sqrt(3) + sqrt(4)) |
| 77 |
self.acut = np.log(9999.0) / (rr - rc) |
| 78 |
for Z in self.numbers: |
| 79 |
if Z not in self.par: |
| 80 |
p = parameters[chemical_symbols[Z]] |
| 81 |
s0 = p[1] * Bohr
|
| 82 |
eta2 = p[3] / Bohr
|
| 83 |
kappa = p[4] / Bohr
|
| 84 |
x = eta2 * beta * s0 |
| 85 |
gamma1 = 0.0
|
| 86 |
gamma2 = 0.0
|
| 87 |
if p[7] == 'fcc': |
| 88 |
for i, n in enumerate([12, 6, 24, 12]): |
| 89 |
r = s0 * beta * sqrt(i + 1)
|
| 90 |
x = n / (12 * (1.0 + exp(self.acut * (r - rc)))) |
| 91 |
gamma1 += x * exp(-eta2 * (r - beta * s0)) |
| 92 |
gamma2 += x * exp(-kappa / beta * (r - beta * s0)) |
| 93 |
elif p[7] == 'dimer': |
| 94 |
r = s0 * beta |
| 95 |
n = 1
|
| 96 |
x = n / (12 * (1.0 + exp(self.acut * (r - rc)))) |
| 97 |
gamma1 += x * exp(-eta2 * (r - beta * s0)) |
| 98 |
gamma2 += x * exp(-kappa / beta * (r - beta * s0)) |
| 99 |
else:
|
| 100 |
raise RuntimeError |
| 101 |
|
| 102 |
self.par[Z] = {'E0': p[0], |
| 103 |
's0': s0,
|
| 104 |
'V0': p[2], |
| 105 |
'eta2': eta2,
|
| 106 |
'kappa': kappa,
|
| 107 |
'lambda': p[5] / Bohr, |
| 108 |
'n0': p[6] / Bohr**3, |
| 109 |
'rc': rc,
|
| 110 |
'gamma1': gamma1,
|
| 111 |
'gamma2': gamma2}
|
| 112 |
#if rc + 0.5 > self.rc:
|
| 113 |
# self.rc = rc + 0.5
|
| 114 |
|
| 115 |
self.ksi = {}
|
| 116 |
for s1, p1 in self.par.items(): |
| 117 |
self.ksi[s1] = {}
|
| 118 |
for s2, p2 in self.par.items(): |
| 119 |
#self.ksi[s1][s2] = (p2['n0'] / p1['n0'] *
|
| 120 |
# exp(eta1 * (p1['s0'] - p2['s0'])))
|
| 121 |
self.ksi[s1][s2] = p2['n0'] / p1['n0'] |
| 122 |
|
| 123 |
self.forces = np.empty((len(atoms), 3)) |
| 124 |
self.sigma1 = np.empty(len(atoms)) |
| 125 |
self.deds = np.empty(len(atoms)) |
| 126 |
|
| 127 |
def update(self, atoms): |
| 128 |
if (self.energy is None or |
| 129 |
len(self.numbers) != len(atoms) or |
| 130 |
(self.numbers != atoms.get_atomic_numbers()).any()):
|
| 131 |
self.initialize(atoms)
|
| 132 |
self.calculate(atoms)
|
| 133 |
elif ((self.positions != atoms.get_positions()).any() or |
| 134 |
(self.pbc != atoms.get_pbc()).any() or |
| 135 |
(self.cell != atoms.get_cell()).any()):
|
| 136 |
self.calculate(atoms)
|
| 137 |
|
| 138 |
def calculation_required(self, atoms, quantities): |
| 139 |
if len(quantities) == 0: |
| 140 |
return False |
| 141 |
|
| 142 |
return (self.energy is None or |
| 143 |
len(self.numbers) != len(atoms) or |
| 144 |
(self.numbers != atoms.get_atomic_numbers()).any() or |
| 145 |
(self.positions != atoms.get_positions()).any() or |
| 146 |
(self.pbc != atoms.get_pbc()).any() or |
| 147 |
(self.cell != atoms.get_cell()).any())
|
| 148 |
|
| 149 |
def get_potential_energy(self, atoms): |
| 150 |
self.update(atoms)
|
| 151 |
return self.energy |
| 152 |
|
| 153 |
def get_numeric_forces(self, atoms): |
| 154 |
self.update(atoms)
|
| 155 |
p = atoms.positions |
| 156 |
p0 = p.copy() |
| 157 |
forces = np.empty_like(p) |
| 158 |
eps = 0.0001
|
| 159 |
for a in range(len(p)): |
| 160 |
for c in range(3): |
| 161 |
p[a, c] += eps |
| 162 |
self.calculate(atoms)
|
| 163 |
de = self.energy
|
| 164 |
p[a, c] -= 2 * eps
|
| 165 |
self.calculate(atoms)
|
| 166 |
de -= self.energy
|
| 167 |
p[a, c] += eps |
| 168 |
forces[a, c] = -de / (2 * eps)
|
| 169 |
p[:] = p0 |
| 170 |
return forces
|
| 171 |
|
| 172 |
def get_forces(self, atoms): |
| 173 |
self.update(atoms)
|
| 174 |
return self.forces.copy() |
| 175 |
|
| 176 |
def get_stress(self, atoms): |
| 177 |
raise NotImplementedError |
| 178 |
|
| 179 |
def calculate(self, atoms): |
| 180 |
self.positions = atoms.get_positions().copy()
|
| 181 |
self.cell = atoms.get_cell().copy()
|
| 182 |
self.pbc = atoms.get_pbc().copy()
|
| 183 |
|
| 184 |
icell = np.linalg.inv(self.cell)
|
| 185 |
scaled = np.dot(self.positions, icell)
|
| 186 |
N = [] |
| 187 |
for i in range(3): |
| 188 |
if self.pbc[i]: |
| 189 |
scaled[:, i] %= 1.0
|
| 190 |
v = icell[:, i] |
| 191 |
h = 1 / sqrt(np.dot(v, v))
|
| 192 |
N.append(int(self.rc / h) + 1) |
| 193 |
else:
|
| 194 |
N.append(0)
|
| 195 |
|
| 196 |
R = np.dot(scaled, self.cell)
|
| 197 |
|
| 198 |
self.energy = 0.0 |
| 199 |
self.sigma1[:] = 0.0 |
| 200 |
self.forces[:] = 0.0 |
| 201 |
|
| 202 |
N1, N2, N3 = N |
| 203 |
natoms = len(atoms)
|
| 204 |
for i1 in range(-N1, N1 + 1): |
| 205 |
for i2 in range(-N2, N2 + 1): |
| 206 |
for i3 in range(-N3, N3 + 1): |
| 207 |
C = np.dot((i1, i2, i3), self.cell)
|
| 208 |
Q = R + C |
| 209 |
c = (i1 == 0 and i2 == 0 and i3 == 0) |
| 210 |
for a1 in range(natoms): |
| 211 |
Z1 = self.numbers[a1]
|
| 212 |
p1 = self.par[Z1]
|
| 213 |
ksi = self.ksi[Z1]
|
| 214 |
for a2 in range(natoms): |
| 215 |
if c and a2 == a1: |
| 216 |
continue
|
| 217 |
d = Q[a2] - R[a1] |
| 218 |
r = sqrt(np.dot(d, d)) |
| 219 |
if r < self.rc + 0.5: |
| 220 |
Z2 = self.numbers[a2]
|
| 221 |
if self.fix_alloys: |
| 222 |
p2 = self.par[Z2]
|
| 223 |
else:
|
| 224 |
p2 = p1 # Old buggy behaviour
|
| 225 |
self.interact1(a1, a2, d, r, p1, p2, ksi[Z2])
|
| 226 |
|
| 227 |
for a in range(natoms): |
| 228 |
Z = self.numbers[a]
|
| 229 |
p = self.par[Z]
|
| 230 |
try:
|
| 231 |
ds = -log(self.sigma1[a] / 12) / (beta * p['eta2']) |
| 232 |
except (OverflowError, ValueError): |
| 233 |
self.deds[a] = 0.0 |
| 234 |
self.energy -= p['E0'] |
| 235 |
continue
|
| 236 |
x = p['lambda'] * ds
|
| 237 |
y = exp(-x) |
| 238 |
z = 6 * p['V0'] * exp(-p['kappa'] * ds) |
| 239 |
self.deds[a] = ((x * y * p['E0'] * p['lambda'] + p['kappa'] * z) / |
| 240 |
(self.sigma1[a] * beta * p['eta2'])) |
| 241 |
e = p['E0'] * ((1 + x) * y - 1) + z |
| 242 |
self.energy += p['E0'] * ((1 + x) * y - 1) + z |
| 243 |
|
| 244 |
for i1 in range(-N1, N1 + 1): |
| 245 |
for i2 in range(-N2, N2 + 1): |
| 246 |
for i3 in range(-N3, N3 + 1): |
| 247 |
C = np.dot((i1, i2, i3), self.cell)
|
| 248 |
Q = R + C |
| 249 |
c = (i1 == 0 and i2 == 0 and i3 == 0) |
| 250 |
for a1 in range(natoms): |
| 251 |
Z1 = self.numbers[a1]
|
| 252 |
p1 = self.par[Z1]
|
| 253 |
ksi = self.ksi[Z1]
|
| 254 |
for a2 in range(natoms): |
| 255 |
if c and a2 == a1: |
| 256 |
continue
|
| 257 |
d = Q[a2] - R[a1] |
| 258 |
r = sqrt(np.dot(d, d)) |
| 259 |
if r < self.rc + 0.5: |
| 260 |
Z2 = self.numbers[a2]
|
| 261 |
if self.fix_alloys: |
| 262 |
p2 = self.par[Z2]
|
| 263 |
else:
|
| 264 |
p2 = p1 # Old buggy behaviour
|
| 265 |
self.interact2(a1, a2, d, r, p1, p2, ksi[Z2])
|
| 266 |
|
| 267 |
def interact1(self, a1, a2, d, r, p1, p2, ksi): |
| 268 |
x = exp(self.acut * (r - self.rc)) |
| 269 |
theta = 1.0 / (1.0 + x) |
| 270 |
y = (0.5 * p1['V0'] * exp(-p2['kappa'] * (r / beta - p2['s0'])) * |
| 271 |
ksi / p1['gamma2'] * theta)
|
| 272 |
self.energy -= y
|
| 273 |
f = y * (p2['kappa'] / beta + self.acut * theta * x) * d / r |
| 274 |
self.forces[a1] += f
|
| 275 |
self.forces[a2] -= f
|
| 276 |
self.sigma1[a1] += (exp(-p2['eta2'] * (r - beta * p2['s0'])) * |
| 277 |
ksi * theta / p1['gamma1'])
|
| 278 |
|
| 279 |
def interact2(self, a1, a2, d, r, p1, p2, ksi): |
| 280 |
x = exp(self.acut * (r - self.rc)) |
| 281 |
theta = 1.0 / (1.0 + x) |
| 282 |
y = (exp(-p2['eta2'] * (r - beta * p2['s0'])) * |
| 283 |
ksi / p1['gamma1'] * theta * self.deds[a1]) |
| 284 |
f = y * (p2['eta2'] + self.acut * theta * x) * d / r |
| 285 |
self.forces[a1] -= f
|
| 286 |
self.forces[a2] += f
|
| 287 |
|