1import warnings 

2 

3from ase.units import kJ 

4 

5import numpy as np 

6 

7from scipy.optimize import curve_fit 

8 

9 

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

11 'pouriertarantola', 'vinet', 'antonschmidt', 'p3'] 

12 

13 

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

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

16 

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

18 return E 

19 

20 

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

22 'From PRB 28,5480 (1983' 

23 

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

25 return E 

26 

27 

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

29 """ 

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

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

32 paper downloaded from Web 

33 

34 case where n=0 

35 """ 

36 

37 E = (E0 + 

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

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

40 return E 

41 

42 

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

44 """ 

45 BirchMurnaghan equation from PRB 70, 224107 

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

47 inversed expression for eta. 

48 """ 

49 

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

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

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

53 return E 

54 

55 

56def check_birchmurnaghan(): 

57 from sympy import symbols, Rational, diff, simplify 

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

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

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

61 print(e) 

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

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

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

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

66 

67 

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

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

70 

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

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

73 

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

75 return E 

76 

77 

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

79 'Vinet equation from PRB 70, 224107' 

80 

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

82 

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

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

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

86 return E 

87 

88 

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

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

91 

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

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

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

95 empirical formula for the volume dependence of pressure, 

96 

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

98 

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

100 is not that meaningful 

101 

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

103 

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

105 as the other equtions do. 

106 """ 

107 

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

109 (1 / (n + 1))) + Einf 

110 return E 

111 

112 

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

114 'polynomial fit' 

115 

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

117 return E 

118 

119 

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

121 """parabola polynomial function 

122 

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

124 the equation of state fits 

125 

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

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

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

129 

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

131 

132 

133class EquationOfState: 

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

135 

136 The following equation is used:: 

137 

138 sjeos (default) 

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

140 

141 :: 

142 

143 2 3 -1/3 

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

145 0 1 2 3 

146 

147 taylor 

148 A third order Taylor series expansion about the minimum volume 

149 

150 murnaghan 

151 PRB 28, 5480 (1983) 

152 

153 birch 

154 Intermetallic compounds: Principles and Practice, 

155 Vol I: Principles. pages 195-210 

156 

157 birchmurnaghan 

158 PRB 70, 224107 

159 

160 pouriertarantola 

161 PRB 70, 224107 

162 

163 vinet 

164 PRB 70, 224107 

165 

166 antonschmidt 

167 Intermetallics 11, 23-32 (2003) 

168 

169 p3 

170 A third order polynomial fit 

171 

172 Use:: 

173 

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

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

176 eos.plot(show=True) 

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 

200 if self.eos_string == 'sj': 

201 return self.fit_sjeos() 

202 

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

204 

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

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

207 

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) 

213 

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 

219 

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 

225 

226 if self.eos_string == 'antonschmidt': 

227 BP = -2 

228 else: 

229 BP = 4 

230 

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

232 

233 # now fit the equation of state 

234 p0 = initial_guess 

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

236 

237 self.eos_parameters = popt 

238 

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 

244 

245 a = 3 * c3 

246 b = 2 * c2 

247 c = c1 

248 

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] 

256 

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

261 

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

263 

264 def getplotdata(self): 

265 if self.v0 is None: 

266 self.fit() 

267 

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) 

273 

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

275 

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

277 """Plot fitted energy curve. 

278 

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.""" 

282 

283 import matplotlib.pyplot as plt 

284 

285 plotdata = self.getplotdata() 

286 

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

288 

289 if show: 

290 plt.show() 

291 if filename is not None: 

292 fig = ax.get_figure() 

293 fig.savefig(filename) 

294 return ax 

295 

296 def fit_sjeos(self): 

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

298 

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

302 

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

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

305 

306 """ 

307 

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) 

311 

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 

317 

318 if self.v0 is None: 

319 raise ValueError('No minimum!') 

320 

321 self.e0 = fit0(t) 

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

323 self.fit0 = fit0 

324 

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

326 

327 

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

332 

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

334 ax.plot(v, e, ls='', marker='o', mec='C0', mfc='C0') # By default blue marker 

335 

336 try: 

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

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

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

340 (eos_string, e0, v0, 

341 B / kJ * 1.e24)) 

342 

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

344 import warnings 

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

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

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

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

349 (eos_string, e0, v0, B)) 

350 

351 return ax 

352 

353 

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

355 """Calculate equation-of-state. 

356 

357 atoms: Atoms object 

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

359 npoints: int 

360 Number of points. 

361 eps: float 

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

363 trajectory: Trjectory object or str 

364 Write configurations to a trajectory file. 

365 callback: function 

366 Called after every energy calculation. 

367 

368 >>> from ase.build import bulk 

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

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

371 >>> a.calc = EMT() 

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

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

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

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

376 3.589825 

377 """ 

378 

379 # Save original positions and cell: 

380 p0 = atoms.get_positions() 

381 c0 = atoms.get_cell() 

382 

383 if isinstance(trajectory, str): 

384 from ase.io import Trajectory 

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

386 

387 if trajectory is not None: 

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

389 'npoints': npoints, 

390 'eps': eps}) 

391 

392 try: 

393 energies = [] 

394 volumes = [] 

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

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

397 volumes.append(atoms.get_volume()) 

398 energies.append(atoms.get_potential_energy()) 

399 if callback: 

400 callback() 

401 if trajectory is not None: 

402 trajectory.write() 

403 return EquationOfState(volumes, energies) 

404 finally: 

405 atoms.cell = c0 

406 atoms.positions = p0 

407 if trajectory is not None: 

408 trajectory.close() 

409 

410 

411class CLICommand: 

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

413 

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

415 more information. 

416 """ 

417 

418 @staticmethod 

419 def add_arguments(parser): 

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

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

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

423 'to write results of fit.') 

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

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

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

427 

428 @staticmethod 

429 def run(args): 

430 from ase.io import read 

431 

432 if not args.plot: 

433 print('# filename ' 

434 'points volume energy bulk modulus') 

435 print('# ' 

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

437 for name in args.trajectories: 

438 if name == '-': 

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

440 import pickle 

441 import sys 

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

443 else: 

444 if '@' in name: 

445 index = None 

446 else: 

447 index = ':' 

448 images = read(name, index=index) 

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

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

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

452 if args.plot: 

453 eos.plot() 

454 else: 

455 try: 

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

457 except ValueError as ex: 

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

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

460 else: 

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

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