1import warnings 

2 

3from ase.units import kJ 

4 

5import numpy as np 

6 

7 

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

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

10 

11 

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

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

14 

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

16 return E 

17 

18 

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

20 'From PRB 28,5480 (1983' 

21 

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

23 return E 

24 

25 

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 

30 paper downloaded from Web 

31 

32 case where n=0 

33 """ 

34 

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 

39 

40 

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

47 

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 

52 

53 

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

64 

65 

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

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

68 

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

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

71 

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

73 return E 

74 

75 

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

77 'Vinet equation from PRB 70, 224107' 

78 

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

80 

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 

85 

86 

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

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

89 

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, 

94 

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

96 

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

98 is not that meaningful 

99 

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

101 

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

103 as the other equtions do. 

104 """ 

105 

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

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

108 return E 

109 

110 

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

112 'polynomial fit' 

113 

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

115 return E 

116 

117 

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

119 """parabola polynomial function 

120 

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

122 the equation of state fits 

123 

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

127 

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

129 

130 

131class EquationOfState: 

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

133 

134 The following equation is used:: 

135 

136 sjeos (default) 

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

138 

139 :: 

140 

141 2 3 -1/3 

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

143 0 1 2 3 

144 

145 taylor 

146 A third order Taylor series expansion about the minimum volume 

147 

148 murnaghan 

149 PRB 28, 5480 (1983) 

150 

151 birch 

152 Intermetallic compounds: Principles and Practice, 

153 Vol I: Principles. pages 195-210 

154 

155 birchmurnaghan 

156 PRB 70, 224107 

157 

158 pouriertarantola 

159 PRB 70, 224107 

160 

161 vinet 

162 PRB 70, 224107 

163 

164 antonschmidt 

165 Intermetallics 11, 23-32 (2003) 

166 

167 p3 

168 A third order polynomial fit 

169 

170 Use:: 

171 

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

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

174 eos.plot(show=True) 

175 

176 """ 

177 

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

179 self.v = np.array(volumes) 

180 self.e = np.array(energies) 

181 

182 if eos == 'sjeos': 

183 eos = 'sj' 

184 self.eos_string = eos 

185 self.v0 = None 

186 

187 def fit(self, warn=True): 

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

189 

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

193 

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

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

196 

197 """ 

198 from scipy.optimize import curve_fit 

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', 

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

336 

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

343 

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

351 

352 return ax 

353 

354 

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

356 """Calculate equation-of-state. 

357 

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. 

368 

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

379 

380 # Save original positions and cell: 

381 p0 = atoms.get_positions() 

382 c0 = atoms.get_cell() 

383 

384 if isinstance(trajectory, str): 

385 from ase.io import Trajectory 

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

387 

388 if trajectory is not None: 

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

390 'npoints': npoints, 

391 'eps': eps}) 

392 

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

410 

411 

412class CLICommand: 

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

414 

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

416 more information. 

417 """ 

418 

419 @staticmethod 

420 def add_arguments(parser): 

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

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

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

424 'to write results of fit.') 

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

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

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

428 

429 @staticmethod 

430 def run(args): 

431 from ase.io import read 

432 

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 

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

444 else: 

445 if '@' in name: 

446 index = None 

447 else: 

448 index = ':' 

449 images = read(name, index=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))