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