Coverage for /builds/ase/ase/ase/thermochemistry.py: 94.23%

1# fmt: off 

2 

3"""Modules for calculating thermochemical information from computational 

4outputs.""" 

5 

6import os 

7import sys 

8from warnings import warn 

9 

10import numpy as np 

11from scipy.integrate import trapezoid 

12 

13from ase import units 

14 

15 

16class ThermoChem: 

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

18 calculations.""" 

19 

20 def get_ZPE_correction(self): 

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

22 return 0.5 * np.sum(self.vib_energies) 

23 

24 def _vibrational_energy_contribution(self, temperature): 

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

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

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

28 in eV.""" 

29 kT = units.kB * temperature 

30 dU = 0. 

31 for energy in self.vib_energies: 

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

33 return dU 

34 

35 def _vibrational_entropy_contribution(self, temperature): 

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

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

38 in eV/K.""" 

39 kT = units.kB * temperature 

40 S_v = 0. 

41 for energy in self.vib_energies: 

42 x = energy / kT 

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

44 S_v *= units.kB 

45 return S_v 

46 

47 def _vprint(self, text): 

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

49 if self.verbose: 

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

51 

52 

53class HarmonicThermo(ThermoChem): 

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

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

56 adsorbates. 

57 

58 Inputs: 

59 

60 vib_energies : list 

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

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

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

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

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

66 number of energies. Units of energies are eV. 

67 potentialenergy : float 

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

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

70 class can be interpreted as the energy corrections) 

71 ignore_imag_modes : bool 

72 If True, any imaginary frequencies will be ignored in the calculation 

73 of the thermochemical properties. If False (default), an error will 

74 be raised if any imaginary frequencies are present. 

75 """ 

76 

77 def __init__(self, vib_energies, potentialenergy=0., 

78 ignore_imag_modes=False): 

79 self.ignore_imag_modes = ignore_imag_modes 

80 

81 # Check for imaginary frequencies. 

82 vib_energies, n_imag = _clean_vib_energies( 

83 vib_energies, ignore_imag_modes=ignore_imag_modes 

84 ) 

85 self.vib_energies = vib_energies 

86 self.n_imag = n_imag 

87 

88 self.potentialenergy = potentialenergy 

89 

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

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

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

93 

94 self.verbose = verbose 

95 write = self._vprint 

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

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

98 write('=' * 31) 

99 

100 U = 0. 

101 

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

103 U += self.potentialenergy 

104 

105 zpe = self.get_ZPE_correction() 

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

107 U += zpe 

108 

109 dU_v = self._vibrational_energy_contribution(temperature) 

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

111 U += dU_v 

112 

113 write('-' * 31) 

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

115 write('=' * 31) 

116 return U 

117 

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

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

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

121 

122 self.verbose = verbose 

123 write = self._vprint 

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

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

126 write('=' * 49) 

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

128 

129 S = 0. 

130 

131 S_v = self._vibrational_entropy_contribution(temperature) 

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

133 S += S_v 

134 

135 write('-' * 49) 

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

137 write('=' * 49) 

138 return S 

139 

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

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

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

143 

144 self.verbose = True 

145 write = self._vprint 

146 

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

148 write('') 

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

150 F = U - temperature * S 

151 

152 write('') 

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

154 write('=' * 23) 

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

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

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

158 write('-' * 23) 

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

160 write('=' * 23) 

161 return F 

162 

163 

164class HinderedThermo(ThermoChem): 

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

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

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

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

169 

170 Inputs: 

171 

172 vib_energies : list 

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

174 ase.vibrations.Vibrations.get_energies). If atoms is not provided, 

175 then the number of energies must match the number of degrees of freedom 

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

177 that this class does not check that the user has supplied the 

178 correct number of energies. 

179 Units of energies are eV. 

180 trans_barrier_energy : float 

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

182 adsorbate to diffuse on the surface. 

183 rot_barrier_energy : float 

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

185 adsorbate to rotate about an axis perpendicular to the surface. 

186 sitedensity : float 

187 density of surface sites in cm^-2 

188 rotationalminima : integer 

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

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

191 potentialenergy : float 

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

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

194 can be interpreted as the energy corrections) 

195 mass : float 

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

197 be calculated from the atoms class) 

198 inertia : float 

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

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

201 atoms class) 

202 atoms : an ASE atoms object 

203 used to calculate rotational moments of inertia and molecular mass 

204 symmetrynumber : integer 

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

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

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

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

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

210 ignore_imag_modes : bool 

211 If True, any imaginary frequencies present after the 3N-3 cut will not 

212 be included in the calculation of the thermochemical properties. If 

213 False (default), an error will be raised if imaginary frequencies are 

214 present after the 3N-3 cut. 

215 """ 

216 

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

218 sitedensity, rotationalminima, potentialenergy=0., 

219 mass=None, inertia=None, atoms=None, symmetrynumber=1, 

220 ignore_imag_modes=False): 

221 

222 self.trans_barrier_energy = trans_barrier_energy * units._e 

223 self.rot_barrier_energy = rot_barrier_energy * units._e 

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

225 self.rotationalminima = rotationalminima 

226 self.potentialenergy = potentialenergy 

227 self.atoms = atoms 

228 self.symmetry = symmetrynumber 

229 self.ignore_imag_modes = ignore_imag_modes 

230 

231 # Sort the vibrations 

232 vib_energies = list(vib_energies) 

233 vib_energies.sort(key=np.abs) 

234 

235 # Keep only the relevant vibrational energies (3N-3) 

236 if atoms: 

237 vib_energies = vib_energies[-(3 * len(atoms) - 3):] 

238 else: 

239 vib_energies = vib_energies[-(len(vib_energies) - 3):] 

240 

241 # Check for imaginary frequencies. 

242 vib_energies, n_imag = _clean_vib_energies( 

243 vib_energies, ignore_imag_modes=ignore_imag_modes 

244 ) 

245 self.vib_energies = vib_energies 

246 self.n_imag = n_imag 

247 

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

249 if mass: 

250 self.mass = mass * units._amu 

251 elif atoms: 

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

253 if inertia: 

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

255 elif atoms: 

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

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

258 else: 

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

260 'adsorbate must be specified or ' 

261 'atoms must be specified.') 

262 

263 # Calculate hindered translational and rotational frequencies 

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

265 self.area)) 

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

267 self.rot_barrier_energy / 

268 (2 * self.inertia)) 

269 

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

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

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

273 specified temperature (K).""" 

274 

275 from scipy.special import iv 

276 

277 self.verbose = verbose 

278 write = self._vprint 

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

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

281 write('=' * 31) 

282 

283 U = 0. 

284 

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

286 U += self.potentialenergy 

287 

288 # Translational Energy 

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

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

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

292 R_t / 2 / T_t * 

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

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

295 dU_t *= units.kB * temperature 

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

297 U += dU_t 

298 

299 # Rotational Energy 

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

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

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

303 R_r / 2 / T_r * 

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

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

306 dU_r *= units.kB * temperature 

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

308 U += dU_r 

309 

310 # Vibrational Energy 

311 dU_v = self._vibrational_energy_contribution(temperature) 

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

313 U += dU_v 

314 

315 # Zero Point Energy 

316 dU_zpe = self.get_zero_point_energy() 

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

318 U += dU_zpe 

319 

320 write('-' * 31) 

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

322 write('=' * 31) 

323 return U 

324 

325 def get_zero_point_energy(self, verbose=True): 

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

327 translator and hindered rotor model""" 

328 

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

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

331 zpe_v = self.get_ZPE_correction() 

332 zpe = zpe_t + zpe_r + zpe_v 

333 return zpe 

334 

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

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

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

338 

339 from scipy.special import iv 

340 

341 self.verbose = verbose 

342 write = self._vprint 

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

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

345 write('=' * 49) 

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

347 

348 S = 0. 

349 

350 # Translational Entropy 

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

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

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

354 R_t / 2 / T_t * 

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

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

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

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

359 S_t *= units.kB 

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

361 S += S_t 

362 

363 # Rotational Entropy 

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

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

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

367 np.log(self.symmetry) - 

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

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

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

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

372 S_r *= units.kB 

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

374 S += S_r 

375 

376 # Vibrational Entropy 

377 S_v = self._vibrational_entropy_contribution(temperature) 

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

379 S += S_v 

380 

381 # Concentration Related Entropy 

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

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

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

385 S_c *= units.kB 

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

387 S += S_c 

388 

389 write('-' * 49) 

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

391 write('=' * 49) 

392 return S 

393 

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

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

396 translator and hindered rotor model at a specified temperature 

397 (K).""" 

398 

399 self.verbose = True 

400 write = self._vprint 

401 

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

403 write('') 

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

405 F = U - temperature * S 

406 

407 write('') 

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

409 write('=' * 23) 

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

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

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

413 write('-' * 23) 

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

415 write('=' * 23) 

416 return F 

417 

418 

419class IdealGasThermo(ThermoChem): 

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

421 based on statistical mechanical treatments in the ideal gas 

422 approximation. 

423 

424 Inputs for enthalpy calculations: 

425 

426 vib_energies : list 

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

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

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

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

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

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

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

434 geometry of the molecule 

435 potentialenergy : float 

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

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

438 class can be interpreted as the energy corrections) 

439 natoms : integer 

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

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

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

443 to be used.) 

444 

445 Extra inputs needed for entropy / free energy calculations: 

446 

447 atoms : an ASE atoms object 

448 used to calculate rotational moments of inertia and molecular mass 

449 symmetrynumber : integer 

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

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

452 2nd Ed. 

453 spin : float 

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

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

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

457 ignore_imag_modes : bool 

458 If True, any imaginary frequencies present after the 3N-5/3N-6 cut 

459 will not be included in the calculation of the thermochemical 

460 properties. If False (default), a ValueError will be raised if 

461 any imaginary frequencies remain after the 3N-5/3N-6 cut. 

462 """ 

463 

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

465 atoms=None, symmetrynumber=None, spin=None, natoms=None, 

466 ignore_imag_modes=False): 

467 self.potentialenergy = potentialenergy 

468 self.geometry = geometry 

469 self.atoms = atoms 

470 self.sigma = symmetrynumber 

471 self.spin = spin 

472 self.ignore_imag_modes = ignore_imag_modes 

473 if natoms is None and atoms: 

474 natoms = len(atoms) 

475 self.natoms = natoms 

476 

477 # Sort the vibrations 

478 vib_energies = list(vib_energies) 

479 vib_energies.sort(key=np.abs) 

480 

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

482 if natoms: 

483 if geometry == 'nonlinear': 

484 vib_energies = vib_energies[-(3 * natoms - 6):] 

485 elif geometry == 'linear': 

486 vib_energies = vib_energies[-(3 * natoms - 5):] 

487 elif geometry == 'monatomic': 

488 vib_energies = [] 

489 else: 

490 raise ValueError(f"Unsupported geometry: {geometry}") 

491 

492 # Check for imaginary frequencies. 

493 vib_energies, n_imag = _clean_vib_energies( 

494 vib_energies, ignore_imag_modes=ignore_imag_modes 

495 ) 

496 self.vib_energies = vib_energies 

497 self.n_imag = n_imag 

498 

499 self.referencepressure = 1.0e5 # Pa 

500 

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

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

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

504 

505 self.verbose = verbose 

506 write = self._vprint 

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

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

509 write('=' * 31) 

510 

511 H = 0. 

512 

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

514 H += self.potentialenergy 

515 

516 zpe = self.get_ZPE_correction() 

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

518 H += zpe 

519 

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

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

522 H += Cv_t * temperature 

523 

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

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

526 elif self.geometry == 'linear': 

527 Cv_r = units.kB 

528 elif self.geometry == 'monatomic': 

529 Cv_r = 0. 

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

531 H += Cv_r * temperature 

532 

533 dH_v = self._vibrational_energy_contribution(temperature) 

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

535 H += dH_v 

536 

537 Cp_corr = units.kB * temperature 

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

539 H += Cp_corr 

540 

541 write('-' * 31) 

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

543 write('=' * 31) 

544 return H 

545 

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

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

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

549 

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

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

552 'specified for entropy and free energy ' 

553 'calculations.') 

554 self.verbose = verbose 

555 write = self._vprint 

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

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

558 (temperature, pressure)) 

559 write('=' * 49) 

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

561 

562 S = 0.0 

563 

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

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

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

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

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

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

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

571 S += S_t 

572 

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

574 if self.geometry == 'monatomic': 

575 S_r = 0.0 

576 elif self.geometry == 'nonlinear': 

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

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

579 S_r = np.sqrt(np.pi * np.prod(inertias)) / self.sigma 

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

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

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

583 elif self.geometry == 'linear': 

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

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

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

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

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

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

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

591 S += S_r 

592 

593 # Electronic entropy. 

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

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

596 S += S_e 

597 

598 # Vibrational entropy. 

599 S_v = self._vibrational_entropy_contribution(temperature) 

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

601 S += S_v 

602 

603 # Pressure correction to translational entropy. 

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

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

606 S += S_p 

607 

608 write('-' * 49) 

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

610 write('=' * 49) 

611 return S 

612 

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

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

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

616 

617 self.verbose = verbose 

618 write = self._vprint 

619 

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

621 write('') 

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

623 G = H - temperature * S 

624 

625 write('') 

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

627 (temperature, pressure)) 

628 write('=' * 23) 

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

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

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

632 write('-' * 23) 

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

634 write('=' * 23) 

635 return G 

636 

637 

638class CrystalThermo(ThermoChem): 

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

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

641 system of 3N independent harmonic oscillators. 

642 

643 Inputs: 

644 

645 phonon_DOS : list 

646 a list of the phonon density of states, 

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

648 value of the corresponding index in phonon_energies. 

649 

650 phonon_energies : list 

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

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

653 the same length as phonon_DOS and integrating phonon_DOS over 

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

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

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

657 phonon_DOS. Units of vibrational energies are eV. 

658 

659 potentialenergy : float 

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

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

662 class can be interpreted as the energy corrections) 

663 

664 formula_units : int 

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

666 thermodynamic quantities calculated will be listed on a 

667 per-unit-cell basis. 

668 """ 

669 

670 def __init__(self, phonon_DOS, phonon_energies, 

671 formula_units=None, potentialenergy=0.): 

672 self.phonon_energies = phonon_energies 

673 self.phonon_DOS = phonon_DOS 

674 

675 if formula_units: 

676 self.formula_units = formula_units 

677 self.potentialenergy = potentialenergy / formula_units 

678 else: 

679 self.formula_units = 0 

680 self.potentialenergy = potentialenergy 

681 

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

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

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

685 

686 self.verbose = verbose 

687 write = self._vprint 

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

689 if self.formula_units == 0: 

690 write('Internal energy components at ' 

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

692 else: 

693 write('Internal energy components at ' 

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

695 write('=' * 31) 

696 

697 U = 0. 

698 

699 omega_e = self.phonon_energies 

700 dos_e = self.phonon_DOS 

701 if omega_e[0] == 0.: 

702 omega_e = np.delete(omega_e, 0) 

703 dos_e = np.delete(dos_e, 0) 

704 

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

706 U += self.potentialenergy 

707 

708 zpe_list = omega_e / 2. 

709 if self.formula_units == 0: 

710 zpe = trapezoid(zpe_list * dos_e, omega_e) 

711 else: 

712 zpe = trapezoid(zpe_list * dos_e, omega_e) / self.formula_units 

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

714 U += zpe 

715 

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

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

718 if self.formula_units == 0: 

719 E_phonon = trapezoid(E_vib * dos_e, omega_e) 

720 else: 

721 E_phonon = trapezoid(E_vib * dos_e, omega_e) / self.formula_units 

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

723 U += E_phonon 

724 

725 write('-' * 31) 

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

727 write('=' * 31) 

728 return U 

729 

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

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

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

733 

734 self.verbose = verbose 

735 write = self._vprint 

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

737 if self.formula_units == 0: 

738 write('Entropy components at ' 

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

740 else: 

741 write('Entropy components at ' 

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

743 write('=' * 49) 

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

745 

746 omega_e = self.phonon_energies 

747 dos_e = self.phonon_DOS 

748 if omega_e[0] == 0.: 

749 omega_e = np.delete(omega_e, 0) 

750 dos_e = np.delete(dos_e, 0) 

751 

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

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

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

755 if self.formula_units == 0: 

756 S = trapezoid(S_vib * dos_e, omega_e) 

757 else: 

758 S = trapezoid(S_vib * dos_e, omega_e) / self.formula_units 

759 

760 write('-' * 49) 

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

762 write('=' * 49) 

763 return S 

764 

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

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

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

768 

769 self.verbose = True 

770 write = self._vprint 

771 

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

773 write('') 

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

775 F = U - temperature * S 

776 

777 write('') 

778 if self.formula_units == 0: 

779 write('Helmholtz free energy components at ' 

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

781 else: 

782 write('Helmholtz free energy components at ' 

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

784 write('=' * 23) 

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

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

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

788 write('-' * 23) 

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

790 write('=' * 23) 

791 return F 

792 

793 

794def _clean_vib_energies(vib_energies, ignore_imag_modes=False): 

795 """Checks for presence of imaginary vibrational modes and removes 

796 them if ignore_imag_modes is True. Also removes +0.j from real 

797 vibrational energies. 

798 

799 Inputs: 

800 

801 vib_energies : list 

802 a list of the vibrational energies 

803 

804 ignore_imag_modes : bool 

805 If True, any imaginary frequencies will be removed. If False, 

806 (default), an error will be raised if imaginary frequencies are 

807 present. 

808 

809 Outputs: 

810 

811 vib_energies : list 

812 a list of the cleaned vibrational energies with imaginary frequencies 

813 removed if ignore_imag_modes is True.