1"""Modules for calculating thermochemical information from computational 

2outputs.""" 

3 

4import os 

5import sys 

6import numpy as np 

7 

8from ase import units 

9 

10 

11class ThermoChem: 

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

13 calculations.""" 

14 

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 

21 

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 

32 

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 

44 

45 def _vprint(self, text): 

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

47 if self.verbose: 

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

49 

50 

51class HarmonicThermo(ThermoChem): 

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

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

54 adsorbates. 

55 

56 Inputs: 

57 

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

70 

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 

78 

79 self.potentialenergy = potentialenergy 

80 

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

84 

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) 

90 

91 U = 0. 

92 

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

94 U += self.potentialenergy 

95 

96 zpe = self.get_ZPE_correction() 

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

98 U += zpe 

99 

100 dU_v = self._vibrational_energy_contribution(temperature) 

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

102 U += dU_v 

103 

104 write('-' * 31) 

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

106 write('=' * 31) 

107 return U 

108 

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

112 

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

119 

120 S = 0. 

121 

122 S_v = self._vibrational_entropy_contribution(temperature) 

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

124 S += S_v 

125 

126 write('-' * 49) 

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

128 write('=' * 49) 

129 return S 

130 

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

134 

135 self.verbose = True 

136 write = self._vprint 

137 

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

139 write('') 

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

141 F = U - temperature * S 

142 

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 

153 

154 

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. 

160 

161 Inputs: 

162 

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 

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

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

201 

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 

213 

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

228 

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 

234 

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

241 

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

246 

247 from scipy.special import iv 

248 

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) 

254 

255 U = 0. 

256 

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

258 U += self.potentialenergy 

259 

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 

270 

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 

281 

282 # Vibrational Energy 

283 dU_v = self._vibrational_energy_contribution(temperature) 

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

285 U += dU_v 

286 

287 # Zero Point Energy 

288 dU_zpe = self.get_zero_point_energy() 

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

290 U += dU_zpe 

291 

292 write('-' * 31) 

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

294 write('=' * 31) 

295 return U 

296 

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

300 

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 

306 

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

310 

311 from scipy.special import iv 

312 

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

319 

320 S = 0. 

321 

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 

334 

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 

347 

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 

352 

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 

360 

361 write('-' * 49) 

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

363 write('=' * 49) 

364 return S 

365 

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

370 

371 self.verbose = True 

372 write = self._vprint 

373 

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

375 write('') 

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

377 F = U - temperature * S 

378 

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 

389 

390 

391class IdealGasThermo(ThermoChem): 

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

393 based on statistical mechanical treatments in the ideal gas 

394 approximation. 

395 

396 Inputs for enthalpy calculations: 

397 

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

416 

417 Extra inputs needed for entropy / free energy calculations: 

418 

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

430 

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 

457 

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

461 

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) 

467 

468 H = 0. 

469 

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

471 H += self.potentialenergy 

472 

473 zpe = self.get_ZPE_correction() 

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

475 H += zpe 

476 

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 

480 

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 

489 

490 dH_v = self._vibrational_energy_contribution(temperature) 

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

492 H += dH_v 

493 

494 Cp_corr = units.kB * temperature 

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

496 H += Cp_corr 

497 

498 write('-' * 31) 

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

500 write('=' * 31) 

501 return H 

502 

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

506 

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

518 

519 S = 0.0 

520 

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 

529 

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 

549 

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 

554 

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 

559 

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 

564 

565 write('-' * 49) 

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

567 write('=' * 49) 

568 return S 

569 

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

573 

574 self.verbose = verbose 

575 write = self._vprint 

576 

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

578 write('') 

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

580 G = H - temperature * S 

581 

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 

593 

594 

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. 

599 

600 Inputs: 

601 

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. 

606 

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. 

615 

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) 

620 

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

626 

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 

631 

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 

638 

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

642 

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) 

653 

654 U = 0. 

655 

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) 

661 

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

663 U += self.potentialenergy 

664 

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 

672 

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 

681 

682 write('-' * 31) 

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

684 write('=' * 31) 

685 return U 

686 

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

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

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

690 

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

702 

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) 

708 

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 

716 

717 write('-' * 49) 

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

719 write('=' * 49) 

720 return S 

721 

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

725 

726 self.verbose = True 

727 write = self._vprint 

728 

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

730 write('') 

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

732 F = U - temperature * S 

733 

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