Coverage for ase / eos.py: 77.16%
197 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-21 15:52 +0000
1# fmt: off
3import warnings
5import numpy as np
7from ase.units import kJ
9eos_names = [
10 'sj',
11 'taylor',
12 'murnaghan',
13 'birch',
14 'birchmurnaghan',
15 'pouriertarantola',
16 'vinet',
17 'anton-schmidt',
18 'p3',
19]
22def taylor(V, E0, beta, alpha, V0):
23 'Taylor Expansion up to 3rd order about V0'
25 E = E0 + beta / 2 * (V - V0)**2 / V0 + alpha / 6 * (V - V0)**3 / V0
26 return E
29def murnaghan(V, E0, B0, BP, V0):
30 'From PRB 28,5480 (1983'
32 E = E0 + B0 * V / BP * (((V0 / V)**BP) / (BP - 1) + 1) - V0 * B0 / (BP - 1)
33 return E
36def birch(V, E0, B0, BP, V0):
37 """
38 From Intermetallic compounds: Principles and Practice, Vol. I: Principles
39 Chapter 9 pages 195-210 by M. Mehl. B. Klein, D. Papaconstantopoulos
40 paper downloaded from Web
42 case where n=0
43 """
45 E = (E0 +
46 9 / 8 * B0 * V0 * ((V0 / V)**(2 / 3) - 1)**2 +
47 9 / 16 * B0 * V0 * (BP - 4) * ((V0 / V)**(2 / 3) - 1)**3)
48 return E
51def birchmurnaghan(V, E0, B0, BP, V0):
52 """
53 BirchMurnaghan equation from PRB 70, 224107
54 Eq. (3) in the paper. Note that there's a typo in the paper and it uses
55 inversed expression for eta.
56 """
58 eta = (V0 / V)**(1 / 3)
59 E = E0 + 9 * B0 * V0 / 16 * (eta**2 - 1)**2 * (
60 6 + BP * (eta**2 - 1) - 4 * eta**2)
61 return E
64def check_birchmurnaghan():
65 from sympy import Rational, diff, simplify, symbols
66 v, b, bp, v0 = symbols('v b bp v0')
67 x = (v0 / v)**Rational(2, 3)
68 e = 9 * b * v0 * (x - 1)**2 * (6 + bp * (x - 1) - 4 * x) / 16
69 print(e)
70 B = diff(e, v, 2) * v
71 BP = -v * diff(B, v) / b
72 print(simplify(B.subs(v, v0)))
73 print(simplify(BP.subs(v, v0)))
76def pouriertarantola(V, E0, B0, BP, V0):
77 'Pourier-Tarantola equation from PRB 70, 224107'
79 eta = (V / V0)**(1 / 3)
80 squiggle = -3 * np.log(eta)
82 E = E0 + B0 * V0 * squiggle**2 / 6 * (3 + squiggle * (BP - 2))
83 return E
86def vinet(V, E0, B0, BP, V0):
87 'Vinet equation from PRB 70, 224107'
89 eta = (V / V0)**(1 / 3)
91 E = (E0 + 2 * B0 * V0 / (BP - 1)**2 *
92 (2 - (5 + 3 * BP * (eta - 1) - 3 * eta) *
93 np.exp(-3 * (BP - 1) * (eta - 1) / 2)))
94 return E
97def anton_schmidt(
98 v: np.ndarray,
99 e0: float,
100 b0: float,
101 bp0: float,
102 v0: float,
103) -> np.ndarray:
104 r"""Anton-Schmidt EOS [1]_.
106 In the original formula, :math:`E_\infty` (energy at infinite separation)
107 is a parameter, which is given by [1]_
109 .. math:: E_\infty = E_0 + \frac{B_0 V_0}{(n + 1)^2}.
111 :math:`n` is related to the Grüneisen parameter [1]_ as
113 .. math:: n = -1/6 - \gamma.
115 Taking the Grüneisen parameters in Ref. [2]_, :math:`n` should be about -2.
117 :math:`n` is also given with :math:`B'_{0}` by [3]_:
119 .. math:: n = -\\frac{1}{2} B'_{0}.
121 References
122 ----------
123 .. [1] H. Anton and P. C. Schmidt, Intermetallics 5, 449 (1997).
124 :doi:`10.1016/s0966-9795(97)00017-4`
125 .. [2] B. Mayer, H. Anton, E. Bott, M. Methfessel, J. Sticht, J. Harris,
126 and P. C. Schmidt, Intermetallics 11, 23 (2003).
127 :doi:`10.1016/s0966-9795(02)00127-9`
128 .. [3] A. Otero-de-la-Roza and V. Luaña,
129 Comput. Phys. Commun. 182, 1708 (2011).
130 :doi:`10.1016/j.cpc.2011.04.016`
132 """
133 n = -0.5 * bp0
134 e_infinity = e0 + b0 * v0 / (n + 1) ** 2
135 e = b0 * v0 / (n + 1) * (v / v0)**(n + 1) * (np.log(v / v0) - 1 / (n + 1))
136 return e + e_infinity
139def p3(V, c0, c1, c2, c3):
140 'polynomial fit'
142 E = c0 + c1 * V + c2 * V**2 + c3 * V**3
143 return E
146def parabola(x, a, b, c):
147 """parabola polynomial function
149 this function is used to fit the data to get good guesses for
150 the equation of state fits
152 a 4th order polynomial fit to get good guesses for
153 was not a good idea because for noisy data the fit is too wiggly
154 2nd order seems to be sufficient, and guarantees a single minimum"""
156 return a + b * x + c * x**2
159class EquationOfState:
160 """Fit equation of state for bulk systems.
162 Parameters
163 ----------
164 eos : str
165 Name of the equation of state to fit. Must be one of:
167 - ``sj`` (default)
168 A third order inverse polynomial fit 10.1103/PhysRevB.67.026103
170 ::
172 2 3 -1/3
173 E(V) = c + c t + c t + c t , t = V
174 0 1 2 3
176 - ``taylor``
177 A third order Taylor series expansion about the minimum volume
179 - ``murnaghan``
180 PRB 28, 5480 (1983)
182 - ``birch``
183 Intermetallic compounds: Principles and Practice,
184 Vol I: Principles. pages 195-210
186 - ``birchmurnaghan``
187 PRB 70, 224107
189 - ``pouriertarantola``
190 PRB 70, 224107
192 - ``vinet``
193 PRB 70, 224107
195 - ``anton-schmidt``
196 H. Anton and P. C. Schmidt, Intermetallics 5, 449 (1997).
197 :doi:`10.1016/s0966-9795(97)00017-4`
199 - ``p3``
200 A third order polynomial fit
202 Examples
203 --------
204 >>> from ase.build import bulk
205 >>> from ase.calculators.emt import EMT
206 >>> atoms = bulk('Al')
207 >>> atoms.calc = EMT()
208 >>> cell = atoms.cell.copy()
209 >>> volumes = []
210 >>> energies = []
211 >>> for x in np.linspace(0.97, 1.01, 5):
212 ... atoms.set_cell(cell * x, scale_atoms=True)
213 ... volumes.append(atoms.get_volume())
214 ... energies.append(atoms.get_potential_energy())
215 >>> eos = EquationOfState(volumes, energies, eos='murnaghan')
216 >>> v0, e0, b0 = eos.fit()
217 >>> ax = eos.plot(show=False)
219 """
221 def __init__(self, volumes, energies, eos='sj'):
222 self.v = np.array(volumes)
223 self.e = np.array(energies)
225 old2new = {
226 'sjeos': 'sj',
227 'antonschmidt': 'anton-schmidt',
228 }
229 for old_name, new_name in old2new.items():
230 if eos == old_name:
231 eos = new_name
232 msg = (
233 f'`{old_name}` is an old EOS name. '
234 f'Use `{new_name}` instead. '
235 'The old EOS names will be removed in ASE 3.30.0.'
236 )
237 warnings.warn(msg, FutureWarning)
238 self.eos_string = eos
239 self.v0 = None
241 def fit(self, warn=True):
242 """Calculate volume, energy, and bulk modulus.
244 Returns the optimal volume, the minimum energy, and the bulk
245 modulus. Notice that the ASE units for the bulk modulus is
246 eV/Angstrom^3 - to get the value in GPa, do this::
248 v0, e0, B = eos.fit()
249 print(B / kJ * 1.0e24, 'GPa')
251 """
252 from scipy.optimize import curve_fit
254 if self.eos_string == 'sj':
255 return self.fit_sjeos()
257 self.func = globals()[self.eos_string.replace('-', '_')]
259 p0 = [min(self.e), 1, 1]
260 popt, _pcov = curve_fit(parabola, self.v, self.e, p0)
262 parabola_parameters = popt
263 # Here I just make sure the minimum is bracketed by the volumes
264 # this if for the solver
265 minvol = min(self.v)
266 maxvol = max(self.v)
268 # the minimum of the parabola is at dE/dV = 0, or 2 * c V +b =0
269 c = parabola_parameters[2]
270 b = parabola_parameters[1]
271 a = parabola_parameters[0]
272 parabola_vmin = -b / 2 / c
274 # evaluate the parabola at the minimum to estimate the groundstate
275 # energy
276 E0 = parabola(parabola_vmin, a, b, c)
277 # estimate the bulk modulus from Vo * E''. E'' = 2 * c
278 B0 = 2 * c * parabola_vmin
280 # (∂B/∂P)_V0 = (∂B/∂V)_V0 / (∂P/∂V)_V0 (4.0 for the 2nd-order BM)
281 bp0 = 4.0
283 initial_guess = [E0, B0, bp0, parabola_vmin]
285 # now fit the equation of state
286 p0 = initial_guess
287 popt, _pcov = curve_fit(self.func, self.v, self.e, p0)
289 self.eos_parameters = popt
291 if self.eos_string == 'p3':
292 c0, c1, c2, c3 = self.eos_parameters
293 # find minimum E in E = c0 + c1 * V + c2 * V**2 + c3 * V**3
294 # dE/dV = c1+ 2 * c2 * V + 3 * c3 * V**2 = 0
295 # solve by quadratic formula with the positive root
297 a = 3 * c3
298 b = 2 * c2
299 c = c1
301 self.v0 = (-b + np.sqrt(b**2 - 4 * a * c)) / (2 * a)
302 self.e0 = p3(self.v0, c0, c1, c2, c3)
303 self.B = (2 * c2 + 6 * c3 * self.v0) * self.v0
304 else:
305 self.v0 = self.eos_parameters[3]
306 self.e0 = self.eos_parameters[0]
307 self.B = self.eos_parameters[1]
309 if warn and not (minvol < self.v0 < maxvol):
310 warnings.warn(
311 'The minimum volume of your fit is not in '
312 'your volumes. You may not have a minimum in your dataset!')
314 return self.v0, self.e0, self.B
316 def getplotdata(self):
317 if self.v0 is None:
318 self.fit()
320 x = np.linspace(min(self.v), max(self.v), 100)
321 if self.eos_string == 'sj':
322 y = self.fit0(x**-(1 / 3))
323 else:
324 y = self.func(x, *self.eos_parameters)
326 return self.eos_string, self.e0, self.v0, self.B, x, y, self.v, self.e
328 def plot(self, filename=None, show=False, ax=None):
329 """Plot fitted energy curve.
331 Uses Matplotlib to plot the energy curve. Use *show=True* to
332 show the figure and *filename='abc.png'* or
333 *filename='abc.eps'* to save the figure to a file."""
335 import matplotlib.pyplot as plt
337 plotdata = self.getplotdata()
339 ax = plot(*plotdata, ax=ax)
341 if show:
342 plt.show()
343 if filename is not None:
344 fig = ax.get_figure()
345 fig.savefig(filename)
346 return ax
348 def fit_sjeos(self):
349 """Calculate volume, energy, and bulk modulus.
351 Returns the optimal volume, the minimum energy, and the bulk
352 modulus. Notice that the ASE units for the bulk modulus is
353 eV/Angstrom^3 - to get the value in GPa, do this::
355 v0, e0, B = eos.fit()
356 print(B / kJ * 1.0e24, 'GPa')
358 """
360 fit0 = np.poly1d(np.polyfit(self.v**-(1 / 3), self.e, 3))
361 fit1 = np.polyder(fit0, 1)
362 fit2 = np.polyder(fit1, 1)
364 self.v0 = None
365 for t in np.roots(fit1):
366 if isinstance(t, float) and t > 0 and fit2(t) > 0:
367 self.v0 = t**-3
368 break
370 if self.v0 is None:
371 raise ValueError('No minimum!')
373 self.e0 = fit0(t)
374 self.B = t**5 * fit2(t) / 9
375 self.fit0 = fit0
377 return self.v0, self.e0, self.B
380def plot(eos_string, e0, v0, B, x, y, v, e, ax=None):
381 if ax is None:
382 import matplotlib.pyplot as plt
383 ax = plt.gca()
385 ax.plot(x, y, ls='-', color='C3') # By default red line
386 ax.plot(v, e, ls='', marker='o', mec='C0',
387 mfc='C0') # By default blue marker
389 try:
390 ax.set_xlabel('volume [Å$^3$]')
391 ax.set_ylabel('energy [eV]')
392 ax.set_title('%s: E: %.3f eV, V: %.3f Å$^3$, B: %.3f GPa' %
393 (eos_string, e0, v0,
394 B / kJ * 1.e24))
396 except ImportError: # XXX what would cause this error? LaTeX?
397 warnings.warn('Could not use LaTeX formatting')
398 ax.set_xlabel('volume [L(length)^3]')
399 ax.set_ylabel('energy [E(energy)]')
400 ax.set_title('%s: E: %.3f E, V: %.3f L^3, B: %.3e E/L^3' %
401 (eos_string, e0, v0, B))
403 return ax
406def calculate_eos(atoms, npoints=5, eps=0.04, trajectory=None, callback=None):
407 """Calculate equation-of-state.
409 atoms: Atoms object
410 System to calculate EOS for. Must have a calculator attached.
411 npoints: int
412 Number of points.
413 eps: float
414 Variation in volume from v0*(1-eps) to v0*(1+eps).
415 trajectory: Trjectory object or str
416 Write configurations to a trajectory file.
417 callback: function
418 Called after every energy calculation.
420 >>> from ase.build import bulk
421 >>> from ase.calculators.emt import EMT
422 >>> from ase.eos import calculate_eos
424 >>> a = bulk('Cu', 'fcc', a=3.6)
425 >>> a.calc = EMT()
426 >>> eos = calculate_eos(a, trajectory='Cu.traj')
427 >>> v, e, B = eos.fit()
428 >>> a = (4 * v)**(1 / 3.0)
429 >>> print('{0:.6f}'.format(a))
430 3.589825
431 """
433 # Save original positions and cell:
434 p0 = atoms.get_positions()
435 c0 = atoms.get_cell()
437 if isinstance(trajectory, str):
438 from ase.io import Trajectory
439 trajectory = Trajectory(trajectory, 'w', atoms)
441 if trajectory is not None:
442 trajectory.set_description({'type': 'eos',
443 'npoints': npoints,
444 'eps': eps})
446 try:
447 energies = []
448 volumes = []
449 for x in np.linspace(1 - eps, 1 + eps, npoints)**(1 / 3):
450 atoms.set_cell(x * c0, scale_atoms=True)
451 volumes.append(atoms.get_volume())
452 energies.append(atoms.get_potential_energy())
453 if callback:
454 callback()
455 if trajectory is not None:
456 trajectory.write()
457 return EquationOfState(volumes, energies)
458 finally:
459 atoms.cell = c0
460 atoms.positions = p0
461 if trajectory is not None:
462 trajectory.close()
465class CLICommand:
466 """Calculate EOS from one or more trajectory files.
468 See https://ase-lib.org/tutorials/eos/eos.html for
469 more information.
470 """
472 @staticmethod
473 def add_arguments(parser):
474 parser.add_argument('trajectories', nargs='+', metavar='trajectory')
475 parser.add_argument('-p', '--plot', action='store_true',
476 help='Plot EOS fit. Default behaviour is '
477 'to write results of fit.')
478 parser.add_argument('-t', '--type', default='sj',
479 help='Type of fit. Must be one of {}.'
480 .format(', '.join(eos_names)))
482 @staticmethod
483 def run(args):
484 from ase.io import read
486 if not args.plot:
487 print('# filename '
488 'points volume energy bulk modulus')
489 print('# '
490 ' [Ang^3] [eV] [GPa]')
491 for name in args.trajectories:
492 if name == '-':
493 # Special case - used by ASE's GUI:
494 import pickle
495 import sys
496 v, e = pickle.load(sys.stdin.buffer)
497 else:
498 if '@' in name:
499 index = None
500 else:
501 index = ':'
502 images = read(name, index=index)
503 v = [atoms.get_volume() for atoms in images]
504 e = [atoms.get_potential_energy() for atoms in images]
505 eos = EquationOfState(v, e, args.type)
506 if args.plot:
507 eos.plot()
508 else:
509 try:
510 v0, e0, B = eos.fit()
511 except ValueError as ex:
512 print('{:30}{:2} {}'
513 .format(name, len(v), ex.message))
514 else:
515 print('{:30}{:2} {:10.3f}{:10.3f}{:14.3f}'
516 .format(name, len(v), v0, e0, B / kJ * 1.0e24))