Coverage for /builds/kinetik161/ase/ase/thermochemistry.py: 95.18%
415 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Modules for calculating thermochemical information from computational
2outputs."""
4import os
5import sys
6from warnings import warn
8import numpy as np
10from ase import units
13class ThermoChem:
14 """Base class containing common methods used in thermochemistry
15 calculations."""
17 def get_ZPE_correction(self):
18 """Returns the zero-point vibrational energy correction in eV."""
19 return 0.5 * np.sum(self.vib_energies)
21 def _vibrational_energy_contribution(self, temperature):
22 """Calculates the change in internal energy due to vibrations from
23 0K to the specified temperature for a set of vibrations given in
24 eV and a temperature given in Kelvin. Returns the energy change
25 in eV."""
26 kT = units.kB * temperature
27 dU = 0.
28 for energy in self.vib_energies:
29 dU += energy / (np.exp(energy / kT) - 1.)
30 return dU
32 def _vibrational_entropy_contribution(self, temperature):
33 """Calculates the entropy due to vibrations for a set of vibrations
34 given in eV and a temperature given in Kelvin. Returns the entropy
35 in eV/K."""
36 kT = units.kB * temperature
37 S_v = 0.
38 for energy in self.vib_energies:
39 x = energy / kT
40 S_v += x / (np.exp(x) - 1.) - np.log(1. - np.exp(-x))
41 S_v *= units.kB
42 return S_v
44 def _vprint(self, text):
45 """Print output if verbose flag True."""
46 if self.verbose:
47 sys.stdout.write(text + os.linesep)
50class HarmonicThermo(ThermoChem):
51 """Class for calculating thermodynamic properties in the approximation
52 that all degrees of freedom are treated harmonically. Often used for
53 adsorbates.
55 Inputs:
57 vib_energies : list
58 a list of the harmonic energies of the adsorbate (e.g., from
59 ase.vibrations.Vibrations.get_energies). The number of
60 energies should match the number of degrees of freedom of the
61 adsorbate; i.e., 3*n, where n is the number of atoms. Note that
62 this class does not check that the user has supplied the correct
63 number of energies. Units of energies are eV.
64 potentialenergy : float
65 the potential energy in eV (e.g., from atoms.get_potential_energy)
66 (if potentialenergy is unspecified, then the methods of this
67 class can be interpreted as the energy corrections)
68 ignore_imag_modes : bool
69 If True, any imaginary frequencies will be ignored in the calculation
70 of the thermochemical properties. If False (default), an error will
71 be raised if any imaginary frequencies are present.
72 """
74 def __init__(self, vib_energies, potentialenergy=0.,
75 ignore_imag_modes=False):
76 self.ignore_imag_modes = ignore_imag_modes
78 # Check for imaginary frequencies.
79 vib_energies, n_imag = _clean_vib_energies(
80 vib_energies, ignore_imag_modes=ignore_imag_modes
81 )
82 self.vib_energies = vib_energies
83 self.n_imag = n_imag
85 self.potentialenergy = potentialenergy
87 def get_internal_energy(self, temperature, verbose=True):
88 """Returns the internal energy, in eV, in the harmonic approximation
89 at a specified temperature (K)."""
91 self.verbose = verbose
92 write = self._vprint
93 fmt = '%-15s%13.3f eV'
94 write('Internal energy components at T = %.2f K:' % temperature)
95 write('=' * 31)
97 U = 0.
99 write(fmt % ('E_pot', self.potentialenergy))
100 U += self.potentialenergy
102 zpe = self.get_ZPE_correction()
103 write(fmt % ('E_ZPE', zpe))
104 U += zpe
106 dU_v = self._vibrational_energy_contribution(temperature)
107 write(fmt % ('Cv_harm (0->T)', dU_v))
108 U += dU_v
110 write('-' * 31)
111 write(fmt % ('U', U))
112 write('=' * 31)
113 return U
115 def get_entropy(self, temperature, verbose=True):
116 """Returns the entropy, in eV/K, in the harmonic approximation
117 at a specified temperature (K)."""
119 self.verbose = verbose
120 write = self._vprint
121 fmt = '%-15s%13.7f eV/K%13.3f eV'
122 write('Entropy components at T = %.2f K:' % temperature)
123 write('=' * 49)
124 write('%15s%13s %13s' % ('', 'S', 'T*S'))
126 S = 0.
128 S_v = self._vibrational_entropy_contribution(temperature)
129 write(fmt % ('S_harm', S_v, S_v * temperature))
130 S += S_v
132 write('-' * 49)
133 write(fmt % ('S', S, S * temperature))
134 write('=' * 49)
135 return S
137 def get_helmholtz_energy(self, temperature, verbose=True):
138 """Returns the Helmholtz free energy, in eV, in the harmonic
139 approximation at a specified temperature (K)."""
141 self.verbose = True
142 write = self._vprint
144 U = self.get_internal_energy(temperature, verbose=verbose)
145 write('')
146 S = self.get_entropy(temperature, verbose=verbose)
147 F = U - temperature * S
149 write('')
150 write('Free energy components at T = %.2f K:' % temperature)
151 write('=' * 23)
152 fmt = '%5s%15.3f eV'
153 write(fmt % ('U', U))
154 write(fmt % ('-T*S', -temperature * S))
155 write('-' * 23)
156 write(fmt % ('F', F))
157 write('=' * 23)
158 return F
161class HinderedThermo(ThermoChem):
162 """Class for calculating thermodynamic properties in the hindered
163 translator and hindered rotor model where all but three degrees of
164 freedom are treated as harmonic vibrations, two are treated as
165 hindered translations, and one is treated as a hindered rotation.
167 Inputs:
169 vib_energies : list
170 a list of all the vibrational energies of the adsorbate (e.g., from
171 ase.vibrations.Vibrations.get_energies). If atoms is not provided,
172 then the number of energies must match the number of degrees of freedom
173 of the adsorbate; i.e., 3*n, where n is the number of atoms. Note
174 that this class does not check that the user has supplied the
175 correct number of energies.
176 Units of energies are eV.
177 trans_barrier_energy : float
178 the translational energy barrier in eV. This is the barrier for an
179 adsorbate to diffuse on the surface.
180 rot_barrier_energy : float
181 the rotational energy barrier in eV. This is the barrier for an
182 adsorbate to rotate about an axis perpendicular to the surface.
183 sitedensity : float
184 density of surface sites in cm^-2
185 rotationalminima : integer
186 the number of equivalent minima for an adsorbate's full rotation.
187 For example, 6 for an adsorbate on an fcc(111) top site
188 potentialenergy : float
189 the potential energy in eV (e.g., from atoms.get_potential_energy)
190 (if potentialenergy is unspecified, then the methods of this class
191 can be interpreted as the energy corrections)
192 mass : float
193 the mass of the adsorbate in amu (if mass is unspecified, then it will
194 be calculated from the atoms class)
195 inertia : float
196 the reduced moment of inertia of the adsorbate in amu*Ang^-2
197 (if inertia is unspecified, then it will be calculated from the
198 atoms class)
199 atoms : an ASE atoms object
200 used to calculate rotational moments of inertia and molecular mass
201 symmetrynumber : integer
202 symmetry number of the adsorbate. This is the number of symmetric arms
203 of the adsorbate and depends upon how it is bound to the surface.
204 For example, propane bound through its end carbon has a symmetry
205 number of 1 but propane bound through its middle carbon has a symmetry
206 number of 2. (if symmetrynumber is unspecified, then the default is 1)
207 ignore_imag_modes : bool
208 If True, any imaginary frequencies present after the 3N-3 cut will not
209 be included in the calculation of the thermochemical properties. If
210 False (default), an error will be raised if imaginary frequencies are
211 present after the 3N-3 cut.
212 """
214 def __init__(self, vib_energies, trans_barrier_energy, rot_barrier_energy,
215 sitedensity, rotationalminima, potentialenergy=0.,
216 mass=None, inertia=None, atoms=None, symmetrynumber=1,
217 ignore_imag_modes=False):
219 self.trans_barrier_energy = trans_barrier_energy * units._e
220 self.rot_barrier_energy = rot_barrier_energy * units._e
221 self.area = 1. / sitedensity / 100.0**2
222 self.rotationalminima = rotationalminima
223 self.potentialenergy = potentialenergy
224 self.atoms = atoms
225 self.symmetry = symmetrynumber
226 self.ignore_imag_modes = ignore_imag_modes
228 # Sort the vibrations
229 vib_energies = list(vib_energies)
230 vib_energies.sort(key=np.abs)
232 # Keep only the relevant vibrational energies (3N-3)
233 if atoms:
234 vib_energies = vib_energies[-(3 * len(atoms) - 3):]
235 else:
236 vib_energies = vib_energies[-(len(vib_energies) - 3):]
238 # Check for imaginary frequencies.
239 vib_energies, n_imag = _clean_vib_energies(
240 vib_energies, ignore_imag_modes=ignore_imag_modes
241 )
242 self.vib_energies = vib_energies
243 self.n_imag = n_imag
245 if (mass or atoms) and (inertia or atoms):
246 if mass:
247 self.mass = mass * units._amu
248 elif atoms:
249 self.mass = np.sum(atoms.get_masses()) * units._amu
250 if inertia:
251 self.inertia = inertia * units._amu / units.m**2
252 elif atoms:
253 self.inertia = (atoms.get_moments_of_inertia()[2] *
254 units._amu / units.m**2)
255 else:
256 raise RuntimeError('Either mass and inertia of the '
257 'adsorbate must be specified or '
258 'atoms must be specified.')
260 # Calculate hindered translational and rotational frequencies
261 self.freq_t = np.sqrt(self.trans_barrier_energy / (2 * self.mass *
262 self.area))
263 self.freq_r = 1. / (2 * np.pi) * np.sqrt(self.rotationalminima**2 *
264 self.rot_barrier_energy /
265 (2 * self.inertia))
267 def get_internal_energy(self, temperature, verbose=True):
268 """Returns the internal energy (including the zero point energy),
269 in eV, in the hindered translator and hindered rotor model at a
270 specified temperature (K)."""
272 from scipy.special import iv
274 self.verbose = verbose
275 write = self._vprint
276 fmt = '%-15s%13.3f eV'
277 write('Internal energy components at T = %.2f K:' % temperature)
278 write('=' * 31)
280 U = 0.
282 write(fmt % ('E_pot', self.potentialenergy))
283 U += self.potentialenergy
285 # Translational Energy
286 T_t = units._k * temperature / (units._hplanck * self.freq_t)
287 R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)
288 dU_t = 2 * (-1. / 2 - 1. / T_t / (2 + 16 * R_t) + R_t / 2 / T_t -
289 R_t / 2 / T_t *
290 iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +
291 1. / T_t / (np.exp(1. / T_t) - 1))
292 dU_t *= units.kB * temperature
293 write(fmt % ('E_trans', dU_t))
294 U += dU_t
296 # Rotational Energy
297 T_r = units._k * temperature / (units._hplanck * self.freq_r)
298 R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)
299 dU_r = (-1. / 2 - 1. / T_r / (2 + 16 * R_r) + R_r / 2 / T_r -
300 R_r / 2 / T_r *
301 iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +
302 1. / T_r / (np.exp(1. / T_r) - 1))
303 dU_r *= units.kB * temperature
304 write(fmt % ('E_rot', dU_r))
305 U += dU_r
307 # Vibrational Energy
308 dU_v = self._vibrational_energy_contribution(temperature)
309 write(fmt % ('E_vib', dU_v))
310 U += dU_v
312 # Zero Point Energy
313 dU_zpe = self.get_zero_point_energy()
314 write(fmt % ('E_ZPE', dU_zpe))
315 U += dU_zpe
317 write('-' * 31)
318 write(fmt % ('U', U))
319 write('=' * 31)
320 return U
322 def get_zero_point_energy(self, verbose=True):
323 """Returns the zero point energy, in eV, in the hindered
324 translator and hindered rotor model"""
326 zpe_t = 2 * (1. / 2 * self.freq_t * units._hplanck / units._e)
327 zpe_r = 1. / 2 * self.freq_r * units._hplanck / units._e
328 zpe_v = self.get_ZPE_correction()
329 zpe = zpe_t + zpe_r + zpe_v
330 return zpe
332 def get_entropy(self, temperature, verbose=True):
333 """Returns the entropy, in eV/K, in the hindered translator
334 and hindered rotor model at a specified temperature (K)."""
336 from scipy.special import iv
338 self.verbose = verbose
339 write = self._vprint
340 fmt = '%-15s%13.7f eV/K%13.3f eV'
341 write('Entropy components at T = %.2f K:' % temperature)
342 write('=' * 49)
343 write('%15s%13s %13s' % ('', 'S', 'T*S'))
345 S = 0.
347 # Translational Entropy
348 T_t = units._k * temperature / (units._hplanck * self.freq_t)
349 R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)
350 S_t = 2 * (-1. / 2 + 1. / 2 * np.log(np.pi * R_t / T_t) -
351 R_t / 2 / T_t *
352 iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +
353 np.log(iv(0, R_t / 2 / T_t)) +
354 1. / T_t / (np.exp(1. / T_t) - 1) -
355 np.log(1 - np.exp(-1. / T_t)))
356 S_t *= units.kB
357 write(fmt % ('S_trans', S_t, S_t * temperature))
358 S += S_t
360 # Rotational Entropy
361 T_r = units._k * temperature / (units._hplanck * self.freq_r)
362 R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)
363 S_r = (-1. / 2 + 1. / 2 * np.log(np.pi * R_r / T_r) -
364 np.log(self.symmetry) -
365 R_r / 2 / T_r * iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +
366 np.log(iv(0, R_r / 2 / T_r)) +
367 1. / T_r / (np.exp(1. / T_r) - 1) -
368 np.log(1 - np.exp(-1. / T_r)))
369 S_r *= units.kB
370 write(fmt % ('S_rot', S_r, S_r * temperature))
371 S += S_r
373 # Vibrational Entropy
374 S_v = self._vibrational_entropy_contribution(temperature)
375 write(fmt % ('S_vib', S_v, S_v * temperature))
376 S += S_v
378 # Concentration Related Entropy
379 N_over_A = np.exp(1. / 3) * (10.0**5 /
380 (units._k * temperature))**(2. / 3)
381 S_c = 1 - np.log(N_over_A) - np.log(self.area)
382 S_c *= units.kB
383 write(fmt % ('S_con', S_c, S_c * temperature))
384 S += S_c
386 write('-' * 49)
387 write(fmt % ('S', S, S * temperature))
388 write('=' * 49)
389 return S
391 def get_helmholtz_energy(self, temperature, verbose=True):
392 """Returns the Helmholtz free energy, in eV, in the hindered
393 translator and hindered rotor model at a specified temperature
394 (K)."""
396 self.verbose = True
397 write = self._vprint
399 U = self.get_internal_energy(temperature, verbose=verbose)
400 write('')
401 S = self.get_entropy(temperature, verbose=verbose)
402 F = U - temperature * S
404 write('')
405 write('Free energy components at T = %.2f K:' % temperature)
406 write('=' * 23)
407 fmt = '%5s%15.3f eV'
408 write(fmt % ('U', U))
409 write(fmt % ('-T*S', -temperature * S))
410 write('-' * 23)
411 write(fmt % ('F', F))
412 write('=' * 23)
413 return F
416class IdealGasThermo(ThermoChem):
417 """Class for calculating thermodynamic properties of a molecule
418 based on statistical mechanical treatments in the ideal gas
419 approximation.
421 Inputs for enthalpy calculations:
423 vib_energies : list
424 a list of the vibrational energies of the molecule (e.g., from
425 ase.vibrations.Vibrations.get_energies). The number of vibrations
426 used is automatically calculated by the geometry and the number of
427 atoms. If more are specified than are needed, then the lowest
428 numbered vibrations are neglected. If either atoms or natoms is
429 unspecified, then uses the entire list. Units are eV.
430 geometry : 'monatomic', 'linear', or 'nonlinear'
431 geometry of the molecule
432 potentialenergy : float
433 the potential energy in eV (e.g., from atoms.get_potential_energy)
434 (if potentialenergy is unspecified, then the methods of this
435 class can be interpreted as the energy corrections)
436 natoms : integer
437 the number of atoms, used along with 'geometry' to determine how
438 many vibrations to use. (Not needed if an atoms object is supplied
439 in 'atoms' or if the user desires the entire list of vibrations
440 to be used.)
442 Extra inputs needed for entropy / free energy calculations:
444 atoms : an ASE atoms object
445 used to calculate rotational moments of inertia and molecular mass
446 symmetrynumber : integer
447 symmetry number of the molecule. See, for example, Table 10.1 and
448 Appendix B of C. Cramer "Essentials of Computational Chemistry",
449 2nd Ed.
450 spin : float
451 the total electronic spin. (0 for molecules in which all electrons
452 are paired, 0.5 for a free radical with a single unpaired electron,
453 1.0 for a triplet with two unpaired electrons, such as O_2.)
454 ignore_imag_modes : bool
455 If True, any imaginary frequencies present after the 3N-5/3N-6 cut
456 will not be included in the calculation of the thermochemical
457 properties. If False (default), a ValueError will be raised if
458 any imaginary frequencies remain after the 3N-5/3N-6 cut.
459 """
461 def __init__(self, vib_energies, geometry, potentialenergy=0.,
462 atoms=None, symmetrynumber=None, spin=None, natoms=None,
463 ignore_imag_modes=False):
464 self.potentialenergy = potentialenergy
465 self.geometry = geometry
466 self.atoms = atoms
467 self.sigma = symmetrynumber
468 self.spin = spin
469 self.ignore_imag_modes = ignore_imag_modes
470 if natoms is None and atoms:
471 natoms = len(atoms)
472 self.natoms = natoms
474 # Sort the vibrations
475 vib_energies = list(vib_energies)
476 vib_energies.sort(key=np.abs)
478 # Cut the vibrations to those needed from the geometry.
479 if natoms:
480 if geometry == 'nonlinear':
481 vib_energies = vib_energies[-(3 * natoms - 6):]
482 elif geometry == 'linear':
483 vib_energies = vib_energies[-(3 * natoms - 5):]
484 elif geometry == 'monatomic':
485 vib_energies = []
486 else:
487 raise ValueError(f"Unsupported geometry: {geometry}")
489 # Check for imaginary frequencies.
490 vib_energies, n_imag = _clean_vib_energies(
491 vib_energies, ignore_imag_modes=ignore_imag_modes
492 )
493 self.vib_energies = vib_energies
494 self.n_imag = n_imag
496 self.referencepressure = 1.0e5 # Pa
498 def get_enthalpy(self, temperature, verbose=True):
499 """Returns the enthalpy, in eV, in the ideal gas approximation
500 at a specified temperature (K)."""
502 self.verbose = verbose
503 write = self._vprint
504 fmt = '%-15s%13.3f eV'
505 write('Enthalpy components at T = %.2f K:' % temperature)
506 write('=' * 31)
508 H = 0.
510 write(fmt % ('E_pot', self.potentialenergy))
511 H += self.potentialenergy
513 zpe = self.get_ZPE_correction()
514 write(fmt % ('E_ZPE', zpe))
515 H += zpe
517 Cv_t = 3. / 2. * units.kB # translational heat capacity (3-d gas)
518 write(fmt % ('Cv_trans (0->T)', Cv_t * temperature))
519 H += Cv_t * temperature
521 if self.geometry == 'nonlinear': # rotational heat capacity
522 Cv_r = 3. / 2. * units.kB
523 elif self.geometry == 'linear':
524 Cv_r = units.kB
525 elif self.geometry == 'monatomic':
526 Cv_r = 0.
527 write(fmt % ('Cv_rot (0->T)', Cv_r * temperature))
528 H += Cv_r * temperature
530 dH_v = self._vibrational_energy_contribution(temperature)
531 write(fmt % ('Cv_vib (0->T)', dH_v))
532 H += dH_v
534 Cp_corr = units.kB * temperature
535 write(fmt % ('(C_v -> C_p)', Cp_corr))
536 H += Cp_corr
538 write('-' * 31)
539 write(fmt % ('H', H))
540 write('=' * 31)
541 return H
543 def get_entropy(self, temperature, pressure, verbose=True):
544 """Returns the entropy, in eV/K, in the ideal gas approximation
545 at a specified temperature (K) and pressure (Pa)."""
547 if self.atoms is None or self.sigma is None or self.spin is None:
548 raise RuntimeError('atoms, symmetrynumber, and spin must be '
549 'specified for entropy and free energy '
550 'calculations.')
551 self.verbose = verbose
552 write = self._vprint
553 fmt = '%-15s%13.7f eV/K%13.3f eV'
554 write('Entropy components at T = %.2f K and P = %.1f Pa:' %
555 (temperature, pressure))
556 write('=' * 49)
557 write('%15s%13s %13s' % ('', 'S', 'T*S'))
559 S = 0.0
561 # Translational entropy (term inside the log is in SI units).
562 mass = sum(self.atoms.get_masses()) * units._amu # kg/molecule
563 S_t = (2 * np.pi * mass * units._k *
564 temperature / units._hplanck**2)**(3.0 / 2)
565 S_t *= units._k * temperature / self.referencepressure
566 S_t = units.kB * (np.log(S_t) + 5.0 / 2.0)
567 write(fmt % ('S_trans (1 bar)', S_t, S_t * temperature))
568 S += S_t
570 # Rotational entropy (term inside the log is in SI units).
571 if self.geometry == 'monatomic':
572 S_r = 0.0
573 elif self.geometry == 'nonlinear':
574 inertias = (self.atoms.get_moments_of_inertia() * units._amu /
575 (10.0**10)**2) # kg m^2
576 S_r = np.sqrt(np.pi * np.prod(inertias)) / self.sigma
577 S_r *= (8.0 * np.pi**2 * units._k * temperature /
578 units._hplanck**2)**(3.0 / 2.0)
579 S_r = units.kB * (np.log(S_r) + 3.0 / 2.0)
580 elif self.geometry == 'linear':
581 inertias = (self.atoms.get_moments_of_inertia() * units._amu /
582 (10.0**10)**2) # kg m^2
583 inertia = max(inertias) # should be two identical and one zero
584 S_r = (8 * np.pi**2 * inertia * units._k * temperature /
585 self.sigma / units._hplanck**2)
586 S_r = units.kB * (np.log(S_r) + 1.)
587 write(fmt % ('S_rot', S_r, S_r * temperature))
588 S += S_r
590 # Electronic entropy.
591 S_e = units.kB * np.log(2 * self.spin + 1)
592 write(fmt % ('S_elec', S_e, S_e * temperature))
593 S += S_e
595 # Vibrational entropy.
596 S_v = self._vibrational_entropy_contribution(temperature)
597 write(fmt % ('S_vib', S_v, S_v * temperature))
598 S += S_v
600 # Pressure correction to translational entropy.
601 S_p = - units.kB * np.log(pressure / self.referencepressure)
602 write(fmt % ('S (1 bar -> P)', S_p, S_p * temperature))
603 S += S_p
605 write('-' * 49)
606 write(fmt % ('S', S, S * temperature))
607 write('=' * 49)
608 return S
610 def get_gibbs_energy(self, temperature, pressure, verbose=True):
611 """Returns the Gibbs free energy, in eV, in the ideal gas
612 approximation at a specified temperature (K) and pressure (Pa)."""
614 self.verbose = verbose
615 write = self._vprint
617 H = self.get_enthalpy(temperature, verbose=verbose)
618 write('')
619 S = self.get_entropy(temperature, pressure, verbose=verbose)
620 G = H - temperature * S
622 write('')
623 write('Free energy components at T = %.2f K and P = %.1f Pa:' %
624 (temperature, pressure))
625 write('=' * 23)
626 fmt = '%5s%15.3f eV'
627 write(fmt % ('H', H))
628 write(fmt % ('-T*S', -temperature * S))
629 write('-' * 23)
630 write(fmt % ('G', G))
631 write('=' * 23)
632 return G
635class CrystalThermo(ThermoChem):
636 """Class for calculating thermodynamic properties of a crystalline
637 solid in the approximation that a lattice of N atoms behaves as a
638 system of 3N independent harmonic oscillators.
640 Inputs:
642 phonon_DOS : list
643 a list of the phonon density of states,
644 where each value represents the phonon DOS at the vibrational energy
645 value of the corresponding index in phonon_energies.
647 phonon_energies : list
648 a list of the range of vibrational energies (hbar*omega) over which
649 the phonon density of states has been evaluated. This list should be
650 the same length as phonon_DOS and integrating phonon_DOS over
651 phonon_energies should yield approximately 3N, where N is the number
652 of atoms per unit cell. If the first element of this list is
653 zero-valued it will be deleted along with the first element of
654 phonon_DOS. Units of vibrational energies are eV.
656 potentialenergy : float
657 the potential energy in eV (e.g., from atoms.get_potential_energy)
658 (if potentialenergy is unspecified, then the methods of this
659 class can be interpreted as the energy corrections)
661 formula_units : int
662 the number of formula units per unit cell. If unspecified, the
663 thermodynamic quantities calculated will be listed on a
664 per-unit-cell basis.
665 """
667 def __init__(self, phonon_DOS, phonon_energies,
668 formula_units=None, potentialenergy=0.):
669 self.phonon_energies = phonon_energies
670 self.phonon_DOS = phonon_DOS
672 if formula_units:
673 self.formula_units = formula_units
674 self.potentialenergy = potentialenergy / formula_units
675 else:
676 self.formula_units = 0
677 self.potentialenergy = potentialenergy
679 def get_internal_energy(self, temperature, verbose=True):
680 """Returns the internal energy, in eV, of crystalline solid
681 at a specified temperature (K)."""
683 self.verbose = verbose
684 write = self._vprint
685 fmt = '%-15s%13.4f eV'
686 if self.formula_units == 0:
687 write('Internal energy components at '
688 'T = %.2f K,\non a per-unit-cell basis:' % temperature)
689 else:
690 write('Internal energy components at '
691 'T = %.2f K,\non a per-formula-unit basis:' % temperature)
692 write('=' * 31)
694 U = 0.
696 omega_e = self.phonon_energies
697 dos_e = self.phonon_DOS
698 if omega_e[0] == 0.:
699 omega_e = np.delete(omega_e, 0)
700 dos_e = np.delete(dos_e, 0)
702 write(fmt % ('E_pot', self.potentialenergy))
703 U += self.potentialenergy
705 zpe_list = omega_e / 2.
706 if self.formula_units == 0:
707 zpe = np.trapz(zpe_list * dos_e, omega_e)
708 else:
709 zpe = np.trapz(zpe_list * dos_e, omega_e) / self.formula_units
710 write(fmt % ('E_ZPE', zpe))
711 U += zpe
713 B = 1. / (units.kB * temperature)
714 E_vib = omega_e / (np.exp(omega_e * B) - 1.)
715 if self.formula_units == 0:
716 E_phonon = np.trapz(E_vib * dos_e, omega_e)
717 else:
718 E_phonon = np.trapz(E_vib * dos_e, omega_e) / self.formula_units
719 write(fmt % ('E_phonon', E_phonon))
720 U += E_phonon
722 write('-' * 31)
723 write(fmt % ('U', U))
724 write('=' * 31)
725 return U
727 def get_entropy(self, temperature, verbose=True):
728 """Returns the entropy, in eV/K, of crystalline solid
729 at a specified temperature (K)."""
731 self.verbose = verbose
732 write = self._vprint
733 fmt = '%-15s%13.7f eV/K%13.4f eV'
734 if self.formula_units == 0:
735 write('Entropy components at '
736 'T = %.2f K,\non a per-unit-cell basis:' % temperature)
737 else:
738 write('Entropy components at '
739 'T = %.2f K,\non a per-formula-unit basis:' % temperature)
740 write('=' * 49)
741 write('%15s%13s %13s' % ('', 'S', 'T*S'))
743 omega_e = self.phonon_energies
744 dos_e = self.phonon_DOS
745 if omega_e[0] == 0.:
746 omega_e = np.delete(omega_e, 0)
747 dos_e = np.delete(dos_e, 0)
749 B = 1. / (units.kB * temperature)
750 S_vib = (omega_e / (temperature * (np.exp(omega_e * B) - 1.)) -
751 units.kB * np.log(1. - np.exp(-omega_e * B)))
752 if self.formula_units == 0:
753 S = np.trapz(S_vib * dos_e, omega_e)
754 else:
755 S = np.trapz(S_vib * dos_e, omega_e) / self.formula_units
757 write('-' * 49)
758 write(fmt % ('S', S, S * temperature))
759 write('=' * 49)
760 return S
762 def get_helmholtz_energy(self, temperature, verbose=True):
763 """Returns the Helmholtz free energy, in eV, of crystalline solid
764 at a specified temperature (K)."""
766 self.verbose = True
767 write = self._vprint
769 U = self.get_internal_energy(temperature, verbose=verbose)
770 write('')
771 S = self.get_entropy(temperature, verbose=verbose)
772 F = U - temperature * S
774 write('')
775 if self.formula_units == 0:
776 write('Helmholtz free energy components at '
777 'T = %.2f K,\non a per-unit-cell basis:' % temperature)
778 else:
779 write('Helmholtz free energy components at '
780 'T = %.2f K,\non a per-formula-unit basis:' % temperature)
781 write('=' * 23)
782 fmt = '%5s%15.4f eV'
783 write(fmt % ('U', U))
784 write(fmt % ('-T*S', -temperature * S))
785 write('-' * 23)
786 write(fmt % ('F', F))
787 write('=' * 23)
788 return F
791def _clean_vib_energies(vib_energies, ignore_imag_modes=False):
792 """Checks for presence of imaginary vibrational modes and removes
793 them if ignore_imag_modes is True. Also removes +0.j from real
794 vibrational energies.
796 Inputs:
798 vib_energies : list
799 a list of the vibrational energies
801 ignore_imag_modes : bool
802 If True, any imaginary frequencies will be removed. If False,
803 (default), an error will be raised if imaginary frequencies are
804 present.
806 Outputs:
808 vib_energies : list
809 a list of the cleaned vibrational energies with imaginary frequencies
810 removed if ignore_imag_modes is True.
812 n_imag : int
813 the number of imaginary frequencies removed
814 """
815 if ignore_imag_modes:
816 n_vib_energies = len(vib_energies)
817 vib_energies = [v for v in vib_energies if np.real(v) > 0]
818 n_imag = n_vib_energies - len(vib_energies)
819 if n_imag > 0:
820 warn(f"{n_imag} imag modes removed", UserWarning)
821 else:
822 if sum(np.iscomplex(vib_energies)):
823 raise ValueError('Imaginary vibrational energies are present.')
824 n_imag = 0
825 vib_energies = np.real(vib_energies) # clear +0.j
827 return vib_energies, n_imag