import warnings
import numpy as np
from ase.units import kJ
eos_names = ['sj', 'taylor', 'murnaghan', 'birch', 'birchmurnaghan',
'pouriertarantola', 'vinet', 'antonschmidt', 'p3']
def taylor(V, E0, beta, alpha, V0):
'Taylor Expansion up to 3rd order about V0'
E = E0 + beta / 2 * (V - V0)**2 / V0 + alpha / 6 * (V - V0)**3 / V0
return E
def murnaghan(V, E0, B0, BP, V0):
'From PRB 28,5480 (1983'
E = E0 + B0 * V / BP * (((V0 / V)**BP) / (BP - 1) + 1) - V0 * B0 / (BP - 1)
return E
def birch(V, E0, B0, BP, V0):
"""
From Intermetallic compounds: Principles and Practice, Vol. I: Principles
Chapter 9 pages 195-210 by M. Mehl. B. Klein, D. Papaconstantopoulos
paper downloaded from Web
case where n=0
"""
E = (E0 +
9 / 8 * B0 * V0 * ((V0 / V)**(2 / 3) - 1)**2 +
9 / 16 * B0 * V0 * (BP - 4) * ((V0 / V)**(2 / 3) - 1)**3)
return E
def birchmurnaghan(V, E0, B0, BP, V0):
"""
BirchMurnaghan equation from PRB 70, 224107
Eq. (3) in the paper. Note that there's a typo in the paper and it uses
inversed expression for eta.
"""
eta = (V0 / V)**(1 / 3)
E = E0 + 9 * B0 * V0 / 16 * (eta**2 - 1)**2 * (
6 + BP * (eta**2 - 1) - 4 * eta**2)
return E
def check_birchmurnaghan():
from sympy import Rational, diff, simplify, symbols
v, b, bp, v0 = symbols('v b bp v0')
x = (v0 / v)**Rational(2, 3)
e = 9 * b * v0 * (x - 1)**2 * (6 + bp * (x - 1) - 4 * x) / 16
print(e)
B = diff(e, v, 2) * v
BP = -v * diff(B, v) / b
print(simplify(B.subs(v, v0)))
print(simplify(BP.subs(v, v0)))
def pouriertarantola(V, E0, B0, BP, V0):
'Pourier-Tarantola equation from PRB 70, 224107'
eta = (V / V0)**(1 / 3)
squiggle = -3 * np.log(eta)
E = E0 + B0 * V0 * squiggle**2 / 6 * (3 + squiggle * (BP - 2))
return E
def vinet(V, E0, B0, BP, V0):
'Vinet equation from PRB 70, 224107'
eta = (V / V0)**(1 / 3)
E = (E0 + 2 * B0 * V0 / (BP - 1)**2 *
(2 - (5 + 3 * BP * (eta - 1) - 3 * eta) *
np.exp(-3 * (BP - 1) * (eta - 1) / 2)))
return E
def antonschmidt(V, Einf, B, n, V0):
"""From Intermetallics 11, 23-32 (2003)
Einf should be E_infinity, i.e. infinite separation, but
according to the paper it does not provide a good estimate
of the cohesive energy. They derive this equation from an
empirical formula for the volume dependence of pressure,
E(vol) = E_inf + int(P dV) from V=vol to V=infinity
but the equation breaks down at large volumes, so E_inf
is not that meaningful
n should be about -2 according to the paper.
I find this equation does not fit volumetric data as well
as the other equtions do.
"""
E = B * V0 / (n + 1) * (V / V0)**(n + 1) * (np.log(V / V0) -
(1 / (n + 1))) + Einf
return E
def p3(V, c0, c1, c2, c3):
'polynomial fit'
E = c0 + c1 * V + c2 * V**2 + c3 * V**3
return E
def parabola(x, a, b, c):
"""parabola polynomial function
this function is used to fit the data to get good guesses for
the equation of state fits
a 4th order polynomial fit to get good guesses for
was not a good idea because for noisy data the fit is too wiggly
2nd order seems to be sufficient, and guarantees a single minimum"""
return a + b * x + c * x**2
[docs]class EquationOfState:
"""Fit equation of state for bulk systems.
The following equation is used::
sjeos (default)
A third order inverse polynomial fit 10.1103/PhysRevB.67.026103
::
2 3 -1/3
E(V) = c + c t + c t + c t , t = V
0 1 2 3
taylor
A third order Taylor series expansion about the minimum volume
murnaghan
PRB 28, 5480 (1983)
birch
Intermetallic compounds: Principles and Practice,
Vol I: Principles. pages 195-210
birchmurnaghan
PRB 70, 224107
pouriertarantola
PRB 70, 224107
vinet
PRB 70, 224107
antonschmidt
Intermetallics 11, 23-32 (2003)
p3
A third order polynomial fit
Use::
eos = EquationOfState(volumes, energies, eos='murnaghan')
v0, e0, B = eos.fit()
eos.plot(show=True)
"""
def __init__(self, volumes, energies, eos='sj'):
self.v = np.array(volumes)
self.e = np.array(energies)
if eos == 'sjeos':
eos = 'sj'
self.eos_string = eos
self.v0 = None
[docs] def fit(self, warn=True):
"""Calculate volume, energy, and bulk modulus.
Returns the optimal volume, the minimum energy, and the bulk
modulus. Notice that the ASE units for the bulk modulus is
eV/Angstrom^3 - to get the value in GPa, do this::
v0, e0, B = eos.fit()
print(B / kJ * 1.0e24, 'GPa')
"""
from scipy.optimize import curve_fit
if self.eos_string == 'sj':
return self.fit_sjeos()
self.func = globals()[self.eos_string]
p0 = [min(self.e), 1, 1]
popt, pcov = curve_fit(parabola, self.v, self.e, p0)
parabola_parameters = popt
# Here I just make sure the minimum is bracketed by the volumes
# this if for the solver
minvol = min(self.v)
maxvol = max(self.v)
# the minimum of the parabola is at dE/dV = 0, or 2 * c V +b =0
c = parabola_parameters[2]
b = parabola_parameters[1]
a = parabola_parameters[0]
parabola_vmin = -b / 2 / c
# evaluate the parabola at the minimum to estimate the groundstate
# energy
E0 = parabola(parabola_vmin, a, b, c)
# estimate the bulk modulus from Vo * E''. E'' = 2 * c
B0 = 2 * c * parabola_vmin
if self.eos_string == 'antonschmidt':
BP = -2
else:
BP = 4
initial_guess = [E0, B0, BP, parabola_vmin]
# now fit the equation of state
p0 = initial_guess
popt, pcov = curve_fit(self.func, self.v, self.e, p0)
self.eos_parameters = popt
if self.eos_string == 'p3':
c0, c1, c2, c3 = self.eos_parameters
# find minimum E in E = c0 + c1 * V + c2 * V**2 + c3 * V**3
# dE/dV = c1+ 2 * c2 * V + 3 * c3 * V**2 = 0
# solve by quadratic formula with the positive root
a = 3 * c3
b = 2 * c2
c = c1
self.v0 = (-b + np.sqrt(b**2 - 4 * a * c)) / (2 * a)
self.e0 = p3(self.v0, c0, c1, c2, c3)
self.B = (2 * c2 + 6 * c3 * self.v0) * self.v0
else:
self.v0 = self.eos_parameters[3]
self.e0 = self.eos_parameters[0]
self.B = self.eos_parameters[1]
if warn and not (minvol < self.v0 < maxvol):
warnings.warn(
'The minimum volume of your fit is not in '
'your volumes. You may not have a minimum in your dataset!')
return self.v0, self.e0, self.B
def getplotdata(self):
if self.v0 is None:
self.fit()
x = np.linspace(min(self.v), max(self.v), 100)
if self.eos_string == 'sj':
y = self.fit0(x**-(1 / 3))
else:
y = self.func(x, *self.eos_parameters)
return self.eos_string, self.e0, self.v0, self.B, x, y, self.v, self.e
[docs] def plot(self, filename=None, show=False, ax=None):
"""Plot fitted energy curve.
Uses Matplotlib to plot the energy curve. Use *show=True* to
show the figure and *filename='abc.png'* or
*filename='abc.eps'* to save the figure to a file."""
import matplotlib.pyplot as plt
plotdata = self.getplotdata()
ax = plot(*plotdata, ax=ax)
if show:
plt.show()
if filename is not None:
fig = ax.get_figure()
fig.savefig(filename)
return ax
def fit_sjeos(self):
"""Calculate volume, energy, and bulk modulus.
Returns the optimal volume, the minimum energy, and the bulk
modulus. Notice that the ASE units for the bulk modulus is
eV/Angstrom^3 - to get the value in GPa, do this::
v0, e0, B = eos.fit()
print(B / kJ * 1.0e24, 'GPa')
"""
fit0 = np.poly1d(np.polyfit(self.v**-(1 / 3), self.e, 3))
fit1 = np.polyder(fit0, 1)
fit2 = np.polyder(fit1, 1)
self.v0 = None
for t in np.roots(fit1):
if isinstance(t, float) and t > 0 and fit2(t) > 0:
self.v0 = t**-3
break
if self.v0 is None:
raise ValueError('No minimum!')
self.e0 = fit0(t)
self.B = t**5 * fit2(t) / 9
self.fit0 = fit0
return self.v0, self.e0, self.B
def plot(eos_string, e0, v0, B, x, y, v, e, ax=None):
if ax is None:
import matplotlib.pyplot as plt
ax = plt.gca()
ax.plot(x, y, ls='-', color='C3') # By default red line
ax.plot(v, e, ls='', marker='o', mec='C0',
mfc='C0') # By default blue marker
try:
ax.set_xlabel('volume [Å$^3$]')
ax.set_ylabel('energy [eV]')
ax.set_title('%s: E: %.3f eV, V: %.3f Å$^3$, B: %.3f GPa' %
(eos_string, e0, v0,
B / kJ * 1.e24))
except ImportError: # XXX what would cause this error? LaTeX?
import warnings
warnings.warn('Could not use LaTeX formatting')
ax.set_xlabel('volume [L(length)^3]')
ax.set_ylabel('energy [E(energy)]')
ax.set_title('%s: E: %.3f E, V: %.3f L^3, B: %.3e E/L^3' %
(eos_string, e0, v0, B))
return ax
[docs]def calculate_eos(atoms, npoints=5, eps=0.04, trajectory=None, callback=None):
"""Calculate equation-of-state.
atoms: Atoms object
System to calculate EOS for. Must have a calculator attached.
npoints: int
Number of points.
eps: float
Variation in volume from v0*(1-eps) to v0*(1+eps).
trajectory: Trjectory object or str
Write configurations to a trajectory file.
callback: function
Called after every energy calculation.
>>> from ase.build import bulk
>>> from ase.calculators.emt import EMT
>>> from ase.eos import calculate_eos
>>> a = bulk('Cu', 'fcc', a=3.6)
>>> a.calc = EMT()
>>> eos = calculate_eos(a, trajectory='Cu.traj')
>>> v, e, B = eos.fit()
>>> a = (4 * v)**(1 / 3.0)
>>> print('{0:.6f}'.format(a))
3.589825
"""
# Save original positions and cell:
p0 = atoms.get_positions()
c0 = atoms.get_cell()
if isinstance(trajectory, str):
from ase.io import Trajectory
trajectory = Trajectory(trajectory, 'w', atoms)
if trajectory is not None:
trajectory.set_description({'type': 'eos',
'npoints': npoints,
'eps': eps})
try:
energies = []
volumes = []
for x in np.linspace(1 - eps, 1 + eps, npoints)**(1 / 3):
atoms.set_cell(x * c0, scale_atoms=True)
volumes.append(atoms.get_volume())
energies.append(atoms.get_potential_energy())
if callback:
callback()
if trajectory is not None:
trajectory.write()
return EquationOfState(volumes, energies)
finally:
atoms.cell = c0
atoms.positions = p0
if trajectory is not None:
trajectory.close()
class CLICommand:
"""Calculate EOS from one or more trajectory files.
See https://wiki.fysik.dtu.dk/ase/tutorials/eos/eos.html for
more information.
"""
@staticmethod
def add_arguments(parser):
parser.add_argument('trajectories', nargs='+', metavar='trajectory')
parser.add_argument('-p', '--plot', action='store_true',
help='Plot EOS fit. Default behaviour is '
'to write results of fit.')
parser.add_argument('-t', '--type', default='sj',
help='Type of fit. Must be one of {}.'
.format(', '.join(eos_names)))
@staticmethod
def run(args):
from ase.io import read
if not args.plot:
print('# filename '
'points volume energy bulk modulus')
print('# '
' [Ang^3] [eV] [GPa]')
for name in args.trajectories:
if name == '-':
# Special case - used by ASE's GUI:
import pickle
import sys
v, e = pickle.load(sys.stdin.buffer)
else:
if '@' in name:
index = None
else:
index = ':'
images = read(name, index=index)
v = [atoms.get_volume() for atoms in images]
e = [atoms.get_potential_energy() for atoms in images]
eos = EquationOfState(v, e, args.type)
if args.plot:
eos.plot()
else:
try:
v0, e0, B = eos.fit()
except ValueError as ex:
print('{:30}{:2} {}'
.format(name, len(v), ex.message))
else:
print('{:30}{:2} {:10.3f}{:10.3f}{:14.3f}'
.format(name, len(v), v0, e0, B / kJ * 1.0e24))