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