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

1import warnings

3import numpy as np

5from ase.units import kJ

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

8 'pouriertarantola', 'vinet', 'antonschmidt', 'p3']

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

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

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

15 return E

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

19 'From PRB 28,5480 (1983'

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

22 return E

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

26 """

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

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

29 paper downloaded from Web

31 case where n=0

32 """

34 E = (E0 +

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

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

37 return E

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

41 """

42 BirchMurnaghan equation from PRB 70, 224107

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

44 inversed expression for eta.

45 """

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

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

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

50 return E

53def check_birchmurnaghan():

54 from sympy import Rational, diff, simplify, symbols

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

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

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

58 print(e)

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

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

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

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

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

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

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

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

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

72 return E

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

76 'Vinet equation from PRB 70, 224107'

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

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

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

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

83 return E

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

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

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

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

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

92 empirical formula for the volume dependence of pressure,

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

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

97 is not that meaningful

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

100

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

102 as the other equtions do.

103 """

104

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

106 (1 / (n + 1))) + Einf

107 return E

108

109

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

111 'polynomial fit'

112

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

114 return E

115

116

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

118 """parabola polynomial function

119

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

121 the equation of state fits

122

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

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

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

126

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

128

129

130class EquationOfState:

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

132

133 The following equation is used::

134

135 sjeos (default)

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

137

138 ::

139

140 2 3 -1/3

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

142 0 1 2 3

143

144 taylor

145 A third order Taylor series expansion about the minimum volume

146

147 murnaghan

148 PRB 28, 5480 (1983)

149

150 birch

151 Intermetallic compounds: Principles and Practice,

152 Vol I: Principles. pages 195-210

153

154 birchmurnaghan

155 PRB 70, 224107

156

157 pouriertarantola

158 PRB 70, 224107

159

160 vinet

161 PRB 70, 224107

162

163 antonschmidt

164 Intermetallics 11, 23-32 (2003)

165

166 p3

167 A third order polynomial fit

168

169 Use::

170

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

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

173 eos.plot(show=True)

174

175 """

176

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

178 self.v = np.array(volumes)

179 self.e = np.array(energies)

180

181 if eos == 'sjeos':

182 eos = 'sj'

183 self.eos_string = eos

184 self.v0 = None

185

186 def fit(self, warn=True):

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

188

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

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

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

192

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

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

195

196 """

197 from scipy.optimize import curve_fit

198

199 if self.eos_string == 'sj':

200 return self.fit_sjeos()

201

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

203

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

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

206

207 parabola_parameters = popt

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

209 # this if for the solver

210 minvol = min(self.v)

211 maxvol = max(self.v)

212

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

214 c = parabola_parameters[2]

215 b = parabola_parameters[1]

216 a = parabola_parameters[0]

217 parabola_vmin = -b / 2 / c

218

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

220 # energy

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

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

223 B0 = 2 * c * parabola_vmin

224

225 if self.eos_string == 'antonschmidt':

226 BP = -2

227 else:

228 BP = 4

229

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

231

232 # now fit the equation of state

233 p0 = initial_guess

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

235

236 self.eos_parameters = popt

237

238 if self.eos_string == 'p3':

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

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

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

242 # solve by quadratic formula with the positive root

243

244 a = 3 * c3

245 b = 2 * c2

246 c = c1

247

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

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

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

251 else:

252 self.v0 = self.eos_parameters[3]

253 self.e0 = self.eos_parameters[0]

254 self.B = self.eos_parameters[1]

255

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

257 warnings.warn(

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

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

260

261 return self.v0, self.e0, self.B

262

263 def getplotdata(self):

264 if self.v0 is None:

265 self.fit()

266

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

268 if self.eos_string == 'sj':

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

270 else:

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

272

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

274

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

276 """Plot fitted energy curve.

277

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

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

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

281

282 import matplotlib.pyplot as plt

283

284 plotdata = self.getplotdata()

285

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

287

288 if show:

289 plt.show()

290 if filename is not None:

291 fig = ax.get_figure()

292 fig.savefig(filename)

293 return ax

294

295 def fit_sjeos(self):

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

297

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

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

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

301

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

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

304

305 """

306

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

308 fit1 = np.polyder(fit0, 1)

309 fit2 = np.polyder(fit1, 1)

310

311 self.v0 = None

312 for t in np.roots(fit1):

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

314 self.v0 = t**-3

315 break

316

317 if self.v0 is None:

318 raise ValueError('No minimum!')

319

320 self.e0 = fit0(t)

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

322 self.fit0 = fit0

323

324 return self.v0, self.e0, self.B

325

326

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

328 if ax is None:

329 import matplotlib.pyplot as plt

330 ax = plt.gca()

331

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

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

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

335

336 try:

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

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

339 ax.set_title('%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('volume [L(length)^3]')

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

348 ax.set_title('%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 >>> from ase.eos import calculate_eos

371

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

373 >>> a.calc = EMT()

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

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

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

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

378 3.589825

379 """

380

381 # Save original positions and cell:

382 p0 = atoms.get_positions()

383 c0 = atoms.get_cell()

384

385 if isinstance(trajectory, str):

386 from ase.io import Trajectory

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

388

389 if trajectory is not None:

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

391 'npoints': npoints,

392 'eps': eps})

393

394 try:

395 energies = []

396 volumes = []

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

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

399 volumes.append(atoms.get_volume())

400 energies.append(atoms.get_potential_energy())

401 if callback:

402 callback()

403 if trajectory is not None:

404 trajectory.write()

405 return EquationOfState(volumes, energies)

406 finally:

407 atoms.cell = c0

408 atoms.positions = p0

409 if trajectory is not None:

410 trajectory.close()

411

412

413class CLICommand:

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

415

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

417 more information.

418 """

419

420 @staticmethod

421 def add_arguments(parser):

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

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

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

425 'to write results of fit.')

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

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

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

429

430 @staticmethod

431 def run(args):

432 from ase.io import read

433

434 if not args.plot:

435 print('# filename '

436 'points volume energy bulk modulus')

437 print('# '

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

439 for name in args.trajectories:

440 if name == '-':

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

442 import pickle

443 import sys

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

445 else:

446 if '@' in name:

447 index = None

448 else:

449 index = ':'

450 images = read(name, index=index)

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

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

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

454 if args.plot:

455 eos.plot()

456 else:

457 try:

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

459 except ValueError as ex:

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

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

462 else:

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

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