Coverage for /builds/ase/ase/ase/eos.py: 63.92%

1# fmt: off

3import warnings

5import numpy as np

7from ase.units import kJ

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

10 'pouriertarantola', 'vinet', 'antonschmidt', 'p3']

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

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

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

17 return E

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

21 'From PRB 28,5480 (1983'

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

24 return E

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

28 """

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

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

31 paper downloaded from Web

33 case where n=0

34 """

36 E = (E0 +

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

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

39 return E

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

43 """

44 BirchMurnaghan equation from PRB 70, 224107

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

46 inversed expression for eta.

47 """

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

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

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

52 return E

55def check_birchmurnaghan():

56 from sympy import Rational, diff, simplify, symbols

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

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

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

60 print(e)

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

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

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

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

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

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

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

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

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

74 return E

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

78 'Vinet equation from PRB 70, 224107'

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

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

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

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

85 return E

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

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

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

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

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

94 empirical formula for the volume dependence of pressure,

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

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

99 is not that meaningful

100

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

102

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

104 as the other equtions do.

105 """

106

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

108 (1 / (n + 1))) + Einf

109 return E

110

111

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

113 'polynomial fit'

114

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

116 return E

117

118

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

120 """parabola polynomial function

121

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

123 the equation of state fits

124

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

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

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

128

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

130

131

132class EquationOfState:

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

134

135 The following equation is used::

136

137 sjeos (default)

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

139

140 ::

141

142 2 3 -1/3

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

144 0 1 2 3

145

146 taylor

147 A third order Taylor series expansion about the minimum volume

148

149 murnaghan

150 PRB 28, 5480 (1983)

151

152 birch

153 Intermetallic compounds: Principles and Practice,

154 Vol I: Principles. pages 195-210

155

156 birchmurnaghan

157 PRB 70, 224107

158

159 pouriertarantola

160 PRB 70, 224107

161

162 vinet

163 PRB 70, 224107

164

165 antonschmidt

166 Intermetallics 11, 23-32 (2003)

167

168 p3

169 A third order polynomial fit

170

171 Use::

172

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

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

175 eos.plot(show=True)

176

177 """

178

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

180 self.v = np.array(volumes)

181 self.e = np.array(energies)

182

183 if eos == 'sjeos':

184 eos = 'sj'

185 self.eos_string = eos

186 self.v0 = None

187

188 def fit(self, warn=True):

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

190

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

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

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

194

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

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

197

198 """

199 from scipy.optimize import curve_fit

200

201 if self.eos_string == 'sj':

202 return self.fit_sjeos()

203

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

205

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

207 popt, _pcov = curve_fit(parabola, self.v, self.e, p0)

208

209 parabola_parameters = popt

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

211 # this if for the solver

212 minvol = min(self.v)

213 maxvol = max(self.v)

214

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

216 c = parabola_parameters[2]

217 b = parabola_parameters[1]

218 a = parabola_parameters[0]

219 parabola_vmin = -b / 2 / c

220

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

222 # energy

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

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

225 B0 = 2 * c * parabola_vmin

226

227 if self.eos_string == 'antonschmidt':

228 BP = -2

229 else:

230 BP = 4

231

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

233

234 # now fit the equation of state

235 p0 = initial_guess

236 popt, _pcov = curve_fit(self.func, self.v, self.e, p0)

237

238 self.eos_parameters = popt

239

240 if self.eos_string == 'p3':

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

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

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

244 # solve by quadratic formula with the positive root

245

246 a = 3 * c3

247 b = 2 * c2

248 c = c1

249

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

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

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

253 else:

254 self.v0 = self.eos_parameters[3]

255 self.e0 = self.eos_parameters[0]

256 self.B = self.eos_parameters[1]

257

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

259 warnings.warn(

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

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

262

263 return self.v0, self.e0, self.B

264

265 def getplotdata(self):

266 if self.v0 is None:

267 self.fit()

268

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

270 if self.eos_string == 'sj':

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

272 else:

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

274

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

276

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

278 """Plot fitted energy curve.

279

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

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

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

283

284 import matplotlib.pyplot as plt

285

286 plotdata = self.getplotdata()

287

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

289

290 if show:

291 plt.show()

292 if filename is not None:

293 fig = ax.get_figure()

294 fig.savefig(filename)

295 return ax

296

297 def fit_sjeos(self):

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

299

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

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

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

303

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

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

306

307 """

308

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

310 fit1 = np.polyder(fit0, 1)

311 fit2 = np.polyder(fit1, 1)

312

313 self.v0 = None

314 for t in np.roots(fit1):

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

316 self.v0 = t**-3

317 break

318

319 if self.v0 is None:

320 raise ValueError('No minimum!')

321

322 self.e0 = fit0(t)

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

324 self.fit0 = fit0

325

326 return self.v0, self.e0, self.B

327

328

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

330 if ax is None:

331 import matplotlib.pyplot as plt

332 ax = plt.gca()

333

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

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

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

337

338 try:

339 ax.set_xlabel('volume [Å$^3$]')

340 ax.set_ylabel('energy [eV]')

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

342 (eos_string, e0, v0,

343 B / kJ * 1.e24))

344

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

346 import warnings

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

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

349 ax.set_ylabel('energy [E(energy)]')

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

351 (eos_string, e0, v0, B))

352

353 return ax

354

355

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

357 """Calculate equation-of-state.

358

359 atoms: Atoms object

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

361 npoints: int

362 Number of points.

363 eps: float

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

365 trajectory: Trjectory object or str

366 Write configurations to a trajectory file.

367 callback: function

368 Called after every energy calculation.

369

370 >>> from ase.build import bulk

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

372 >>> from ase.eos import calculate_eos

373

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

375 >>> a.calc = EMT()

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

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

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

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

380 3.589825

381 """

382

383 # Save original positions and cell:

384 p0 = atoms.get_positions()

385 c0 = atoms.get_cell()

386

387 if isinstance(trajectory, str):

388 from ase.io import Trajectory

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

390

391 if trajectory is not None:

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

393 'npoints': npoints,

394 'eps': eps})

395

396 try:

397 energies = []

398 volumes = []

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

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

401 volumes.append(atoms.get_volume())

402 energies.append(atoms.get_potential_energy())

403 if callback:

404 callback()

405 if trajectory is not None:

406 trajectory.write()

407 return EquationOfState(volumes, energies)

408 finally:

409 atoms.cell = c0

410 atoms.positions = p0

411 if trajectory is not None:

412 trajectory.close()

413

414

415class CLICommand:

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

417

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

419 more information.

420 """

421

422 @staticmethod

423 def add_arguments(parser):

424 parser.add_argument('trajectories', nargs='+', metavar='trajectory')

425 parser.add_argument('-p', '--plot', action='store_true',

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

427 'to write results of fit.')

428 parser.add_argument('-t', '--type', default='sj',

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

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

431

432 @staticmethod

433 def run(args):

434 from ase.io import read

435

436 if not args.plot:

437 print('# filename '

438 'points volume energy bulk modulus')

439 print('# '

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

441 for name in args.trajectories:

442 if name == '-':

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

444 import pickle

445 import sys

446 v, e = pickle.load(sys.stdin.buffer)

447 else:

448 if '@' in name:

449 index = None

450 else:

451 index = ':'

452 images = read(name, index=index)

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

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

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

456 if args.plot:

457 eos.plot()

458 else:

459 try:

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

461 except ValueError as ex:

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

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

464 else:

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

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