r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import warnings

3from ase.units import kJ

5import numpy as np

8eos_names = ['sj', 'taylor', 'murnaghan', 'birch', 'birchmurnaghan',

9 'pouriertarantola', 'vinet', 'antonschmidt', 'p3']

12def taylor(V, E0, beta, alpha, V0):

13 'Taylor Expansion up to 3rd order about V0'

15 E = E0 + beta / 2 * (V - V0)**2 / V0 + alpha / 6 * (V - V0)**3 / V0

16 return E

19def murnaghan(V, E0, B0, BP, V0):

20 'From PRB 28,5480 (1983'

22 E = E0 + B0 * V / BP * (((V0 / V)**BP) / (BP - 1) + 1) - V0 * B0 / (BP - 1)

23 return E

26def birch(V, E0, B0, BP, V0):

27 """

28 From Intermetallic compounds: Principles and Practice, Vol. I: Principles

29 Chapter 9 pages 195-210 by M. Mehl. B. Klein, D. Papaconstantopoulos

32 case where n=0

33 """

35 E = (E0 +

36 9 / 8 * B0 * V0 * ((V0 / V)**(2 / 3) - 1)**2 +

37 9 / 16 * B0 * V0 * (BP - 4) * ((V0 / V)**(2 / 3) - 1)**3)

38 return E

41def birchmurnaghan(V, E0, B0, BP, V0):

42 """

43 BirchMurnaghan equation from PRB 70, 224107

44 Eq. (3) in the paper. Note that there's a typo in the paper and it uses

45 inversed expression for eta.

46 """

48 eta = (V0 / V)**(1 / 3)

49 E = E0 + 9 * B0 * V0 / 16 * (eta**2 - 1)**2 * (

50 6 + BP * (eta**2 - 1) - 4 * eta**2)

51 return E

54def check_birchmurnaghan():

55 from sympy import symbols, Rational, diff, simplify

56 v, b, bp, v0 = symbols('v b bp v0')

57 x = (v0 / v)**Rational(2, 3)

58 e = 9 * b * v0 * (x - 1)**2 * (6 + bp * (x - 1) - 4 * x) / 16

59 print(e)

60 B = diff(e, v, 2) * v

61 BP = -v * diff(B, v) / b

62 print(simplify(B.subs(v, v0)))

63 print(simplify(BP.subs(v, v0)))

66def pouriertarantola(V, E0, B0, BP, V0):

67 'Pourier-Tarantola equation from PRB 70, 224107'

69 eta = (V / V0)**(1 / 3)

70 squiggle = -3 * np.log(eta)

72 E = E0 + B0 * V0 * squiggle**2 / 6 * (3 + squiggle * (BP - 2))

73 return E

76def vinet(V, E0, B0, BP, V0):

77 'Vinet equation from PRB 70, 224107'

79 eta = (V / V0)**(1 / 3)

81 E = (E0 + 2 * B0 * V0 / (BP - 1)**2 *

82 (2 - (5 + 3 * BP * (eta - 1) - 3 * eta) *

83 np.exp(-3 * (BP - 1) * (eta - 1) / 2)))

84 return E

87def antonschmidt(V, Einf, B, n, V0):

88 """From Intermetallics 11, 23-32 (2003)

90 Einf should be E_infinity, i.e. infinite separation, but

91 according to the paper it does not provide a good estimate

92 of the cohesive energy. They derive this equation from an

93 empirical formula for the volume dependence of pressure,

95 E(vol) = E_inf + int(P dV) from V=vol to V=infinity

97 but the equation breaks down at large volumes, so E_inf

98 is not that meaningful

100 n should be about -2 according to the paper.

102 I find this equation does not fit volumetric data as well

103 as the other equtions do.

104 """

106 E = B * V0 / (n + 1) * (V / V0)**(n + 1) * (np.log(V / V0) -

107 (1 / (n + 1))) + Einf

108 return E

111def p3(V, c0, c1, c2, c3):

112 'polynomial fit'

114 E = c0 + c1 * V + c2 * V**2 + c3 * V**3

115 return E

118def parabola(x, a, b, c):

119 """parabola polynomial function

121 this function is used to fit the data to get good guesses for

122 the equation of state fits

124 a 4th order polynomial fit to get good guesses for

125 was not a good idea because for noisy data the fit is too wiggly

126 2nd order seems to be sufficient, and guarantees a single minimum"""

128 return a + b * x + c * x**2

131class EquationOfState:

132 """Fit equation of state for bulk systems.

134 The following equation is used::

136 sjeos (default)

137 A third order inverse polynomial fit 10.1103/PhysRevB.67.026103

139 ::

141 2 3 -1/3

142 E(V) = c + c t + c t + c t , t = V

143 0 1 2 3

145 taylor

146 A third order Taylor series expansion about the minimum volume

148 murnaghan

149 PRB 28, 5480 (1983)

151 birch

152 Intermetallic compounds: Principles and Practice,

153 Vol I: Principles. pages 195-210

155 birchmurnaghan

156 PRB 70, 224107

158 pouriertarantola

159 PRB 70, 224107

161 vinet

162 PRB 70, 224107

164 antonschmidt

165 Intermetallics 11, 23-32 (2003)

167 p3

168 A third order polynomial fit

170 Use::

172 eos = EquationOfState(volumes, energies, eos='murnaghan')

173 v0, e0, B = eos.fit()

174 eos.plot(show=True)

176 """

178 def __init__(self, volumes, energies, eos='sj'):

179 self.v = np.array(volumes)

180 self.e = np.array(energies)

182 if eos == 'sjeos':

183 eos = 'sj'

184 self.eos_string = eos

185 self.v0 = None

187 def fit(self, warn=True):

188 """Calculate volume, energy, and bulk modulus.

190 Returns the optimal volume, the minimum energy, and the bulk

191 modulus. Notice that the ASE units for the bulk modulus is

192 eV/Angstrom^3 - to get the value in GPa, do this::

194 v0, e0, B = eos.fit()

195 print(B / kJ * 1.0e24, 'GPa')

197 """

198 from scipy.optimize import curve_fit

200 if self.eos_string == 'sj':

201 return self.fit_sjeos()

203 self.func = globals()[self.eos_string]

205 p0 = [min(self.e), 1, 1]

206 popt, pcov = curve_fit(parabola, self.v, self.e, p0)

208 parabola_parameters = popt

209 # Here I just make sure the minimum is bracketed by the volumes

210 # this if for the solver

211 minvol = min(self.v)

212 maxvol = max(self.v)

214 # the minimum of the parabola is at dE/dV = 0, or 2 * c V +b =0

215 c = parabola_parameters[2]

216 b = parabola_parameters[1]

217 a = parabola_parameters[0]

218 parabola_vmin = -b / 2 / c

220 # evaluate the parabola at the minimum to estimate the groundstate

221 # energy

222 E0 = parabola(parabola_vmin, a, b, c)

223 # estimate the bulk modulus from Vo * E''. E'' = 2 * c

224 B0 = 2 * c * parabola_vmin

226 if self.eos_string == 'antonschmidt':

227 BP = -2

228 else:

229 BP = 4

231 initial_guess = [E0, B0, BP, parabola_vmin]

233 # now fit the equation of state

234 p0 = initial_guess

235 popt, pcov = curve_fit(self.func, self.v, self.e, p0)

237 self.eos_parameters = popt

239 if self.eos_string == 'p3':

240 c0, c1, c2, c3 = self.eos_parameters

241 # find minimum E in E = c0 + c1 * V + c2 * V**2 + c3 * V**3

242 # dE/dV = c1+ 2 * c2 * V + 3 * c3 * V**2 = 0

243 # solve by quadratic formula with the positive root

245 a = 3 * c3

246 b = 2 * c2

247 c = c1

249 self.v0 = (-b + np.sqrt(b**2 - 4 * a * c)) / (2 * a)

250 self.e0 = p3(self.v0, c0, c1, c2, c3)

251 self.B = (2 * c2 + 6 * c3 * self.v0) * self.v0

252 else:

253 self.v0 = self.eos_parameters[3]

254 self.e0 = self.eos_parameters[0]

255 self.B = self.eos_parameters[1]

257 if warn and not (minvol < self.v0 < maxvol):

258 warnings.warn(

259 'The minimum volume of your fit is not in '

260 'your volumes. You may not have a minimum in your dataset!')

262 return self.v0, self.e0, self.B

264 def getplotdata(self):

265 if self.v0 is None:

266 self.fit()

268 x = np.linspace(min(self.v), max(self.v), 100)

269 if self.eos_string == 'sj':

270 y = self.fit0(x**-(1 / 3))

271 else:

272 y = self.func(x, *self.eos_parameters)

274 return self.eos_string, self.e0, self.v0, self.B, x, y, self.v, self.e

276 def plot(self, filename=None, show=False, ax=None):

277 """Plot fitted energy curve.

279 Uses Matplotlib to plot the energy curve. Use *show=True* to

280 show the figure and *filename='abc.png'* or

281 *filename='abc.eps'* to save the figure to a file."""

283 import matplotlib.pyplot as plt

285 plotdata = self.getplotdata()

287 ax = plot(*plotdata, ax=ax)

289 if show:

290 plt.show()

291 if filename is not None:

292 fig = ax.get_figure()

293 fig.savefig(filename)

294 return ax

296 def fit_sjeos(self):

297 """Calculate volume, energy, and bulk modulus.

299 Returns the optimal volume, the minimum energy, and the bulk

300 modulus. Notice that the ASE units for the bulk modulus is

301 eV/Angstrom^3 - to get the value in GPa, do this::

303 v0, e0, B = eos.fit()

304 print(B / kJ * 1.0e24, 'GPa')

306 """

308 fit0 = np.poly1d(np.polyfit(self.v**-(1 / 3), self.e, 3))

309 fit1 = np.polyder(fit0, 1)

310 fit2 = np.polyder(fit1, 1)

312 self.v0 = None

313 for t in np.roots(fit1):

314 if isinstance(t, float) and t > 0 and fit2(t) > 0:

315 self.v0 = t**-3

316 break

318 if self.v0 is None:

319 raise ValueError('No minimum!')

321 self.e0 = fit0(t)

322 self.B = t**5 * fit2(t) / 9

323 self.fit0 = fit0

325 return self.v0, self.e0, self.B

328def plot(eos_string, e0, v0, B, x, y, v, e, ax=None):

329 if ax is None:

330 import matplotlib.pyplot as plt

331 ax = plt.gca()

333 ax.plot(x, y, ls='-', color='C3') # By default red line

334 ax.plot(v, e, ls='', marker='o', mec='C0',

335 mfc='C0') # By default blue marker

337 try:

338 ax.set_xlabel(u'volume [Å\$^3\$]')

339 ax.set_ylabel(u'energy [eV]')

340 ax.set_title(u'%s: E: %.3f eV, V: %.3f Å\$^3\$, B: %.3f GPa' %

341 (eos_string, e0, v0,

342 B / kJ * 1.e24))

344 except ImportError: # XXX what would cause this error? LaTeX?

345 import warnings

346 warnings.warn('Could not use LaTeX formatting')

347 ax.set_xlabel(u'volume [L(length)^3]')

348 ax.set_ylabel(u'energy [E(energy)]')

349 ax.set_title(u'%s: E: %.3f E, V: %.3f L^3, B: %.3e E/L^3' %

350 (eos_string, e0, v0, B))

352 return ax

355def calculate_eos(atoms, npoints=5, eps=0.04, trajectory=None, callback=None):

356 """Calculate equation-of-state.

358 atoms: Atoms object

359 System to calculate EOS for. Must have a calculator attached.

360 npoints: int

361 Number of points.

362 eps: float

363 Variation in volume from v0*(1-eps) to v0*(1+eps).

364 trajectory: Trjectory object or str

365 Write configurations to a trajectory file.

366 callback: function

367 Called after every energy calculation.

369 >>> from ase.build import bulk

370 >>> from ase.calculators.emt import EMT

371 >>> a = bulk('Cu', 'fcc', a=3.6)

372 >>> a.calc = EMT()

373 >>> eos = calculate_eos(a, trajectory='Cu.traj')

374 >>> v, e, B = eos.fit()

375 >>> a = (4 * v)**(1 / 3.0)

376 >>> print('{0:.6f}'.format(a))

377 3.589825

378 """

380 # Save original positions and cell:

381 p0 = atoms.get_positions()

382 c0 = atoms.get_cell()

384 if isinstance(trajectory, str):

385 from ase.io import Trajectory

386 trajectory = Trajectory(trajectory, 'w', atoms)

388 if trajectory is not None:

389 trajectory.set_description({'type': 'eos',

390 'npoints': npoints,

391 'eps': eps})

393 try:

394 energies = []

395 volumes = []

396 for x in np.linspace(1 - eps, 1 + eps, npoints)**(1 / 3):

397 atoms.set_cell(x * c0, scale_atoms=True)

398 volumes.append(atoms.get_volume())

399 energies.append(atoms.get_potential_energy())

400 if callback:

401 callback()

402 if trajectory is not None:

403 trajectory.write()

404 return EquationOfState(volumes, energies)

405 finally:

406 atoms.cell = c0

407 atoms.positions = p0

408 if trajectory is not None:

409 trajectory.close()

412class CLICommand:

413 """Calculate EOS from one or more trajectory files.

415 See https://wiki.fysik.dtu.dk/ase/tutorials/eos/eos.html for

417 """

419 @staticmethod

423 help='Plot EOS fit. Default behaviour is '

424 'to write results of fit.')

426 help='Type of fit. Must be one of {}.'

427 .format(', '.join(eos_names)))

429 @staticmethod

430 def run(args):

433 if not args.plot:

434 print('# filename '

435 'points volume energy bulk modulus')

436 print('# '

437 ' [Ang^3] [eV] [GPa]')

438 for name in args.trajectories:

439 if name == '-':

440 # Special case - used by ASE's GUI:

441 import pickle

442 import sys

444 else:

445 if '@' in name:

446 index = None

447 else:

448 index = ':'

450 v = [atoms.get_volume() for atoms in images]

451 e = [atoms.get_potential_energy() for atoms in images]

452 eos = EquationOfState(v, e, args.type)

453 if args.plot:

454 eos.plot()

455 else:

456 try:

457 v0, e0, B = eos.fit()

458 except ValueError as ex:

459 print('{:30}{:2} {}'

460 .format(name, len(v), ex.message))

461 else:

462 print('{:30}{:2} {:10.3f}{:10.3f}{:14.3f}'

463 .format(name, len(v), v0, e0, B / kJ * 1.0e24))