r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Modules for calculating thermochemical information from computational

2outputs."""

4import os

5import sys

6import numpy as np

8from ase import units

11class ThermoChem:

12 """Base class containing common methods used in thermochemistry

13 calculations."""

15 def get_ZPE_correction(self):

16 """Returns the zero-point vibrational energy correction in eV."""

17 zpe = 0.

18 for energy in self.vib_energies:

19 zpe += 0.5 * energy

20 return zpe

22 def _vibrational_energy_contribution(self, temperature):

23 """Calculates the change in internal energy due to vibrations from

24 0K to the specified temperature for a set of vibrations given in

25 eV and a temperature given in Kelvin. Returns the energy change

26 in eV."""

27 kT = units.kB * temperature

28 dU = 0.

29 for energy in self.vib_energies:

30 dU += energy / (np.exp(energy / kT) - 1.)

31 return dU

33 def _vibrational_entropy_contribution(self, temperature):

34 """Calculates the entropy due to vibrations for a set of vibrations

35 given in eV and a temperature given in Kelvin. Returns the entropy

36 in eV/K."""

37 kT = units.kB * temperature

38 S_v = 0.

39 for energy in self.vib_energies:

40 x = energy / kT

41 S_v += x / (np.exp(x) - 1.) - np.log(1. - np.exp(-x))

42 S_v *= units.kB

43 return S_v

45 def _vprint(self, text):

46 """Print output if verbose flag True."""

47 if self.verbose:

48 sys.stdout.write(text + os.linesep)

51class HarmonicThermo(ThermoChem):

52 """Class for calculating thermodynamic properties in the approximation

53 that all degrees of freedom are treated harmonically. Often used for

56 Inputs:

58 vib_energies : list

59 a list of the harmonic energies of the adsorbate (e.g., from

60 ase.vibrations.Vibrations.get_energies). The number of

61 energies should match the number of degrees of freedom of the

62 adsorbate; i.e., 3*n, where n is the number of atoms. Note that

63 this class does not check that the user has supplied the correct

64 number of energies. Units of energies are eV.

65 potentialenergy : float

66 the potential energy in eV (e.g., from atoms.get_potential_energy)

67 (if potentialenergy is unspecified, then the methods of this

68 class can be interpreted as the energy corrections)

69 """

71 def __init__(self, vib_energies, potentialenergy=0.):

72 self.vib_energies = vib_energies

73 # Check for imaginary frequencies.

74 if sum(np.iscomplex(self.vib_energies)):

75 raise ValueError('Imaginary vibrational energies are present.')

76 else:

77 self.vib_energies = np.real(self.vib_energies) # clear +0.j

79 self.potentialenergy = potentialenergy

81 def get_internal_energy(self, temperature, verbose=True):

82 """Returns the internal energy, in eV, in the harmonic approximation

83 at a specified temperature (K)."""

85 self.verbose = verbose

86 write = self._vprint

87 fmt = '%-15s%13.3f eV'

88 write('Internal energy components at T = %.2f K:' % temperature)

89 write('=' * 31)

91 U = 0.

93 write(fmt % ('E_pot', self.potentialenergy))

94 U += self.potentialenergy

96 zpe = self.get_ZPE_correction()

97 write(fmt % ('E_ZPE', zpe))

98 U += zpe

100 dU_v = self._vibrational_energy_contribution(temperature)

101 write(fmt % ('Cv_harm (0->T)', dU_v))

102 U += dU_v

104 write('-' * 31)

105 write(fmt % ('U', U))

106 write('=' * 31)

107 return U

109 def get_entropy(self, temperature, verbose=True):

110 """Returns the entropy, in eV/K, in the harmonic approximation

111 at a specified temperature (K)."""

113 self.verbose = verbose

114 write = self._vprint

115 fmt = '%-15s%13.7f eV/K%13.3f eV'

116 write('Entropy components at T = %.2f K:' % temperature)

117 write('=' * 49)

118 write('%15s%13s %13s' % ('', 'S', 'T*S'))

120 S = 0.

122 S_v = self._vibrational_entropy_contribution(temperature)

123 write(fmt % ('S_harm', S_v, S_v * temperature))

124 S += S_v

126 write('-' * 49)

127 write(fmt % ('S', S, S * temperature))

128 write('=' * 49)

129 return S

131 def get_helmholtz_energy(self, temperature, verbose=True):

132 """Returns the Helmholtz free energy, in eV, in the harmonic

133 approximation at a specified temperature (K)."""

135 self.verbose = True

136 write = self._vprint

138 U = self.get_internal_energy(temperature, verbose=verbose)

139 write('')

140 S = self.get_entropy(temperature, verbose=verbose)

141 F = U - temperature * S

143 write('')

144 write('Free energy components at T = %.2f K:' % temperature)

145 write('=' * 23)

146 fmt = '%5s%15.3f eV'

147 write(fmt % ('U', U))

148 write(fmt % ('-T*S', -temperature * S))

149 write('-' * 23)

150 write(fmt % ('F', F))

151 write('=' * 23)

152 return F

155class HinderedThermo(ThermoChem):

156 """Class for calculating thermodynamic properties in the hindered

157 translator and hindered rotor model where all but three degrees of

158 freedom are treated as harmonic vibrations, two are treated as

159 hindered translations, and one is treated as a hindered rotation.

161 Inputs:

163 vib_energies : list

164 a list of all the vibrational energies of the adsorbate (e.g., from

165 ase.vibrations.Vibrations.get_energies). The number of energies

166 should match the number of degrees of freedom of the adsorbate;

167 i.e., 3*n, where n is the number of atoms. Note that this class does

168 not check that the user has supplied the correct number of energies.

169 Units of energies are eV.

170 trans_barrier_energy : float

171 the translational energy barrier in eV. This is the barrier for an

172 adsorbate to diffuse on the surface.

173 rot_barrier_energy : float

174 the rotational energy barrier in eV. This is the barrier for an

176 sitedensity : float

177 density of surface sites in cm^-2

178 rotationalminima : integer

179 the number of equivalent minima for an adsorbate's full rotation.

180 For example, 6 for an adsorbate on an fcc(111) top site

181 potentialenergy : float

182 the potential energy in eV (e.g., from atoms.get_potential_energy)

183 (if potentialenergy is unspecified, then the methods of this class

184 can be interpreted as the energy corrections)

185 mass : float

186 the mass of the adsorbate in amu (if mass is unspecified, then it will

187 be calculated from the atoms class)

188 inertia : float

189 the reduced moment of inertia of the adsorbate in amu*Ang^-2

190 (if inertia is unspecified, then it will be calculated from the

191 atoms class)

192 atoms : an ASE atoms object

193 used to calculate rotational moments of inertia and molecular mass

194 symmetrynumber : integer

195 symmetry number of the adsorbate. This is the number of symmetric arms

196 of the adsorbate and depends upon how it is bound to the surface.

197 For example, propane bound through its end carbon has a symmetry

198 number of 1 but propane bound through its middle carbon has a symmetry

199 number of 2. (if symmetrynumber is unspecified, then the default is 1)

200 """

202 def __init__(self, vib_energies, trans_barrier_energy, rot_barrier_energy,

203 sitedensity, rotationalminima, potentialenergy=0.,

204 mass=None, inertia=None, atoms=None, symmetrynumber=1):

205 self.vib_energies = sorted(vib_energies, reverse=True)[:-3]

206 self.trans_barrier_energy = trans_barrier_energy * units._e

207 self.rot_barrier_energy = rot_barrier_energy * units._e

208 self.area = 1. / sitedensity / 100.0**2

209 self.rotationalminima = rotationalminima

210 self.potentialenergy = potentialenergy

211 self.atoms = atoms

212 self.symmetry = symmetrynumber

214 if (mass or atoms) and (inertia or atoms):

215 if mass:

216 self.mass = mass * units._amu

217 elif atoms:

218 self.mass = np.sum(atoms.get_masses()) * units._amu

219 if inertia:

220 self.inertia = inertia * units._amu / units.m**2

221 elif atoms:

222 self.inertia = (atoms.get_moments_of_inertia()[2] *

223 units._amu / units.m**2)

224 else:

225 raise RuntimeError('Either mass and inertia of the '

226 'adsorbate must be specified or '

227 'atoms must be specified.')

229 # Make sure no imaginary frequencies remain.

230 if sum(np.iscomplex(self.vib_energies)):

231 raise ValueError('Imaginary frequencies are present.')

232 else:

233 self.vib_energies = np.real(self.vib_energies) # clear +0.j

235 # Calculate hindered translational and rotational frequencies

236 self.freq_t = np.sqrt(self.trans_barrier_energy / (2 * self.mass *

237 self.area))

238 self.freq_r = 1. / (2 * np.pi) * np.sqrt(self.rotationalminima**2 *

239 self.rot_barrier_energy /

240 (2 * self.inertia))

242 def get_internal_energy(self, temperature, verbose=True):

243 """Returns the internal energy (including the zero point energy),

244 in eV, in the hindered translator and hindered rotor model at a

245 specified temperature (K)."""

247 from scipy.special import iv

249 self.verbose = verbose

250 write = self._vprint

251 fmt = '%-15s%13.3f eV'

252 write('Internal energy components at T = %.2f K:' % temperature)

253 write('=' * 31)

255 U = 0.

257 write(fmt % ('E_pot', self.potentialenergy))

258 U += self.potentialenergy

260 # Translational Energy

261 T_t = units._k * temperature / (units._hplanck * self.freq_t)

262 R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)

263 dU_t = 2 * (-1. / 2 - 1. / T_t / (2 + 16 * R_t) + R_t / 2 / T_t -

264 R_t / 2 / T_t *

265 iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +

266 1. / T_t / (np.exp(1. / T_t) - 1))

267 dU_t *= units.kB * temperature

268 write(fmt % ('E_trans', dU_t))

269 U += dU_t

271 # Rotational Energy

272 T_r = units._k * temperature / (units._hplanck * self.freq_r)

273 R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)

274 dU_r = (-1. / 2 - 1. / T_r / (2 + 16 * R_r) + R_r / 2 / T_r -

275 R_r / 2 / T_r *

276 iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +

277 1. / T_r / (np.exp(1. / T_r) - 1))

278 dU_r *= units.kB * temperature

279 write(fmt % ('E_rot', dU_r))

280 U += dU_r

282 # Vibrational Energy

283 dU_v = self._vibrational_energy_contribution(temperature)

284 write(fmt % ('E_vib', dU_v))

285 U += dU_v

287 # Zero Point Energy

288 dU_zpe = self.get_zero_point_energy()

289 write(fmt % ('E_ZPE', dU_zpe))

290 U += dU_zpe

292 write('-' * 31)

293 write(fmt % ('U', U))

294 write('=' * 31)

295 return U

297 def get_zero_point_energy(self, verbose=True):

298 """Returns the zero point energy, in eV, in the hindered

299 translator and hindered rotor model"""

301 zpe_t = 2 * (1. / 2 * self.freq_t * units._hplanck / units._e)

302 zpe_r = 1. / 2 * self.freq_r * units._hplanck / units._e

303 zpe_v = self.get_ZPE_correction()

304 zpe = zpe_t + zpe_r + zpe_v

305 return zpe

307 def get_entropy(self, temperature, verbose=True):

308 """Returns the entropy, in eV/K, in the hindered translator

309 and hindered rotor model at a specified temperature (K)."""

311 from scipy.special import iv

313 self.verbose = verbose

314 write = self._vprint

315 fmt = '%-15s%13.7f eV/K%13.3f eV'

316 write('Entropy components at T = %.2f K:' % temperature)

317 write('=' * 49)

318 write('%15s%13s %13s' % ('', 'S', 'T*S'))

320 S = 0.

322 # Translational Entropy

323 T_t = units._k * temperature / (units._hplanck * self.freq_t)

324 R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)

325 S_t = 2 * (-1. / 2 + 1. / 2 * np.log(np.pi * R_t / T_t) -

326 R_t / 2 / T_t *

327 iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +

328 np.log(iv(0, R_t / 2 / T_t)) +

329 1. / T_t / (np.exp(1. / T_t) - 1) -

330 np.log(1 - np.exp(-1. / T_t)))

331 S_t *= units.kB

332 write(fmt % ('S_trans', S_t, S_t * temperature))

333 S += S_t

335 # Rotational Entropy

336 T_r = units._k * temperature / (units._hplanck * self.freq_r)

337 R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)

338 S_r = (-1. / 2 + 1. / 2 * np.log(np.pi * R_r / T_r) -

339 np.log(self.symmetry) -

340 R_r / 2 / T_r * iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +

341 np.log(iv(0, R_r / 2 / T_r)) +

342 1. / T_r / (np.exp(1. / T_r) - 1) -

343 np.log(1 - np.exp(-1. / T_r)))

344 S_r *= units.kB

345 write(fmt % ('S_rot', S_r, S_r * temperature))

346 S += S_r

348 # Vibrational Entropy

349 S_v = self._vibrational_entropy_contribution(temperature)

350 write(fmt % ('S_vib', S_v, S_v * temperature))

351 S += S_v

353 # Concentration Related Entropy

354 N_over_A = np.exp(1. / 3) * (10.0**5 /

355 (units._k * temperature))**(2. / 3)

356 S_c = 1 - np.log(N_over_A) - np.log(self.area)

357 S_c *= units.kB

358 write(fmt % ('S_con', S_c, S_c * temperature))

359 S += S_c

361 write('-' * 49)

362 write(fmt % ('S', S, S * temperature))

363 write('=' * 49)

364 return S

366 def get_helmholtz_energy(self, temperature, verbose=True):

367 """Returns the Helmholtz free energy, in eV, in the hindered

368 translator and hindered rotor model at a specified temperature

369 (K)."""

371 self.verbose = True

372 write = self._vprint

374 U = self.get_internal_energy(temperature, verbose=verbose)

375 write('')

376 S = self.get_entropy(temperature, verbose=verbose)

377 F = U - temperature * S

379 write('')

380 write('Free energy components at T = %.2f K:' % temperature)

381 write('=' * 23)

382 fmt = '%5s%15.3f eV'

383 write(fmt % ('U', U))

384 write(fmt % ('-T*S', -temperature * S))

385 write('-' * 23)

386 write(fmt % ('F', F))

387 write('=' * 23)

388 return F

391class IdealGasThermo(ThermoChem):

392 """Class for calculating thermodynamic properties of a molecule

393 based on statistical mechanical treatments in the ideal gas

394 approximation.

396 Inputs for enthalpy calculations:

398 vib_energies : list

399 a list of the vibrational energies of the molecule (e.g., from

400 ase.vibrations.Vibrations.get_energies). The number of vibrations

401 used is automatically calculated by the geometry and the number of

402 atoms. If more are specified than are needed, then the lowest

403 numbered vibrations are neglected. If either atoms or natoms is

404 unspecified, then uses the entire list. Units are eV.

405 geometry : 'monatomic', 'linear', or 'nonlinear'

406 geometry of the molecule

407 potentialenergy : float

408 the potential energy in eV (e.g., from atoms.get_potential_energy)

409 (if potentialenergy is unspecified, then the methods of this

410 class can be interpreted as the energy corrections)

411 natoms : integer

412 the number of atoms, used along with 'geometry' to determine how

413 many vibrations to use. (Not needed if an atoms object is supplied

414 in 'atoms' or if the user desires the entire list of vibrations

415 to be used.)

417 Extra inputs needed for entropy / free energy calculations:

419 atoms : an ASE atoms object

420 used to calculate rotational moments of inertia and molecular mass

421 symmetrynumber : integer

422 symmetry number of the molecule. See, for example, Table 10.1 and

423 Appendix B of C. Cramer "Essentials of Computational Chemistry",

424 2nd Ed.

425 spin : float

426 the total electronic spin. (0 for molecules in which all electrons

427 are paired, 0.5 for a free radical with a single unpaired electron,

428 1.0 for a triplet with two unpaired electrons, such as O_2.)

429 """

431 def __init__(self, vib_energies, geometry, potentialenergy=0.,

432 atoms=None, symmetrynumber=None, spin=None, natoms=None):

433 self.potentialenergy = potentialenergy

434 self.geometry = geometry

435 self.atoms = atoms

436 self.sigma = symmetrynumber

437 self.spin = spin

438 if natoms is None:

439 if atoms:

440 natoms = len(atoms)

441 # Cut the vibrations to those needed from the geometry.

442 if natoms:

443 if geometry == 'nonlinear':

444 self.vib_energies = vib_energies[-(3 * natoms - 6):]

445 elif geometry == 'linear':

446 self.vib_energies = vib_energies[-(3 * natoms - 5):]

447 elif geometry == 'monatomic':

448 self.vib_energies = []

449 else:

450 self.vib_energies = vib_energies

451 # Make sure no imaginary frequencies remain.

452 if sum(np.iscomplex(self.vib_energies)):

453 raise ValueError('Imaginary frequencies are present.')

454 else:

455 self.vib_energies = np.real(self.vib_energies) # clear +0.j

456 self.referencepressure = 1.0e5 # Pa

458 def get_enthalpy(self, temperature, verbose=True):

459 """Returns the enthalpy, in eV, in the ideal gas approximation

460 at a specified temperature (K)."""

462 self.verbose = verbose

463 write = self._vprint

464 fmt = '%-15s%13.3f eV'

465 write('Enthalpy components at T = %.2f K:' % temperature)

466 write('=' * 31)

468 H = 0.

470 write(fmt % ('E_pot', self.potentialenergy))

471 H += self.potentialenergy

473 zpe = self.get_ZPE_correction()

474 write(fmt % ('E_ZPE', zpe))

475 H += zpe

477 Cv_t = 3. / 2. * units.kB # translational heat capacity (3-d gas)

478 write(fmt % ('Cv_trans (0->T)', Cv_t * temperature))

479 H += Cv_t * temperature

481 if self.geometry == 'nonlinear': # rotational heat capacity

482 Cv_r = 3. / 2. * units.kB

483 elif self.geometry == 'linear':

484 Cv_r = units.kB

485 elif self.geometry == 'monatomic':

486 Cv_r = 0.

487 write(fmt % ('Cv_rot (0->T)', Cv_r * temperature))

488 H += Cv_r * temperature

490 dH_v = self._vibrational_energy_contribution(temperature)

491 write(fmt % ('Cv_vib (0->T)', dH_v))

492 H += dH_v

494 Cp_corr = units.kB * temperature

495 write(fmt % ('(C_v -> C_p)', Cp_corr))

496 H += Cp_corr

498 write('-' * 31)

499 write(fmt % ('H', H))

500 write('=' * 31)

501 return H

503 def get_entropy(self, temperature, pressure, verbose=True):

504 """Returns the entropy, in eV/K, in the ideal gas approximation

505 at a specified temperature (K) and pressure (Pa)."""

507 if self.atoms is None or self.sigma is None or self.spin is None:

508 raise RuntimeError('atoms, symmetrynumber, and spin must be '

509 'specified for entropy and free energy '

510 'calculations.')

511 self.verbose = verbose

512 write = self._vprint

513 fmt = '%-15s%13.7f eV/K%13.3f eV'

514 write('Entropy components at T = %.2f K and P = %.1f Pa:' %

515 (temperature, pressure))

516 write('=' * 49)

517 write('%15s%13s %13s' % ('', 'S', 'T*S'))

519 S = 0.0

521 # Translational entropy (term inside the log is in SI units).

522 mass = sum(self.atoms.get_masses()) * units._amu # kg/molecule

523 S_t = (2 * np.pi * mass * units._k *

524 temperature / units._hplanck**2)**(3.0 / 2)

525 S_t *= units._k * temperature / self.referencepressure

526 S_t = units.kB * (np.log(S_t) + 5.0 / 2.0)

527 write(fmt % ('S_trans (1 bar)', S_t, S_t * temperature))

528 S += S_t

530 # Rotational entropy (term inside the log is in SI units).

531 if self.geometry == 'monatomic':

532 S_r = 0.0

533 elif self.geometry == 'nonlinear':

534 inertias = (self.atoms.get_moments_of_inertia() * units._amu /

535 (10.0**10)**2) # kg m^2

536 S_r = np.sqrt(np.pi * np.product(inertias)) / self.sigma

537 S_r *= (8.0 * np.pi**2 * units._k * temperature /

538 units._hplanck**2)**(3.0 / 2.0)

539 S_r = units.kB * (np.log(S_r) + 3.0 / 2.0)

540 elif self.geometry == 'linear':

541 inertias = (self.atoms.get_moments_of_inertia() * units._amu /

542 (10.0**10)**2) # kg m^2

543 inertia = max(inertias) # should be two identical and one zero

544 S_r = (8 * np.pi**2 * inertia * units._k * temperature /

545 self.sigma / units._hplanck**2)

546 S_r = units.kB * (np.log(S_r) + 1.)

547 write(fmt % ('S_rot', S_r, S_r * temperature))

548 S += S_r

550 # Electronic entropy.

551 S_e = units.kB * np.log(2 * self.spin + 1)

552 write(fmt % ('S_elec', S_e, S_e * temperature))

553 S += S_e

555 # Vibrational entropy.

556 S_v = self._vibrational_entropy_contribution(temperature)

557 write(fmt % ('S_vib', S_v, S_v * temperature))

558 S += S_v

560 # Pressure correction to translational entropy.

561 S_p = - units.kB * np.log(pressure / self.referencepressure)

562 write(fmt % ('S (1 bar -> P)', S_p, S_p * temperature))

563 S += S_p

565 write('-' * 49)

566 write(fmt % ('S', S, S * temperature))

567 write('=' * 49)

568 return S

570 def get_gibbs_energy(self, temperature, pressure, verbose=True):

571 """Returns the Gibbs free energy, in eV, in the ideal gas

572 approximation at a specified temperature (K) and pressure (Pa)."""

574 self.verbose = verbose

575 write = self._vprint

577 H = self.get_enthalpy(temperature, verbose=verbose)

578 write('')

579 S = self.get_entropy(temperature, pressure, verbose=verbose)

580 G = H - temperature * S

582 write('')

583 write('Free energy components at T = %.2f K and P = %.1f Pa:' %

584 (temperature, pressure))

585 write('=' * 23)

586 fmt = '%5s%15.3f eV'

587 write(fmt % ('H', H))

588 write(fmt % ('-T*S', -temperature * S))

589 write('-' * 23)

590 write(fmt % ('G', G))

591 write('=' * 23)

592 return G

595class CrystalThermo(ThermoChem):

596 """Class for calculating thermodynamic properties of a crystalline

597 solid in the approximation that a lattice of N atoms behaves as a

598 system of 3N independent harmonic oscillators.

600 Inputs:

602 phonon_DOS : list

603 a list of the phonon density of states,

604 where each value represents the phonon DOS at the vibrational energy

605 value of the corresponding index in phonon_energies.

607 phonon_energies : list

608 a list of the range of vibrational energies (hbar*omega) over which

609 the phonon density of states has been evaluated. This list should be

610 the same length as phonon_DOS and integrating phonon_DOS over

611 phonon_energies should yield approximately 3N, where N is the number

612 of atoms per unit cell. If the first element of this list is

613 zero-valued it will be deleted along with the first element of

614 phonon_DOS. Units of vibrational energies are eV.

616 potentialenergy : float

617 the potential energy in eV (e.g., from atoms.get_potential_energy)

618 (if potentialenergy is unspecified, then the methods of this

619 class can be interpreted as the energy corrections)

621 formula_units : int

622 the number of formula units per unit cell. If unspecified, the

623 thermodynamic quantities calculated will be listed on a

624 per-unit-cell basis.

625 """

627 def __init__(self, phonon_DOS, phonon_energies,

628 formula_units=None, potentialenergy=0.):

629 self.phonon_energies = phonon_energies

630 self.phonon_DOS = phonon_DOS

632 if formula_units:

633 self.formula_units = formula_units

634 self.potentialenergy = potentialenergy / formula_units

635 else:

636 self.formula_units = 0

637 self.potentialenergy = potentialenergy

639 def get_internal_energy(self, temperature, verbose=True):

640 """Returns the internal energy, in eV, of crystalline solid

641 at a specified temperature (K)."""

643 self.verbose = verbose

644 write = self._vprint

645 fmt = '%-15s%13.4f eV'

646 if self.formula_units == 0:

647 write('Internal energy components at '

648 'T = %.2f K,\non a per-unit-cell basis:' % temperature)

649 else:

650 write('Internal energy components at '

651 'T = %.2f K,\non a per-formula-unit basis:' % temperature)

652 write('=' * 31)

654 U = 0.

656 omega_e = self.phonon_energies

657 dos_e = self.phonon_DOS

658 if omega_e[0] == 0.:

659 omega_e = np.delete(omega_e, 0)

660 dos_e = np.delete(dos_e, 0)

662 write(fmt % ('E_pot', self.potentialenergy))

663 U += self.potentialenergy

665 zpe_list = omega_e / 2.

666 if self.formula_units == 0:

667 zpe = np.trapz(zpe_list * dos_e, omega_e)

668 else:

669 zpe = np.trapz(zpe_list * dos_e, omega_e) / self.formula_units

670 write(fmt % ('E_ZPE', zpe))

671 U += zpe

673 B = 1. / (units.kB * temperature)

674 E_vib = omega_e / (np.exp(omega_e * B) - 1.)

675 if self.formula_units == 0:

676 E_phonon = np.trapz(E_vib * dos_e, omega_e)

677 else:

678 E_phonon = np.trapz(E_vib * dos_e, omega_e) / self.formula_units

679 write(fmt % ('E_phonon', E_phonon))

680 U += E_phonon

682 write('-' * 31)

683 write(fmt % ('U', U))

684 write('=' * 31)

685 return U

687 def get_entropy(self, temperature, verbose=True):

688 """Returns the entropy, in eV/K, of crystalline solid

689 at a specified temperature (K)."""

691 self.verbose = verbose

692 write = self._vprint

693 fmt = '%-15s%13.7f eV/K%13.4f eV'

694 if self.formula_units == 0:

695 write('Entropy components at '

696 'T = %.2f K,\non a per-unit-cell basis:' % temperature)

697 else:

698 write('Entropy components at '

699 'T = %.2f K,\non a per-formula-unit basis:' % temperature)

700 write('=' * 49)

701 write('%15s%13s %13s' % ('', 'S', 'T*S'))

703 omega_e = self.phonon_energies

704 dos_e = self.phonon_DOS

705 if omega_e[0] == 0.:

706 omega_e = np.delete(omega_e, 0)

707 dos_e = np.delete(dos_e, 0)

709 B = 1. / (units.kB * temperature)

710 S_vib = (omega_e / (temperature * (np.exp(omega_e * B) - 1.)) -

711 units.kB * np.log(1. - np.exp(-omega_e * B)))

712 if self.formula_units == 0:

713 S = np.trapz(S_vib * dos_e, omega_e)

714 else:

715 S = np.trapz(S_vib * dos_e, omega_e) / self.formula_units

717 write('-' * 49)

718 write(fmt % ('S', S, S * temperature))

719 write('=' * 49)

720 return S

722 def get_helmholtz_energy(self, temperature, verbose=True):

723 """Returns the Helmholtz free energy, in eV, of crystalline solid

724 at a specified temperature (K)."""

726 self.verbose = True

727 write = self._vprint

729 U = self.get_internal_energy(temperature, verbose=verbose)

730 write('')

731 S = self.get_entropy(temperature, verbose=verbose)

732 F = U - temperature * S

734 write('')

735 if self.formula_units == 0:

736 write('Helmholtz free energy components at '

737 'T = %.2f K,\non a per-unit-cell basis:' % temperature)

738 else:

739 write('Helmholtz free energy components at '

740 'T = %.2f K,\non a per-formula-unit basis:' % temperature)

741 write('=' * 23)

742 fmt = '%5s%15.4f eV'

743 write(fmt % ('U', U))

744 write(fmt % ('-T*S', -temperature * S))

745 write('-' * 23)

746 write(fmt % ('F', F))

747 write('=' * 23)

748 return F