Coverage for ase / thermochemistry.py: 92.86%
658 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3"""Modules for calculating thermochemical information from computational
4outputs."""
6import os
7import sys
8from abc import ABC, abstractmethod
9from collections.abc import Sequence
10from typing import Literal
11from warnings import catch_warnings, simplefilter, warn
13import numpy as np
14from scipy.integrate import trapezoid
16from ase import Atoms, units
18_GEOMETRY_OPTIONS = Literal['linear', 'nonlinear', 'monatomic']
19_VIB_SELECTION_OPTIONS = Literal['exact', 'all', 'highest', 'abs_highest']
20_FLOAT_OR_FLOATWITHDICT = float | tuple[float, dict[str, float]]
21_FLOATWITHDICT = tuple[float, dict[str, float]]
24def _sum_contributions(contrib_dicts: dict[str, float]) -> float:
25 """Combine a Dict of floats to their sum.
27 Ommits keys starting with an underscore.
28 """
29 return np.sum([value for key, value in contrib_dicts.items()
30 if not key.startswith('_')])
33class AbstractMode(ABC):
34 """Abstract base class for mode objects.
36 Input:
38 energy: float
39 Initialize the mode with its energy in eV.
40 Note: This energy refers to the vibrational energy of the mode.
41 Get it e.g. from ase.vibrations.Vibrations.get_energies.
42 """
44 def __init__(self, energy: float) -> None:
45 self.energy = energy
47 @abstractmethod
48 def get_internal_energy(
49 self,
50 temperature: float,
51 contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
52 raise NotImplementedError
54 @abstractmethod
55 def get_entropy(self, temperature: float,
56 contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
57 raise NotImplementedError
60class HarmonicMode(AbstractMode):
61 """Class for a single harmonic mode."""
63 def get_ZPE_correction(self) -> float:
64 """Returns the zero-point vibrational energy correction in eV."""
65 return 0.5 * self.energy
67 def get_vib_energy_contribution(self, temperature: float) -> float:
68 """Calculates the change in internal energy due to vibrations from
69 0K to the specified temperature for a set of vibrations given in
70 eV and a temperature given in Kelvin.
71 Returns the energy change in eV."""
72 kT = units.kB * temperature
73 return self.energy / (np.exp(self.energy / kT) - 1.)
75 def get_vib_entropy_contribution(self, temperature: float) -> float:
76 """Calculates the entropy due to vibrations given in eV and a
77 temperature given in Kelvin. Returns the entropy in eV/K."""
78 e = self.energy / units.kB / temperature
79 S_v = e / (np.exp(e) - 1.) - np.log(1. - np.exp(-e))
80 S_v *= units.kB
81 return S_v
83 def get_internal_energy(self,
84 temperature: float,
85 contributions: bool
86 ) -> _FLOAT_OR_FLOATWITHDICT:
87 """Returns the internal energy in the harmonic approximation at a
88 specified temperature (K)."""
89 ret = {}
90 ret['ZPE'] = self.get_ZPE_correction()
91 ret['dU_v'] = self.get_vib_energy_contribution(temperature)
92 ret_sum = _sum_contributions(ret)
94 if contributions:
95 return ret_sum, ret
96 else:
97 return ret_sum
99 def get_entropy(self,
100 temperature: float,
101 contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
102 """Returns the entropy in the harmonic approximation at a specified
103 temperature (K).
105 Parameters
106 ----------
107 temperature : float
108 The temperature in Kelvin.
109 contributions : bool
110 If True, return the contributions to the entropy as a dict too.
111 Here it will only contain the vibrational entropy.
113 Returns
114 -------
115 float or Tuple[float, Dict[str, float]]
116 """
117 ret = {}
118 ret['S_v'] = self.get_vib_entropy_contribution(temperature)
119 if contributions:
120 return ret['S_v'], ret
121 else:
122 return ret['S_v']
125class RRHOMode(HarmonicMode):
126 """Class for a single RRHO mode, including Grimme's scaling method
127 based on :doi:`10.1002/chem.201200497` and :doi:`10.1039/D1SC00621E`.
129 Inputs:
131 mean_inertia : float
132 The mean moment of inertia in amu*angstrom^2. Use
133 `np.mean(ase.Atoms.get_moments_of_inertia())` to calculate.
134 tau : float
135 the vibrational energy threshold in :math:`cm^{-1}`, named
136 :math:`\\tau` in :doi:`10.1039/D1SC00621E`.
137 Values close or equal to 0 will result in the standard harmonic
138 approximation. Defaults to :math:`35cm^{-1}`.
139 treat_int_energy : bool
140 Extend the msRRHO treatement to the internal thermal energy.
141 If False, only the entropy contribution as in Grimmmes paper is
142 modified according to the RRHO scheme. If True, the approach of
143 Otlyotov and Minenkov :doi:`10.1002/jcc.27129` is used to modify
144 the internal energy contribution.
145 """
147 def __init__(self, energy: float,
148 mean_inertia: float,
149 tau: float = 35.0,
150 treat_int_energy: bool = False) -> None:
151 if np.iscomplex(energy):
152 raise ValueError(
153 "Imaginary frequencies are not allowed in RRHO mode.")
154 super().__init__(energy)
155 self._mean_inertia = mean_inertia
156 self._tau = tau
157 self._alpha = 4 # from paper 10.1002/chem.201200497
158 self.treat_int_energy = treat_int_energy
160 @property
161 def frequency(self) -> float:
162 return self.energy / units.invcm
164 def _head_gordon_damp(self, freq: float) -> float:
165 """Head-Gordon damping function.
167 Equation 8 from :doi:`10.1002/chem.201200497`
169 Parameters
170 ----------
171 freq : float
172 The frequency in the same unit as tau.
174 Returns
175 -------
176 float
177 """
178 return 1 / (1 + (self._tau / freq)**self._alpha)
180 def _apply_head_gordon_damp(self, freq: float,
181 large_part: float,
182 small_part: float) -> float:
183 """Apply the head-gordon damping scheme to two contributions.
185 Equation 7 from :doi:`10.1002/chem.201200497`
187 Returns the damped sum of the two contributions."""
188 part_one = self._head_gordon_damp(freq) * large_part
189 part_two = (1 - self._head_gordon_damp(freq)) * small_part
190 return part_one + part_two
192 def get_RRHO_entropy_r(self, temperature: float) -> float:
193 """Calculates the rotation of a rigid rotor for low frequency modes.
195 Equation numbering from :doi:`10.1002/chem.201200497`
197 Returns the entropy contribution in eV/K."""
198 kT = units._k * temperature
199 R = units._k * units._Nav # J / K / mol
200 B_av = (self._mean_inertia /
201 (units.kg * units.m**2)) # from amu/A^2 to kg m^2
202 # note, some codes use B_av = 1e-44 as in 10.1002/chem.201200497
203 # eq 4
204 omega = units._c * self.frequency * 1e2 # s^-1
205 mu = units._hplanck / (8 * np.pi**2 * omega) # kg m^2
206 # eq 5
207 mu_prime = (mu * B_av) / (mu + B_av) # kg m^2
208 # eq 6
209 x = np.sqrt(8 * np.pi**3 * mu_prime * kT / (units._hplanck)**2)
210 # filter zeros out and set them to zero
211 if x == 0.0:
212 log_x = 0.0
213 else:
214 log_x = np.log(x)
215 S_r = R * (0.5 + log_x) # J/(Js)^2
216 S_r *= units.J / units._Nav # J/K/mol to eV/K
217 return S_r
219 def get_entropy(self,
220 temperature: float,
221 contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
222 ret = {}
223 ret['_S_vib_v'] = self.get_vib_entropy_contribution(temperature)
224 ret['_S_vib_r'] = self.get_RRHO_entropy_r(temperature)
225 ret['S_vib_damped'] = self._apply_head_gordon_damp(
226 self.frequency, ret['_S_vib_v'], ret['_S_vib_r'])
227 if contributions:
228 return ret['S_vib_damped'], ret
229 else:
230 return ret['S_vib_damped']
232 def get_rrho_internal_energy_v_contribution(
233 self, temperature: float) -> float:
234 """RRHO Vibrational Internal Energy Contribution from
235 :doi:`10.1002/jcc.27129`.
237 Returns the internal energy contribution in eV."""
238 # equation numbering from :doi:`10.1002/jcc.27129`
239 # eq 4
240 # hv = self.energy
241 theta = self.energy / units.kB
242 E_v = 0.5 + 1 / (np.exp(theta / temperature) - 1)
243 E_v *= self.energy # = theta * units.kB
244 return E_v
246 @staticmethod
247 def get_rrho_internal_energy_r_contribution(temperature: float) -> float:
248 """Calculates the rotation of a rigid rotor contribution.
250 Equation numbering from :doi:`10.1002/jcc.27129`
252 Returns the internal energy contribution in eV."""
253 # eq 5
254 R = units._k * units._Nav
255 E_r = R * temperature / 2
256 E_r *= units.J / units._Nav
257 return E_r
259 def get_internal_energy(self,
260 temperature: float,
261 contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
262 """Returns the internal energy in the msRRHO approximation at a
263 specified temperature (K).
265 If self.treat_int_energy is True, the approach of Otlyotov
266 and Minenkov :doi:`10.1002/jcc.27129` is used. Otherwise, the approach
267 of Grimme :doi:`10.1002/chem.201200497` is used.
268 """
269 if self.treat_int_energy:
270 # Otlyotov and Minenkov approach with damping between vibrational
271 # and rotational contributions to the internal energy
272 # Note: The ZPE is not needed here, as the formula in the paper
273 # uses the "bottom of the well" as reference. See
274 # https://gaussian.com/wp-content/uploads/dl/thermo.pdf
275 # for more formulas
276 ret = {}
277 ret['_dU_vib_v'] = self.get_rrho_internal_energy_v_contribution(
278 temperature)
279 ret['_dU_vib_r'] = self.get_rrho_internal_energy_r_contribution(
280 temperature)
281 ret['dU_vib_damped'] = self._apply_head_gordon_damp(
282 self.frequency, ret['_dU_vib_v'], ret['_dU_vib_r'])
283 if contributions:
284 return ret['dU_vib_damped'], ret
285 else:
286 return ret['dU_vib_damped']
287 else:
288 # Grimme uses the Harmonic approach for the internal energy
289 return super().get_internal_energy(temperature, contributions)
292class BaseThermoChem(ABC):
293 """Abstract base class containing common methods used in thermochemistry
294 calculations."""
296 def __init__(self,
297 vib_energies: Sequence[float] | None = None,
298 atoms: Atoms | None = None,
299 modes: Sequence[AbstractMode] | None = None,
300 spin: float | None = None) -> None:
301 if vib_energies is None and modes is None:
302 raise ValueError("Either vib_energies or modes must be provided.")
303 elif modes:
304 self.modes = modes
305 # in the monoatomic case, the vib_energies can be an empty list
306 else:
307 self._vib_energies = vib_energies
308 self.referencepressure = 1.0e5 # Pa
309 if atoms:
310 # check that atoms has no periodic boundary conditions
311 # moments of inertia are wrong otherwise.
312 if atoms.pbc.any():
313 raise ValueError("Atoms object should not have periodic "
314 "boundary conditions for BaseThermoChem.")
315 self.atoms = atoms
316 self.spin = spin
318 @classmethod
319 def from_transition_state(cls, vib_energies: Sequence[complex],
320 *args, **kwargs) -> "BaseThermoChem":
321 """Create a new instance for a transition state.
323 This will work just as the standard constructor, but will remove
324 one imaginary frequency from the given vib_energies first.
325 If there is more than one imaginary frequency, an error will be raised.
327 Returns
328 -------
329 BaseThermoChem instance
330 """
332 if sum(np.iscomplex(vib_energies)):
333 # supress user warning
334 with catch_warnings():
335 simplefilter("ignore")
336 vib_e, n_imag = _clean_vib_energies(vib_energies, True)
337 if n_imag != 1:
338 raise ValueError("Not exactly one imaginary frequency found.")
339 else:
340 raise ValueError("No imaginary frequencies found in vib_energies.")
342 thermo = cls(vib_e, *args, **kwargs)
344 return thermo
346 @staticmethod
347 def combine_contributions(
348 contrib_dicts: Sequence[dict[str, float]]) -> dict[str, float]:
349 """Combine the contributions from multiple modes."""
350 ret: dict[str, float] = {}
351 for contrib_dict in contrib_dicts:
352 for key, value in contrib_dict.items():
353 if key in ret:
354 ret[key] += value
355 else:
356 ret[key] = value
357 return ret
359 def print_contributions(
360 self, contributions: dict[str, float], verbose: bool) -> None:
361 """Print the contributions."""
362 if verbose:
363 fmt = "{:<15s}{:13.3f} eV"
364 for key, value in contributions.items():
365 # subvalues start with _
366 if key.startswith('_'):
367 key.replace('_', '... ')
368 self._vprint(fmt.format(key, value))
370 @abstractmethod
371 def get_internal_energy(self, temperature: float, verbose: bool) -> float:
372 raise NotImplementedError
374 @abstractmethod
375 def get_entropy(self, temperature: float,
376 pressure: float = units.bar,
377 verbose: bool = True) -> float:
378 raise NotImplementedError
380 @property
381 def vib_energies(self) -> np.ndarray:
382 """Get the vibrational energies in eV.
384 For backwards compatibility, this will return the modes' energies if
385 they are available. Otherwise, it will return the stored vib_energies.
386 """
387 # build from modes if available
388 if hasattr(self, 'modes'):
389 return np.array([mode.energy for mode in self.modes])
390 elif hasattr(self, '_vib_energies'):
391 return np.array(self._vib_energies)
392 else:
393 raise AttributeError("No vibrational energies available.")
395 @vib_energies.setter
396 def vib_energies(self, value: Sequence[float]) -> None:
397 """For backwards compatibility, raise a deprecation warning."""
398 warn(
399 "The vib_energies attribute is deprecated and will be removed in a"
400 "future release. Please use the modes attribute instead."
401 "Setting this outside the constructor will not update the modes",
402 DeprecationWarning)
403 self._vib_energies = value
405 def get_ZPE_correction(self) -> float:
406 """Returns the zero-point vibrational energy correction in eV."""
407 return 0.5 * np.sum(self.vib_energies)
409 @staticmethod
410 def get_ideal_translational_energy(temperature: float) -> float:
411 """Returns the translational heat capacity times T in eV.
413 Parameters
414 ----------
415 temperature : float
416 The temperature in Kelvin.
418 Returns
419 -------
420 float
421 """
422 return 3. / 2. * units.kB * \
423 temperature # translational heat capacity (3-d gas)
425 @staticmethod
426 def get_ideal_rotational_energy(geometry: _GEOMETRY_OPTIONS,
427 temperature: float) -> float:
428 """Returns the rotational heat capacity times T in eV.
430 Parameters
431 ----------
432 geometry : str
433 The geometry of the molecule. Options are 'nonlinear',
434 'linear', and 'monatomic'.
435 temperature : float
436 The temperature in Kelvin.
438 Returns
439 -------
440 float
441 """
442 if geometry == 'nonlinear': # rotational heat capacity
443 Cv_r = 3. / 2. * units.kB
444 elif geometry == 'linear':
445 Cv_r = units.kB
446 elif geometry == 'monatomic':
447 Cv_r = 0.
448 else:
449 raise ValueError('Invalid geometry: %s' % geometry)
450 return Cv_r * temperature
452 def get_ideal_trans_entropy(
453 self,
454 atoms: Atoms,
455 temperature: float) -> float:
456 """Returns the translational entropy in eV/K.
458 Parameters
459 ----------
460 atoms : ase.Atoms
461 The atoms object.
462 temperature : float
463 The temperature in Kelvin.
465 Returns
466 -------
467 float
468 """
469 # Translational entropy (term inside the log is in SI units).
470 mass = sum(atoms.get_masses()) * units._amu # kg/molecule
471 S_t = (2 * np.pi * mass * units._k *
472 temperature / units._hplanck**2)**(3.0 / 2)
473 S_t *= units._k * temperature / self.referencepressure
474 S_t = units.kB * (np.log(S_t) + 5.0 / 2.0)
475 return S_t
477 def get_vib_energy_contribution(self, temperature: float) -> float:
478 """Calculates the change in internal energy due to vibrations from
479 0K to the specified temperature for a set of vibrations given in
480 eV and a temperature given in Kelvin.
481 Returns the energy change in eV."""
482 # Leftover for HinderedThermo class, delete once that is converted
483 # to use modes
484 kT = units.kB * temperature
485 dU = 0.
487 for energy in self.vib_energies:
488 dU += energy / (np.exp(energy / kT) - 1.)
489 return dU
491 def get_vib_entropy_contribution(self,
492 temperature: float,
493 return_list: bool = False
494 ) -> float | Sequence[float]:
495 """Calculates the entropy due to vibrations for a set of vibrations
496 given in eV and a temperature given in Kelvin. Returns the entropy
497 in eV/K."""
498 kT = units.kB * temperature
499 S_v = 0.
500 energies = np.array(self.vib_energies)
501 energies /= kT # eV/ eV/K*K
502 S_v = energies / (np.exp(energies) - 1.) - \
503 np.log(1. - np.exp(-energies))
504 S_v *= units.kB
505 if return_list:
506 return S_v
507 else:
508 return np.sum(S_v)
510 def _vprint(self, text):
511 """Print output if verbose flag True."""
512 if self.verbose:
513 sys.stdout.write(text + os.linesep)
515 def get_ideal_entropy(
516 self,
517 temperature: float,
518 translation: bool = False,
519 vibration: bool = False,
520 rotation: bool = False,
521 geometry: _GEOMETRY_OPTIONS | None = None,
522 electronic: bool = False,
523 pressure: float | None = None,
524 symmetrynumber: int | None = None) -> _FLOATWITHDICT:
525 """Returns the entropy, in eV/K and a dict of the contributions.
527 Parameters
528 ----------
529 temperature : float
530 The temperature in Kelvin.
531 translation : bool
532 Include translational entropy.
533 vibration : bool
534 Include vibrational entropy.
535 rotation : bool
536 Include rotational entropy.
537 geometry : str
538 The geometry of the molecule. Options are 'nonlinear',
539 'linear', and 'monatomic'.
540 electronic : bool
541 Include electronic entropy.
542 pressure : float
543 The pressure in Pa. Only needed for the translational entropy.
544 symmetrynumber : int
545 The symmetry number of the molecule. Only needed for linear and
546 nonlinear molecules.
548 Returns
549 -------
550 Tuple of one float and one dict
551 The float is the total entropy in eV/K.
552 The dict contains the contributions to the entropy.
553 """
555 if (geometry in ['linear', 'nonlinear']) and (symmetrynumber is None):
556 raise ValueError(
557 'Symmetry number required for linear and nonlinear molecules.')
559 if not hasattr(self, 'atoms'):
560 raise ValueError(
561 'Atoms object required for ideal entropy calculation.')
563 if electronic and (self.spin is None):
564 raise ValueError(
565 'Spin value required for electronic entropy calculation.')
567 S: float = 0.0
568 ret = {}
570 if translation:
571 S_t = self.get_ideal_trans_entropy(self.atoms, temperature)
572 ret['S_t'] = S_t
573 S += S_t
574 if pressure:
575 # Pressure correction to translational entropy.
576 S_p = - units.kB * np.log(pressure / self.referencepressure)
577 S += S_p
578 ret['S_p'] = S_p
580 if rotation:
581 # Rotational entropy (term inside the log is in SI units).
582 if geometry == 'monatomic':
583 S_r = 0.0
584 elif geometry == 'nonlinear':
585 inertias = (self.atoms.get_moments_of_inertia() * units._amu /
586 1e10**2) # kg m^2
587 S_r = np.sqrt(np.pi * np.prod(inertias)) / symmetrynumber
588 S_r *= (8.0 * np.pi**2 * units._k * temperature /
589 units._hplanck**2)**(3.0 / 2.0)
590 S_r = units.kB * (np.log(S_r) + 3.0 / 2.0)
591 elif geometry == 'linear':
592 inertias = (self.atoms.get_moments_of_inertia() * units._amu /
593 (10.0**10)**2) # kg m^2
594 inertia = max(inertias) # should be two identical and one zero
595 S_r = (8 * np.pi**2 * inertia * units._k * temperature /
596 symmetrynumber / units._hplanck**2)
597 S_r = units.kB * (np.log(S_r) + 1.)
598 else:
599 raise RuntimeError(f"Invalid geometry: {geometry}")
600 S += S_r
601 ret['S_r'] = S_r
603 # Electronic entropy.
604 if electronic:
605 assert self.spin is not None # for mypy, error is raised above
606 S_e = units.kB * np.log(2 * self.spin + 1)
607 S += S_e
608 ret['S_e'] = S_e
610 # Vibrational entropy
611 if vibration:
612 S_v = self.get_vib_entropy_contribution(temperature,
613 return_list=False)
614 assert isinstance(S_v, float) # make mypy happy
615 S += S_v
616 ret['S_v'] = S_v
618 return S, ret
621class HarmonicThermo(BaseThermoChem):
622 """Class for calculating thermodynamic properties in the approximation
623 that all degrees of freedom are treated harmonically. Often used for
624 adsorbates.
626 Note: This class does not include the translational and rotational
627 contributions to the entropy by default. Use the get_ideal_entropy method
628 for that and add them manually.
630 Inputs:
632 vib_energies : list
633 a list of the harmonic energies of the adsorbate (e.g., from
634 ase.vibrations.Vibrations.get_energies). The number of
635 energies should match the number of degrees of freedom of the
636 adsorbate; i.e., 3*n, where n is the number of atoms. Note that
637 this class does not check that the user has supplied the correct
638 number of energies. Units of energies are eV. Note, that if 'modes'
639 are given, these energies are not used.
640 potentialenergy : float
641 the potential energy in eV (e.g., from atoms.get_potential_energy)
642 (if potentialenergy is unspecified, then the methods of this
643 class can be interpreted as the energy corrections)
644 ignore_imag_modes : bool
645 If 'True', any imaginary frequencies will be removed in the
646 calculation of the thermochemical properties.
647 If 'False' (default), an error will be raised if any imaginary
648 frequencies are present.
649 modes : list of AbstractMode
650 A list of mode objects. If not provided, :class:`HarmonicMode` objects
651 will be created from the vib_energies. This is useful if you want to
652 replace individual modes with non-harmonic modes.
653 """
655 def __init__(self, vib_energies: Sequence[complex],
656 potentialenergy: float = 0.,
657 ignore_imag_modes: bool = False,
658 modes: Sequence[AbstractMode] | None = None) -> None:
660 # Check for imaginary frequencies.
661 vib_energies, n_imag = _clean_vib_energies(
662 vib_energies, ignore_imag_modes)
663 if modes is None:
664 modes = [HarmonicMode(energy) for energy in vib_energies]
665 super().__init__(modes=modes)
667 self.n_imag = n_imag
669 self.potentialenergy = potentialenergy
671 def get_internal_energy(self, temperature: float,
672 verbose: bool = True) -> float:
673 """Returns the internal energy, in eV, in the harmonic approximation
674 at a specified temperature (K)."""
676 self.verbose = verbose
677 vprint = self._vprint
678 fmt = '%-15s%13.3f eV'
679 vprint('Internal energy components at T = %.2f K:' % temperature)
680 vprint('=' * 31)
682 vprint(fmt % ('E_pot', self.potentialenergy))
684 U, contribs = zip(*[mode.get_internal_energy(
685 temperature, contributions=True) for mode in self.modes])
686 U = np.sum(U)
687 U += self.potentialenergy
689 self.print_contributions(self.combine_contributions(contribs), verbose)
690 vprint('-' * 31)
691 vprint(fmt % ('U', U))
692 vprint('=' * 31)
694 return U
696 def get_entropy(self, temperature: float,
697 pressure: float = units.bar,
698 verbose: bool = True) -> float:
699 """Returns the entropy, in eV/K at a specified temperature (K).
701 Note: This does not include the translational and rotational
702 contributions to the entropy. Use the get_ideal_entropy method
703 for that.
705 Parameters
706 ----------
707 temperature : float
708 The temperature in Kelvin.
709 pressure : float
710 Not used, but kept for compatibility with other classes.
711 verbose : bool
712 If True, print the contributions to the entropy.
714 Returns
715 -------
716 float
717 """
719 self.verbose = verbose
720 vprint = self._vprint
721 fmt = '%-15s%13.7f eV/K%13.3f eV'
722 vprint('Entropy components at T = %.2f K:' % temperature)
723 vprint('=' * 49)
724 vprint('%15s%13s %13s' % ('', 'S', 'T*S'))
726 S, contribs = zip(*[mode.get_entropy(
727 temperature, contributions=True) for mode in self.modes])
728 S = np.sum(S)
730 self.print_contributions(self.combine_contributions(contribs), verbose)
731 vprint('-' * 49)
732 vprint(fmt % ('S', S, S * temperature))
733 vprint('=' * 49)
735 return S
737 def get_helmholtz_energy(self, temperature: float,
738 verbose: bool = True) -> float:
739 """Returns the Helmholtz free energy, in eV, in the harmonic
740 approximation at a specified temperature (K)."""
742 self.verbose = True
743 vprint = self._vprint
745 U = self.get_internal_energy(temperature, verbose=verbose)
746 vprint('')
747 S = self.get_entropy(temperature, verbose=verbose)
748 F = U - temperature * S
750 vprint('')
751 vprint('Free energy components at T = %.2f K:' % temperature)
752 vprint('=' * 23)
753 fmt = '%5s%15.3f eV'
754 vprint(fmt % ('U', U))
755 vprint(fmt % ('-T*S', -temperature * S))
756 vprint('-' * 23)
757 vprint(fmt % ('F', F))
758 vprint('=' * 23)
759 return F
762class QuasiHarmonicThermo(HarmonicThermo):
763 """Subclass of :class:`HarmonicThermo`, including the quasi-harmonic
764 approximation of Cramer, Truhlar and coworkers :doi:`10.1021/jp205508z`.
766 Note: This class not include the translational and rotational
767 contributions to the entropy by default. Use the get_ideal_entropy method
768 for that and add them manually.
770 Inputs:
772 vib_energies : list
773 a list of the energies of the vibrations (e.g., from
774 ase.vibrations.Vibrations.get_energies). The number of
775 energies should match the number of degrees of freedom of the
776 adsorbate; i.e., 3*n, where n is the number of atoms. Note that
777 this class does not check that the user has supplied the correct
778 number of energies. Units of energies are eV.
779 potentialenergy : float
780 the potential energy in eV (e.g., from atoms.get_potential_energy)
781 (if potentialenergy is unspecified, then the methods of this
782 class can be interpreted as the energy corrections)
783 ignore_imag_modes : bool
784 If 'True', any imaginary frequencies will be removed in the
785 calculation of the thermochemical properties.
786 If 'False' (default), an error will be raised if any imaginary
787 frequencies are present.
788 modes : list of AbstractMode
789 A list of mode objects. If not provided, :class:`HarmonicMode` objects
790 will be created from the raised vib_energies. This is useful if you want
791 to replace individual modes with non-harmonic modes.
792 raise_to : float
793 The value to which all frequencies smaller than this value will be
794 raised. Unit is eV. Defaults to :math:`100cm^{-1} = 0.012398 eV`.
796 """
798 @staticmethod
799 def _raise(input: Sequence[float], raise_to: float) -> Sequence[float]:
800 return [raise_to if x < raise_to else x for x in input]
802 def __init__(self, vib_energies: Sequence[complex],
803 potentialenergy: float = 0.,
804 ignore_imag_modes: bool = False,
805 modes: Sequence[AbstractMode] | None = None,
806 raise_to: float = 100 * units.invcm) -> None:
808 # Check for imaginary frequencies.
809 vib_energies, _ = _clean_vib_energies(
810 vib_energies, ignore_imag_modes)
811 # raise the low frequencies to a certain value
812 vib_energies = self._raise(vib_energies, raise_to)
813 if modes is None:
814 modes = [HarmonicMode(energy) for energy in vib_energies]
815 super().__init__(vib_energies,
816 potentialenergy=potentialenergy,
817 ignore_imag_modes=ignore_imag_modes,
818 modes=modes)
821class MSRRHOThermo(QuasiHarmonicThermo):
822 """Subclass of :class:`QuasiHarmonicThermo`,
823 including Grimme's scaling method based on
824 :doi:`10.1002/chem.201200497` and :doi:`10.1039/D1SC00621E`.
826 Note: This class not include the translational and rotational
827 contributions to the entropy by default. Use the get_ideal_entropy method
828 for that and add them manually.
830 We enforce treating imaginary modes as Grimme suggests (converting
831 them to real by multiplying them with :math:`-i`). So make sure to check
832 your input energies.
834 Inputs:
836 vib_energies : list
837 a list of the energies of the vibrations (e.g., from
838 ase.vibrations.Vibrations.get_energies). The number of
839 energies should match the number of degrees of freedom of the
840 atoms object; i.e., 3*n, where n is the number of atoms. Note that
841 this class does not check that the user has supplied the correct
842 number of energies. Units of energies are eV. Note, that if 'modes'
843 are given, these energies are not used.
844 atoms: an ASE atoms object
845 used to calculate rotational moments of inertia and molecular mass.
846 Should not contain periodic boundary conditions.
847 tau : float
848 the vibrational energy threshold in :math:`cm^{-1}`, namcomplexed
849 :math:`\\tau` in :doi:`10.1039/D1SC00621E`.
850 Values close or equal to 0 will result in the standard harmonic
851 approximation. Defaults to :math:`35cm^{-1}`.
852 nu_scal : float
853 Linear scaling factor for the vibrational frequencies. Named
854 :math:`\\nu_{scal}` in :doi:`10.1039/D1SC00621E`.
855 Defaults to 1.0, check the `Truhlar group database
856 <https://comp.chem.umn.edu/freqscale/index.html>`_
857 for values corresponding to your level of theory.
858 Note that for :math:`\\nu_{scal}=1.0` this method is equivalent to
859 the quasi-RRHO method in :doi:`10.1002/chem.201200497`.
860 treat_int_energy : bool
861 Extend the msRRHO treatement to the internal energy. If False, only
862 the entropy contribution as in Grimmmes paper is considered.
863 If true, the approach of Otlyotov and Minenkov
864 :doi:`10.1002/jcc.27129` is used.
865 modes : list of AbstractMode
866 A list of mode objects. If not provided, :class:`RRHOMode` objects will
867 be created from the raised vib_energies. This is useful if you want to
868 replace individual modes with non-harmonic modes.
870 """
872 def __init__(self, vib_energies: Sequence[complex], atoms: Atoms,
873 potentialenergy: float = 0.,
874 tau: float = 35., nu_scal: float = 1.0,
875 treat_int_energy: bool = False,
876 modes: Sequence[AbstractMode] | None = None) -> None:
878 # check that atoms has no periodic boundary conditions
879 # moments of inertia are wrong otherwise.
880 if atoms.pbc.any():
881 raise ValueError("Atoms object should not have periodic "
882 "boundary conditions for MSRRHOThermo.")
883 inertia = np.mean(atoms.get_moments_of_inertia())
884 self.atoms = atoms
886 # clean the energies by converting imaginary to real
887 # i.e. multiply with -i
888 vib_e: Sequence[float] # make mypy happy
889 vib_e = [np.imag(v) if np.iscomplex(v)
890 else np.real(v) for v in vib_energies]
892 self.nu_scal = nu_scal
893 # scale the frequencies (i.e. energies) before passing them on
894 vib_e_scaled = [float(v) for v in np.multiply(vib_e, nu_scal)]
896 if modes is None:
897 modes = [RRHOMode(energy, inertia,
898 tau=tau,
899 treat_int_energy=treat_int_energy
900 ) for energy in vib_e_scaled]
902 super().__init__(vib_e_scaled,
903 potentialenergy=potentialenergy,
904 ignore_imag_modes=False,
905 modes=modes,
906 raise_to=0.0)
907 self.treat_int_energy = treat_int_energy
910class HinderedThermo(BaseThermoChem):
911 """Class for calculating thermodynamic properties in the hindered
912 translator and hindered rotor model where all but three degrees of
913 freedom are treated as harmonic vibrations, two are treated as
914 hindered translations, and one is treated as a hindered rotation.
916 Inputs:
918 vib_energies : list
919 a list of all the vibrational energies of the adsorbate (e.g., from
920 ase.vibrations.Vibrations.get_energies). If atoms is not provided,
921 then the number of energies must match the number of degrees of freedom
922 of the adsorbate; i.e., 3*n, where n is the number of atoms. Note
923 that this class does not check that the user has supplied the
924 correct number of energies.
925 Units of energies are eV.
926 trans_barrier_energy : float
927 the translational energy barrier in eV. This is the barrier for an
928 adsorbate to diffuse on the surface.
929 rot_barrier_energy : float
930 the rotational energy barrier in eV. This is the barrier for an
931 adsorbate to rotate about an axis perpendicular to the surface.
932 sitedensity : float
933 density of surface sites in cm^-2
934 rotationalminima : integer
935 the number of equivalent minima for an adsorbate's full rotation.
936 For example, 6 for an adsorbate on an fcc(111) top site
937 potentialenergy : float
938 the potential energy in eV (e.g., from atoms.get_potential_energy)
939 (if potentialenergy is unspecified, then the methods of this class
940 can be interpreted as the energy corrections)
941 mass : float
942 the mass of the adsorbate in amu (if mass is unspecified, then it will
943 be calculated from the atoms class)
944 inertia : float
945 the reduced moment of inertia of the adsorbate in amu*Ang^-2
946 (if inertia is unspecified, then it will be calculated from the
947 atoms class)
948 atoms : an ASE atoms object
949 used to calculate rotational moments of inertia and molecular mass
950 symmetrynumber : integer
951 symmetry number of the adsorbate. This is the number of symmetric arms
952 of the adsorbate and depends upon how it is bound to the surface.
953 For example, propane bound through its end carbon has a symmetry
954 number of 1 but propane bound through its middle carbon has a symmetry
955 number of 2. (if symmetrynumber is unspecified, then the default is 1)
956 ignore_imag_modes : bool
957 If 'True', any imaginary frequencies present after the 3N-3 cut will
958 be removed in the calculation of the thermochemical properties.
959 If 'False' (default), an error will be raised if imaginary frequencies
960 are present after the 3N-3 cut.
961 """
963 def __init__(self, vib_energies, trans_barrier_energy, rot_barrier_energy,
964 sitedensity, rotationalminima, potentialenergy=0.,
965 mass=None, inertia=None, atoms=None, symmetrynumber=1,
966 ignore_imag_modes: bool = False) -> None:
968 self.trans_barrier_energy = trans_barrier_energy * units._e
969 self.rot_barrier_energy = rot_barrier_energy * units._e
970 self.area = 1. / sitedensity / 100.0**2
971 self.rotationalminima = rotationalminima
972 self.potentialenergy = potentialenergy
973 self.atoms = atoms
974 self.symmetry = symmetrynumber
976 # Sort the vibrations
977 vib_energies = list(vib_energies)
978 vib_energies.sort(key=np.abs)
980 # Keep only the relevant vibrational energies (3N-3)
981 if atoms:
982 vib_energies = vib_energies[-(3 * len(atoms) - 3):]
983 else:
984 vib_energies = vib_energies[-(len(vib_energies) - 3):]
986 # Check for imaginary frequencies.
987 vib_energies, n_imag = _clean_vib_energies(
988 vib_energies, ignore_imag_modes)
989 super().__init__(vib_energies=vib_energies)
990 self.n_imag = n_imag
992 if (mass or atoms) and (inertia or atoms):
993 if mass:
994 self.mass = mass * units._amu
995 elif atoms:
996 self.mass = np.sum(atoms.get_masses()) * units._amu
997 if inertia:
998 self.inertia = inertia * units._amu / units.m**2
999 elif atoms:
1000 # This is tricky because the moments of inertia
1001 # only work for non-periodic systems.
1002 # Check it here and raise an error if periodic
1003 if atoms.get_pbc().any():
1004 raise RuntimeError('Atoms object with periodic '
1005 'boundary conditions cannot be '
1006 'used to calculate moments of '
1007 'inertia. Please provide '
1008 'a non-periodic Atoms object.')
1009 self.inertia = (atoms.get_moments_of_inertia()[2] *
1010 units._amu / units.m**2)
1011 else:
1012 raise RuntimeError('Either mass and inertia of the '
1013 'adsorbate must be specified or '
1014 'atoms must be specified.')
1016 # Calculate hindered translational and rotational frequencies
1017 self.freq_t = np.sqrt(self.trans_barrier_energy /
1018 (2 * self.mass * self.area))
1019 self.freq_r = 1. / (2 * np.pi) * np.sqrt(self.rotationalminima**2 *
1020 self.rot_barrier_energy /
1021 (2 * self.inertia))
1023 def get_internal_energy(self, temperature, verbose=True):
1024 """Returns the internal energy (including the zero point energy),
1025 in eV, in the hindered translator and hindered rotor model at a
1026 specified temperature (K)."""
1028 from scipy.special import iv
1030 self.verbose = verbose
1031 vprint = self._vprint
1032 fmt = '%-15s%13.3f eV'
1033 vprint('Internal energy components at T = %.2f K:' % temperature)
1034 vprint('=' * 31)
1036 U = 0.
1038 vprint(fmt % ('E_pot', self.potentialenergy))
1039 U += self.potentialenergy
1041 # Translational Energy
1042 T_t = units._k * temperature / (units._hplanck * self.freq_t)
1043 R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)
1044 dU_t = 2 * (-1. / 2 - 1. / T_t / (2 + 16 * R_t) + R_t / 2 / T_t -
1045 R_t / 2 / T_t *
1046 iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +
1047 1. / T_t / (np.exp(1. / T_t) - 1))
1048 dU_t *= units.kB * temperature
1049 vprint(fmt % ('E_trans', dU_t))
1050 U += dU_t
1052 # Rotational Energy
1053 T_r = units._k * temperature / (units._hplanck * self.freq_r)
1054 R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)
1055 dU_r = (-1. / 2 - 1. / T_r / (2 + 16 * R_r) + R_r / 2 / T_r -
1056 R_r / 2 / T_r *
1057 iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +
1058 1. / T_r / (np.exp(1. / T_r) - 1))
1059 dU_r *= units.kB * temperature
1060 vprint(fmt % ('E_rot', dU_r))
1061 U += dU_r
1063 # Vibrational Energy
1064 dU_v = self.get_vib_energy_contribution(temperature)
1065 vprint(fmt % ('E_vib', dU_v))
1066 U += dU_v
1068 # Zero Point Energy
1069 dU_zpe = self.get_zero_point_energy()
1070 vprint(fmt % ('E_ZPE', dU_zpe))
1071 U += dU_zpe
1073 vprint('-' * 31)
1074 vprint(fmt % ('U', U))
1075 vprint('=' * 31)
1076 return U
1078 def get_zero_point_energy(self, verbose=True):
1079 """Returns the zero point energy, in eV, in the hindered
1080 translator and hindered rotor model"""
1082 zpe_t = 2 * (1. / 2 * self.freq_t * units._hplanck / units._e)
1083 zpe_r = 1. / 2 * self.freq_r * units._hplanck / units._e
1084 zpe_v = self.get_ZPE_correction()
1085 zpe = zpe_t + zpe_r + zpe_v
1086 return zpe
1088 def get_entropy(self, temperature,
1089 pressure=units.bar,
1090 verbose=True):
1091 """Returns the entropy, in eV/K, in the hindered translator
1092 and hindered rotor model at a specified temperature (K)."""
1094 from scipy.special import iv
1096 self.verbose = verbose
1097 vrpint = self._vprint
1098 fmt = '%-15s%13.7f eV/K%13.3f eV'
1099 vrpint('Entropy components at T = %.2f K:' % temperature)
1100 vrpint('=' * 49)
1101 vrpint('%15s%13s %13s' % ('', 'S', 'T*S'))
1103 S = 0.
1105 # Translational Entropy
1106 T_t = units._k * temperature / (units._hplanck * self.freq_t)
1107 R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)
1108 S_t = 2 * (-1. / 2 + 1. / 2 * np.log(np.pi * R_t / T_t) -
1109 R_t / 2 / T_t *
1110 iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +
1111 np.log(iv(0, R_t / 2 / T_t)) +
1112 1. / T_t / (np.exp(1. / T_t) - 1) -
1113 np.log(1 - np.exp(-1. / T_t)))
1114 S_t *= units.kB
1115 vrpint(fmt % ('S_trans', S_t, S_t * temperature))
1116 S += S_t
1118 # Rotational Entropy
1119 T_r = units._k * temperature / (units._hplanck * self.freq_r)
1120 R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)
1121 S_r = (-1. / 2 + 1. / 2 * np.log(np.pi * R_r / T_r) -
1122 np.log(self.symmetry) -
1123 R_r / 2 / T_r * iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +
1124 np.log(iv(0, R_r / 2 / T_r)) +
1125 1. / T_r / (np.exp(1. / T_r) - 1) -
1126 np.log(1 - np.exp(-1. / T_r)))
1127 S_r *= units.kB
1128 vrpint(fmt % ('S_rot', S_r, S_r * temperature))
1129 S += S_r
1131 # Vibrational Entropy
1132 S_v = self.get_vib_entropy_contribution(temperature)
1133 vrpint(fmt % ('S_vib', S_v, S_v * temperature))
1134 S += S_v
1136 # Concentration Related Entropy
1137 N_over_A = np.exp(1. / 3) * (10.0**5 /
1138 (units._k * temperature))**(2. / 3)
1139 S_c = 1 - np.log(N_over_A) - np.log(self.area)
1140 S_c *= units.kB
1141 vrpint(fmt % ('S_con', S_c, S_c * temperature))
1142 S += S_c
1144 vrpint('-' * 49)
1145 vrpint(fmt % ('S', S, S * temperature))
1146 vrpint('=' * 49)
1147 return S
1149 def get_helmholtz_energy(self, temperature, verbose=True):
1150 """Returns the Helmholtz free energy, in eV, in the hindered
1151 translator and hindered rotor model at a specified temperature
1152 (K)."""
1154 self.verbose = True
1155 vprint = self._vprint
1157 U = self.get_internal_energy(temperature, verbose=verbose)
1158 vprint('')
1159 S = self.get_entropy(temperature, verbose=verbose)
1160 F = U - temperature * S
1162 vprint('')
1163 vprint('Free energy components at T = %.2f K:' % temperature)
1164 vprint('=' * 23)
1165 fmt = '%5s%15.3f eV'
1166 vprint(fmt % ('U', U))
1167 vprint(fmt % ('-T*S', -temperature * S))
1168 vprint('-' * 23)
1169 vprint(fmt % ('F', F))
1170 vprint('=' * 23)
1171 return F
1174class IdealGasThermo(BaseThermoChem):
1175 """Class for calculating thermodynamic properties of a molecule
1176 based on statistical mechanical treatments in the ideal gas
1177 approximation.
1179 Inputs for enthalpy calculations:
1181 vib_energies : list
1182 a list of the vibrational energies of the molecule (e.g., from
1183 ase.vibrations.Vibrations.get_energies). The number of expected
1184 vibrations is calculated by the geometry and the number of atoms.
1185 By default, the lowest energies are cut to enforce the expected
1186 length. By setting vib_selection, this behaviour can be changed.
1187 If either atoms and natoms is unspecified, the full list is used.
1188 Units are eV.
1189 geometry : 'monatomic', 'linear', or 'nonlinear'
1190 geometry of the molecule
1191 potentialenergy : float
1192 the potential energy in eV (e.g., from atoms.get_potential_energy)
1193 (if potentialenergy is unspecified, then the methods of this
1194 class can be interpreted as the energy corrections)
1195 natoms : integer
1196 the number of atoms, used along with 'geometry' to determine how
1197 many vibrations to use. (Not needed if an atoms object is supplied
1198 in 'atoms' or if the user desires the entire list of vibrations
1199 to be used.)
1200 vib_selection : 'exact', 'highest' (default), 'abs_highest', 'all'
1201 selection of input vibrational energies considered to be true
1202 vibrations (excluding translations and rotations) implied by the
1203 geometry and number of atoms. Only applied if number of atoms
1204 is provided, either with natoms or atoms.
1205 - 'exact': no selection; the number of input vibrational energies must
1206 be equal to the number of true vibrations
1207 - 'highest': select vibrational energies whose square is the highest
1208 , i.e. large real energies followed by small real energies, small
1209 imaginary energies and large imaginary energies. Will not catch
1210 unrealistically large imaginary frequencies, resulting from e.g.
1211 unrelaxed molecules. This is the default.
1212 - 'abs_highest': select vibrational energies whose absolute number are
1213 highest, real or imaginary. Will fail intentionally for
1214 unrealistically large imaginary frequencies, but also unintentionally
1215 for small real energies, corresponding e.g. to frustrated rotations
1216 in a larger molecule.
1217 - 'all': Use all input vibrational energies without checking the
1218 number of them.
1219 ignore_imag_modes : bool
1220 If 'True', any imaginary frequencies present after the 3N-3 cut will
1221 be removed in the calculation of the thermochemical properties.
1222 If 'False' (default), an error will be raised if imaginary frequencies
1223 are present after the 3N-3 cut.
1225 Extra inputs needed for entropy / free energy calculations:
1227 atoms : an ASE atoms object
1228 used to calculate rotational moments of inertia and molecular mass
1229 symmetrynumber : integer
1230 symmetry number of the molecule. See, for example, Table 10.1 and
1231 Appendix B of C. Cramer "Essentials of Computational Chemistry",
1232 2nd Ed.
1233 spin : float
1234 the total electronic spin. (0 for molecules in which all electrons
1235 are paired, 0.5 for a free radical with a single unpaired electron,
1236 1.0 for a triplet with two unpaired electrons, such as O_2.)
1238 """
1240 def __init__(self, vib_energies: Sequence[complex],
1241 geometry: _GEOMETRY_OPTIONS,
1242 potentialenergy: float = 0.,
1243 atoms: Atoms | None = None,
1244 symmetrynumber: int | None = None,
1245 spin: float | None = None,
1246 natoms: int | None = None,
1247 vib_selection: _VIB_SELECTION_OPTIONS | None = 'highest',
1248 ignore_imag_modes: bool = False,
1249 modes: Sequence[AbstractMode] | None = None) -> None:
1250 self.potentialenergy = potentialenergy
1251 self.geometry = geometry
1252 self.sigma = symmetrynumber
1253 if natoms is None and atoms:
1254 natoms = len(atoms)
1256 # Preliminary sorting vibrations by square
1257 vib_energies = list(vib_energies)
1258 vib_energies.sort(key=lambda f: (f ** 2).real)
1260 if natoms and vib_selection != 'all':
1261 # Determine number of true vibrations from geometry
1262 if geometry == 'nonlinear':
1263 num_vibs = 3 * natoms - 6
1264 elif geometry == 'linear':
1265 num_vibs = 3 * natoms - 5
1266 elif geometry == 'monatomic':
1267 num_vibs = 0
1268 else:
1269 raise ValueError(f"Unsupported geometry: {geometry}")
1270 if num_vibs < 0:
1271 raise ValueError(
1272 f"Number of atoms ({natoms}) too small for " +
1273 f"geometry '{geometry}'."
1274 )
1276 if vib_selection == 'exact':
1277 # Demand the expected number of true vibrations
1278 if len(vib_energies) != num_vibs:
1279 raise ValueError(
1280 f"{num_vibs} vibrational energies expected, " +
1281 f"{len(vib_energies)} received.\n" +
1282 "To select true vibrational energies automatically," +
1283 " set vib_selection."
1284 )
1286 elif vib_selection in ('highest', 'abs_highest'):
1287 # Cut vibrations to get the expected number of true vibrations
1288 if num_vibs == 0:
1289 vib_energies = []
1290 else:
1291 if vib_selection == 'abs_highest':
1292 vib_energies.sort(key=np.abs)
1293 vib_energies = vib_energies[-num_vibs:]
1294 if len(vib_energies) != num_vibs:
1295 raise ValueError(
1296 f"Too few vibration modes ({len(vib_energies)}) " +
1297 "after selection. {num_vibs} expected"
1298 )
1299 else:
1300 raise ValueError(
1301 f"(Unsupported vib_selection: {vib_selection})"
1302 )
1303 elif not natoms and not atoms and vib_selection != 'all':
1304 raise ValueError(
1305 "Either natoms or atoms must be specified to select " +
1306 "true vibrational energies. Else, set vib_selection to 'all'."
1307 )
1309 # Check for imaginary frequencies.
1310 vib_energies, n_imag = _clean_vib_energies(
1311 vib_energies, ignore_imag_modes)
1312 super().__init__(vib_energies=vib_energies,
1313 atoms=atoms,
1314 spin=spin)
1315 self.n_imag = n_imag
1317 def get_internal_energy(self, temperature: float,
1318 verbose: bool = True) -> float:
1319 """Returns the internal energy, in eV, in the ideal gas approximation
1320 at a specified temperature (K)."""
1322 self.verbose = verbose
1323 vprint = self._vprint
1324 fmt = '%-15s%13.3f eV'
1325 vprint('Enthalpy components at T = %.2f K:' % temperature)
1326 vprint('=' * 31)
1328 U = 0.
1330 vprint(fmt % ('E_pot', self.potentialenergy))
1331 U += self.potentialenergy
1333 zpe = self.get_ZPE_correction()
1334 vprint(fmt % ('E_ZPE', zpe))
1335 U += zpe
1337 Cv_tT = self.get_ideal_translational_energy(temperature)
1338 vprint(fmt % ('Cv_trans (0->T)', Cv_tT))
1339 U += Cv_tT
1341 Cv_rT = self.get_ideal_rotational_energy(self.geometry, temperature)
1342 vprint(fmt % ('Cv_rot (0->T)', Cv_rT))
1343 U += Cv_rT
1345 dU_v = self.get_vib_energy_contribution(temperature)
1346 vprint(fmt % ('Cv_vib (0->T)', dU_v))
1347 U += dU_v
1349 vprint('-' * 31)
1350 vprint(fmt % ('U', U))
1351 vprint('=' * 31)
1352 return U
1354 def get_enthalpy(self, temperature: float,
1355 verbose: bool = True) -> float:
1356 """Returns the enthalpy, in eV, in the ideal gas approximation
1357 at a specified temperature (K)."""
1359 self.verbose = verbose
1360 vprint = self._vprint
1361 fmt = '%-15s%13.3f eV'
1362 vprint('Enthalpy components at T = %.2f K:' % temperature)
1363 vprint('=' * 31)
1365 H = 0.
1366 H += self.get_internal_energy(temperature, verbose=verbose)
1368 Cp_corr = units.kB * temperature
1369 vprint(fmt % ('(C_v -> C_p)', Cp_corr))
1370 H += Cp_corr
1372 vprint('-' * 31)
1373 vprint(fmt % ('H', H))
1374 vprint('=' * 31)
1375 return H
1377 def get_entropy(self, temperature: float,
1378 pressure: float = units.bar,
1379 verbose: bool = True) -> float:
1380 """Returns the entropy, in eV/K, in the ideal gas approximation
1381 at a specified temperature (K) and pressure (Pa)."""
1383 if self.atoms is None or self.sigma is None or self.spin is None:
1384 raise RuntimeError('atoms, symmetrynumber, and spin must be '
1385 'specified for entropy and free energy '
1386 'calculations.')
1387 self.verbose = verbose
1388 vprint = self._vprint
1389 fmt = '%-15s%13.7f eV/K%13.3f eV'
1390 vprint(f'Entropy components at T = {temperature:.2f} K and'
1391 f' P = {pressure:.1f} Pa:')
1392 vprint('=' * 49)
1393 vprint('{"":15s}{"S":13s} {"T*S:13s}')
1394 S, S_dict = self.get_ideal_entropy(temperature,
1395 translation=True,
1396 vibration=True,
1397 rotation=True,
1398 geometry=self.geometry,
1399 electronic=True,
1400 pressure=pressure,
1401 symmetrynumber=self.sigma)
1403 vprint(
1404 fmt %
1405 ('S_trans (1 bar)',
1406 S_dict['S_t'],
1407 S_dict['S_t'] *
1408 temperature))
1409 vprint(fmt % ('S_rot', S_dict['S_r'], S_dict['S_r'] * temperature))
1410 vprint(fmt % ('S_elec', S_dict['S_e'], S_dict['S_e'] * temperature))
1411 vprint(fmt % ('S_vib', S_dict['S_v'], S_dict['S_v'] * temperature))
1412 vprint(
1413 fmt %
1414 ('S (1 bar -> P)',
1415 S_dict['S_p'],
1416 S_dict['S_p'] * temperature))
1417 vprint('-' * 49)
1418 vprint(fmt % ('S', S, S * temperature))
1419 vprint('=' * 49)
1420 return S
1422 def get_gibbs_energy(self, temperature: float,
1423 pressure: float,
1424 verbose: bool = True) -> float:
1425 """Returns the Gibbs free energy, in eV, in the ideal gas
1426 approximation at a specified temperature (K) and pressure (Pa)."""
1428 self.verbose = verbose
1429 vprint = self._vprint
1431 H = self.get_enthalpy(temperature, verbose=verbose)
1432 vprint('')
1433 S = self.get_entropy(temperature, pressure=pressure, verbose=verbose)
1434 G = H - temperature * S
1436 vprint('')
1437 vprint('Free energy components at T = %.2f K and P = %.1f Pa:' %
1438 (temperature, pressure))
1439 vprint('=' * 23)
1440 fmt = '%5s%15.3f eV'
1441 vprint(fmt % ('H', H))
1442 vprint(fmt % ('-T*S', -temperature * S))
1443 vprint('-' * 23)
1444 vprint(fmt % ('G', G))
1445 vprint('=' * 23)
1446 return G
1449class CrystalThermo(BaseThermoChem):
1450 """Class for calculating thermodynamic properties of a crystalline
1451 solid in the approximation that a lattice of N atoms behaves as a
1452 system of 3N independent harmonic oscillators.
1454 Inputs:
1456 phonon_DOS : list
1457 a list of the phonon density of states,
1458 where each value represents the phonon DOS at the vibrational energy
1459 value of the corresponding index in phonon_energies.
1461 phonon_energies : list
1462 a list of the range of vibrational energies (hbar*omega) over which
1463 the phonon density of states has been evaluated. This list should be
1464 the same length as phonon_DOS and integrating phonon_DOS over
1465 phonon_energies should yield approximately 3N, where N is the number
1466 of atoms per unit cell. If the first element of this list is
1467 zero-valued it will be deleted along with the first element of
1468 phonon_DOS. Units of vibrational energies are eV.
1470 potentialenergy : float
1471 the potential energy in eV (e.g., from atoms.get_potential_energy)
1472 (if potentialenergy is unspecified, then the methods of this
1473 class can be interpreted as the energy corrections)
1475 formula_units : int
1476 the number of formula units per unit cell. If unspecified, the
1477 thermodynamic quantities calculated will be listed on a
1478 per-unit-cell basis.
1479 """
1481 @classmethod
1482 def from_transition_state(cls, *args, **kwargs) -> "CrystalThermo":
1483 """
1484 Not yet implemented for CrystalThermo but needed to overwrite
1485 BaseThermoChem method.
1486 """
1487 raise NotImplementedError(
1488 "from_transition_state is not implemented for CrystalThermo."
1489 )
1491 def __init__(self, phonon_DOS, phonon_energies,
1492 formula_units=None, potentialenergy=0.):
1493 self.phonon_energies = phonon_energies
1494 self.phonon_DOS = phonon_DOS
1496 if formula_units:
1497 self.formula_units = formula_units
1498 self.potentialenergy = potentialenergy / formula_units
1499 else:
1500 self.formula_units = 0
1501 self.potentialenergy = potentialenergy
1503 def get_internal_energy(self, temperature, verbose=True):
1504 """Returns the internal energy, in eV, of crystalline solid
1505 at a specified temperature (K)."""
1507 self.verbose = verbose
1508 vprint = self._vprint
1509 fmt = '%-15s%13.4f eV'
1510 if self.formula_units == 0:
1511 vprint('Internal energy components at '
1512 'T = %.2f K,\non a per-unit-cell basis:' % temperature)
1513 else:
1514 vprint('Internal energy components at '
1515 'T = %.2f K,\non a per-formula-unit basis:' % temperature)
1516 vprint('=' * 31)
1518 U = 0.
1520 omega_e = self.phonon_energies
1521 dos_e = self.phonon_DOS
1522 if omega_e[0] == 0.:
1523 omega_e = np.delete(omega_e, 0)
1524 dos_e = np.delete(dos_e, 0)
1526 vprint(fmt % ('E_pot', self.potentialenergy))
1527 U += self.potentialenergy
1529 zpe_list = omega_e / 2.
1530 if self.formula_units == 0:
1531 zpe = trapezoid(zpe_list * dos_e, omega_e)
1532 else:
1533 zpe = trapezoid(zpe_list * dos_e, omega_e) / self.formula_units
1534 vprint(fmt % ('E_ZPE', zpe))
1535 U += zpe
1537 B = 1. / (units.kB * temperature)
1538 E_vib = omega_e / (np.exp(omega_e * B) - 1.)
1539 if self.formula_units == 0:
1540 E_phonon = trapezoid(E_vib * dos_e, omega_e)
1541 else:
1542 E_phonon = trapezoid(E_vib * dos_e, omega_e) / self.formula_units
1543 vprint(fmt % ('E_phonon', E_phonon))
1544 U += E_phonon
1546 vprint('-' * 31)
1547 vprint(fmt % ('U', U))
1548 vprint('=' * 31)
1549 return U
1551 def get_entropy(self, temperature, verbose=True):
1552 """Returns the entropy, in eV/K, of crystalline solid
1553 at a specified temperature (K)."""
1555 self.verbose = verbose
1556 vprint = self._vprint
1557 fmt = '%-15s%13.7f eV/K%13.4f eV'
1558 if self.formula_units == 0:
1559 vprint('Entropy components at '
1560 'T = %.2f K,\non a per-unit-cell basis:' % temperature)
1561 else:
1562 vprint('Entropy components at '
1563 'T = %.2f K,\non a per-formula-unit basis:' % temperature)
1564 vprint('=' * 49)
1565 vprint('%15s%13s %13s' % ('', 'S', 'T*S'))
1567 omega_e = self.phonon_energies
1568 dos_e = self.phonon_DOS
1569 if omega_e[0] == 0.:
1570 omega_e = np.delete(omega_e, 0)
1571 dos_e = np.delete(dos_e, 0)
1573 B = 1. / (units.kB * temperature)
1574 S_vib = (omega_e / (temperature * (np.exp(omega_e * B) - 1.)) -
1575 units.kB * np.log(1. - np.exp(-omega_e * B)))
1576 if self.formula_units == 0:
1577 S = trapezoid(S_vib * dos_e, omega_e)
1578 else:
1579 S = trapezoid(S_vib * dos_e, omega_e) / self.formula_units
1581 vprint('-' * 49)
1582 vprint(fmt % ('S', S, S * temperature))
1583 vprint('=' * 49)
1584 return S
1586 def get_helmholtz_energy(self, temperature, verbose=True):
1587 """Returns the Helmholtz free energy, in eV, of crystalline solid
1588 at a specified temperature (K)."""
1590 self.verbose = True
1591 vprint = self._vprint
1593 U = self.get_internal_energy(temperature, verbose=verbose)
1594 vprint('')
1595 S = self.get_entropy(temperature, verbose=verbose)
1596 F = U - temperature * S
1598 vprint('')
1599 if self.formula_units == 0:
1600 vprint('Helmholtz free energy components at '
1601 'T = %.2f K,\non a per-unit-cell basis:' % temperature)
1602 else:
1603 vprint('Helmholtz free energy components at '
1604 'T = %.2f K,\non a per-formula-unit basis:' % temperature)
1605 vprint('=' * 23)
1606 fmt = '%5s%15.4f eV'
1607 vprint(fmt % ('U', U))
1608 vprint(fmt % ('-T*S', -temperature * S))
1609 vprint('-' * 23)
1610 vprint(fmt % ('F', F))
1611 vprint('=' * 23)
1612 return F
1615def _clean_vib_energies(vib_energies: Sequence[complex],
1616 ignore_imag_modes: bool = False
1617 ) -> tuple[Sequence[float], int]:
1618 """Checks and deal with the presence of imaginary vibrational modes
1620 Also removes +0.j from real vibrational energies.
1622 Inputs:
1624 vib_energies : list
1625 a list of the vibrational energies
1627 ignore_imag_modes : bool
1628 If 'True', any imaginary frequencies will be removed.
1629 If 'False' (default), an error will be raised if imaginary
1630 frequencies are present.
1632 Outputs:
1634 vib_energies : list
1635 a list of the real vibrational energies.
1637 n_imag : int
1638 the number of imaginary frequencies treated.
1639 """
1640 # raise to a value, accept int and float in contrast to documentation
1641 if ignore_imag_modes:
1642 n_vib_energies = len(vib_energies)
1643 vib_energies = [np.real(v) for v in vib_energies if np.real(v) > 0]
1644 n_imag = n_vib_energies - len(vib_energies)
1645 if n_imag > 0:
1646 warn(f"{n_imag} imag modes removed", UserWarning)
1647 else:
1648 if sum(np.iscomplex(vib_energies)):
1649 raise ValueError('Imaginary vibrational energies are present.')
1650 n_imag = 0
1651 ret = [np.real(v) for v in vib_energies] # clear +0.j
1653 return ret, n_imag