# fmt: off
"""Modules for calculating thermochemical information from computational
outputs."""
import os
import sys
from abc import ABC, abstractmethod
from typing import Dict, Literal, Optional, Sequence, Tuple, Union
from warnings import catch_warnings, simplefilter, warn
import numpy as np
from scipy.integrate import trapezoid
from ase import Atoms, units
_GEOMETRY_OPTIONS = Literal['linear', 'nonlinear', 'monatomic']
_VIB_SELECTION_OPTIONS = Literal['exact', 'all', 'highest', 'abs_highest']
_FLOAT_OR_FLOATWITHDICT = Union[float, Tuple[float, Dict[str, float]]]
_FLOATWITHDICT = Tuple[float, Dict[str, float]]
def _sum_contributions(contrib_dicts: Dict[str, float]) -> float:
"""Combine a Dict of floats to their sum.
Ommits keys starting with an underscore.
"""
return np.sum([value for key, value in contrib_dicts.items()
if not key.startswith('_')])
[docs]
class AbstractMode(ABC):
"""Abstract base class for mode objects.
Input:
energy: float
Initialize the mode with its energy in eV.
Note: This energy refers to the vibrational energy of the mode.
Get it e.g. from ase.vibrations.Vibrations.get_energies.
"""
def __init__(self, energy: float) -> None:
self.energy = energy
@abstractmethod
def get_internal_energy(
self,
temperature: float,
contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
raise NotImplementedError
@abstractmethod
def get_entropy(self, temperature: float,
contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
raise NotImplementedError
[docs]
class HarmonicMode(AbstractMode):
"""Class for a single harmonic mode."""
[docs]
def get_ZPE_correction(self) -> float:
"""Returns the zero-point vibrational energy correction in eV."""
return 0.5 * self.energy
[docs]
def get_vib_energy_contribution(self, temperature: float) -> float:
"""Calculates the change in internal energy due to vibrations from
0K to the specified temperature for a set of vibrations given in
eV and a temperature given in Kelvin.
Returns the energy change in eV."""
kT = units.kB * temperature
return self.energy / (np.exp(self.energy / kT) - 1.)
[docs]
def get_vib_entropy_contribution(self, temperature: float) -> float:
"""Calculates the entropy due to vibrations given in eV and a
temperature given in Kelvin. Returns the entropy in eV/K."""
e = self.energy / units.kB / temperature
S_v = e / (np.exp(e) - 1.) - np.log(1. - np.exp(-e))
S_v *= units.kB
return S_v
[docs]
def get_internal_energy(self,
temperature: float,
contributions: bool
) -> _FLOAT_OR_FLOATWITHDICT:
"""Returns the internal energy in the harmonic approximation at a
specified temperature (K)."""
ret = {}
ret['ZPE'] = self.get_ZPE_correction()
ret['dU_v'] = self.get_vib_energy_contribution(temperature)
ret_sum = _sum_contributions(ret)
if contributions:
return ret_sum, ret
else:
return ret_sum
[docs]
def get_entropy(self,
temperature: float,
contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
"""Returns the entropy in the harmonic approximation at a specified
temperature (K).
Parameters
----------
temperature : float
The temperature in Kelvin.
contributions : bool
If True, return the contributions to the entropy as a dict too.
Here it will only contain the vibrational entropy.
Returns
-------
float or Tuple[float, Dict[str, float]]
"""
ret = {}
ret['S_v'] = self.get_vib_entropy_contribution(temperature)
if contributions:
return ret['S_v'], ret
else:
return ret['S_v']
[docs]
class RRHOMode(HarmonicMode):
"""Class for a single RRHO mode, including Grimme's scaling method
based on :doi:`10.1002/chem.201200497` and :doi:`10.1039/D1SC00621E`.
Inputs:
mean_inertia : float
The mean moment of inertia in amu*angstrom^2. Use
`np.mean(ase.Atoms.get_moments_of_inertia())` to calculate.
tau : float
the vibrational energy threshold in :math:`cm^{-1}`, named
:math:`\\tau` in :doi:`10.1039/D1SC00621E`.
Values close or equal to 0 will result in the standard harmonic
approximation. Defaults to :math:`35cm^{-1}`.
treat_int_energy : bool
Extend the msRRHO treatement to the internal thermal energy.
If False, only the entropy contribution as in Grimmmes paper is
modified according to the RRHO scheme. If True, the approach of
Otlyotov and Minenkov :doi:`10.1002/jcc.27129` is used to modify
the internal energy contribution.
"""
def __init__(self, energy: float,
mean_inertia: float,
tau: float = 35.0,
treat_int_energy: bool = False) -> None:
if np.iscomplex(energy):
raise ValueError(
"Imaginary frequencies are not allowed in RRHO mode.")
super().__init__(energy)
self._mean_inertia = mean_inertia
self._tau = tau
self._alpha = 4 # from paper 10.1002/chem.201200497
self.treat_int_energy = treat_int_energy
@property
def frequency(self) -> float:
return self.energy / units.invcm
def _head_gordon_damp(self, freq: float) -> float:
"""Head-Gordon damping function.
Equation 8 from :doi:`10.1002/chem.201200497`
Parameters
----------
freq : float
The frequency in the same unit as tau.
Returns
-------
float
"""
return 1 / (1 + (self._tau / freq)**self._alpha)
def _apply_head_gordon_damp(self, freq: float,
large_part: float,
small_part: float) -> float:
"""Apply the head-gordon damping scheme to two contributions.
Equation 7 from :doi:`10.1002/chem.201200497`
Returns the damped sum of the two contributions."""
part_one = self._head_gordon_damp(freq) * large_part
part_two = (1 - self._head_gordon_damp(freq)) * small_part
return part_one + part_two
[docs]
def get_RRHO_entropy_r(self, temperature: float) -> float:
"""Calculates the rotation of a rigid rotor for low frequency modes.
Equation numbering from :doi:`10.1002/chem.201200497`
Returns the entropy contribution in eV/K."""
kT = units._k * temperature
R = units._k * units._Nav # J / K / mol
B_av = (self._mean_inertia /
(units.kg * units.m**2)) # from amu/A^2 to kg m^2
# note, some codes use B_av = 1e-44 as in 10.1002/chem.201200497
# eq 4
omega = units._c * self.frequency * 1e2 # s^-1
mu = units._hplanck / (8 * np.pi**2 * omega) # kg m^2
# eq 5
mu_prime = (mu * B_av) / (mu + B_av) # kg m^2
# eq 6
x = np.sqrt(8 * np.pi**3 * mu_prime * kT / (units._hplanck)**2)
# filter zeros out and set them to zero
if x == 0.0:
log_x = 0.0
else:
log_x = np.log(x)
S_r = R * (0.5 + log_x) # J/(Js)^2
S_r *= units.J / units._Nav # J/K/mol to eV/K
return S_r
[docs]
def get_entropy(self,
temperature: float,
contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
ret = {}
ret['_S_vib_v'] = self.get_vib_entropy_contribution(temperature)
ret['_S_vib_r'] = self.get_RRHO_entropy_r(temperature)
ret['S_vib_damped'] = self._apply_head_gordon_damp(
self.frequency, ret['_S_vib_v'], ret['_S_vib_r'])
if contributions:
return ret['S_vib_damped'], ret
else:
return ret['S_vib_damped']
[docs]
def get_rrho_internal_energy_v_contribution(
self, temperature: float) -> float:
"""RRHO Vibrational Internal Energy Contribution from
:doi:`10.1002/jcc.27129`.
Returns the internal energy contribution in eV."""
# equation numbering from :doi:`10.1002/jcc.27129`
# eq 4
# hv = self.energy
theta = self.energy / units.kB
E_v = 0.5 + 1 / (np.exp(theta / temperature) - 1)
E_v *= self.energy # = theta * units.kB
return E_v
[docs]
@staticmethod
def get_rrho_internal_energy_r_contribution(temperature: float) -> float:
"""Calculates the rotation of a rigid rotor contribution.
Equation numbering from :doi:`10.1002/jcc.27129`
Returns the internal energy contribution in eV."""
# eq 5
R = units._k * units._Nav
E_r = R * temperature / 2
E_r *= units.J / units._Nav
return E_r
[docs]
def get_internal_energy(self,
temperature: float,
contributions: bool) -> _FLOAT_OR_FLOATWITHDICT:
"""Returns the internal energy in the msRRHO approximation at a
specified temperature (K).
If self.treat_int_energy is True, the approach of Otlyotov
and Minenkov :doi:`10.1002/jcc.27129` is used. Otherwise, the approach
of Grimme :doi:`10.1002/chem.201200497` is used.
"""
if self.treat_int_energy:
# Otlyotov and Minenkov approach with damping between vibrational
# and rotational contributions to the internal energy
# Note: The ZPE is not needed here, as the formula in the paper
# uses the "bottom of the well" as reference. See
# https://gaussian.com/wp-content/uploads/dl/thermo.pdf
# for more formulas
ret = {}
ret['_dU_vib_v'] = self.get_rrho_internal_energy_v_contribution(
temperature)
ret['_dU_vib_r'] = self.get_rrho_internal_energy_r_contribution(
temperature)
ret['dU_vib_damped'] = self._apply_head_gordon_damp(
self.frequency, ret['_dU_vib_v'], ret['_dU_vib_r'])
if contributions:
return ret['dU_vib_damped'], ret
else:
return ret['dU_vib_damped']
else:
# Grimme uses the Harmonic approach for the internal energy
return super().get_internal_energy(temperature, contributions)
[docs]
class BaseThermoChem(ABC):
"""Abstract base class containing common methods used in thermochemistry
calculations."""
def __init__(self,
vib_energies: Optional[Sequence[float]] = None,
atoms: Optional[Atoms] = None,
modes: Optional[Sequence[AbstractMode]] = None,
spin: Optional[float] = None) -> None:
if vib_energies is None and modes is None:
raise ValueError("Either vib_energies or modes must be provided.")
elif modes:
self.modes = modes
# in the monoatomic case, the vib_energies can be an empty list
else:
self._vib_energies = vib_energies
self.referencepressure = 1.0e5 # Pa
if atoms:
self.atoms = atoms
self.spin = spin
[docs]
@classmethod
def from_transition_state(cls, vib_energies: Sequence[complex],
*args, **kwargs) -> "BaseThermoChem":
"""Create a new instance for a transition state.
This will work just as the standard constructor, but will remove
one imaginary frequency from the given vib_energies first.
If there is more than one imaginary frequency, an error will be raised.
Returns
-------
BaseThermoChem instance
"""
if sum(np.iscomplex(vib_energies)):
# supress user warning
with catch_warnings():
simplefilter("ignore")
vib_e, n_imag = _clean_vib_energies(vib_energies, True)
if n_imag != 1:
raise ValueError("Not exactly one imaginary frequency found.")
else:
raise ValueError("No imaginary frequencies found in vib_energies.")
thermo = cls(vib_e, *args, **kwargs)
return thermo
[docs]
@staticmethod
def combine_contributions(
contrib_dicts: Sequence[Dict[str, float]]) -> Dict[str, float]:
"""Combine the contributions from multiple modes."""
ret: dict[str, float] = {}
for contrib_dict in contrib_dicts:
for key, value in contrib_dict.items():
if key in ret:
ret[key] += value
else:
ret[key] = value
return ret
[docs]
def print_contributions(
self, contributions: Dict[str, float], verbose: bool) -> None:
"""Print the contributions."""
if verbose:
fmt = "{:<15s}{:13.3f} eV"
for key, value in contributions.items():
# subvalues start with _
if key.startswith('_'):
key.replace('_', '... ')
self._vprint(fmt.format(key, value))
@abstractmethod
def get_internal_energy(self, temperature: float, verbose: bool) -> float:
raise NotImplementedError
@abstractmethod
def get_entropy(self, temperature: float,
pressure: float = units.bar,
verbose: bool = True) -> float:
raise NotImplementedError
@property
def vib_energies(self) -> np.ndarray:
"""Get the vibrational energies in eV.
For backwards compatibility, this will return the modes' energies if
they are available. Otherwise, it will return the stored vib_energies.
"""
# build from modes if available
if hasattr(self, 'modes'):
return np.array([mode.energy for mode in self.modes])
elif hasattr(self, '_vib_energies'):
return np.array(self._vib_energies)
else:
raise AttributeError("No vibrational energies available.")
@vib_energies.setter
def vib_energies(self, value: Sequence[float]) -> None:
"""For backwards compatibility, raise a deprecation warning."""
warn(
"The vib_energies attribute is deprecated and will be removed in a"
"future release. Please use the modes attribute instead."
"Setting this outside the constructor will not update the modes",
DeprecationWarning)
self._vib_energies = value
[docs]
def get_ZPE_correction(self) -> float:
"""Returns the zero-point vibrational energy correction in eV."""
return 0.5 * np.sum(self.vib_energies)
[docs]
@staticmethod
def get_ideal_translational_energy(temperature: float) -> float:
"""Returns the translational heat capacity times T in eV.
Parameters
----------
temperature : float
The temperature in Kelvin.
Returns
-------
float
"""
return 3. / 2. * units.kB * \
temperature # translational heat capacity (3-d gas)
[docs]
@staticmethod
def get_ideal_rotational_energy(geometry: _GEOMETRY_OPTIONS,
temperature: float) -> float:
"""Returns the rotational heat capacity times T in eV.
Parameters
----------
geometry : str
The geometry of the molecule. Options are 'nonlinear',
'linear', and 'monatomic'.
temperature : float
The temperature in Kelvin.
Returns
-------
float
"""
if geometry == 'nonlinear': # rotational heat capacity
Cv_r = 3. / 2. * units.kB
elif geometry == 'linear':
Cv_r = units.kB
elif geometry == 'monatomic':
Cv_r = 0.
else:
raise ValueError('Invalid geometry: %s' % geometry)
return Cv_r * temperature
[docs]
def get_ideal_trans_entropy(
self,
atoms: Atoms,
temperature: float) -> float:
"""Returns the translational entropy in eV/K.
Parameters
----------
atoms : ase.Atoms
The atoms object.
temperature : float
The temperature in Kelvin.
Returns
-------
float
"""
# Translational entropy (term inside the log is in SI units).
mass = sum(atoms.get_masses()) * units._amu # kg/molecule
S_t = (2 * np.pi * mass * units._k *
temperature / units._hplanck**2)**(3.0 / 2)
S_t *= units._k * temperature / self.referencepressure
S_t = units.kB * (np.log(S_t) + 5.0 / 2.0)
return S_t
[docs]
def get_vib_energy_contribution(self, temperature: float) -> float:
"""Calculates the change in internal energy due to vibrations from
0K to the specified temperature for a set of vibrations given in
eV and a temperature given in Kelvin.
Returns the energy change in eV."""
# Leftover for HinderedThermo class, delete once that is converted
# to use modes
kT = units.kB * temperature
dU = 0.
for energy in self.vib_energies:
dU += energy / (np.exp(energy / kT) - 1.)
return dU
[docs]
def get_vib_entropy_contribution(self,
temperature: float,
return_list: bool = False
) -> Union[float, Sequence[float]]:
"""Calculates the entropy due to vibrations for a set of vibrations
given in eV and a temperature given in Kelvin. Returns the entropy
in eV/K."""
kT = units.kB * temperature
S_v = 0.
energies = np.array(self.vib_energies)
energies /= kT # eV/ eV/K*K
S_v = energies / (np.exp(energies) - 1.) - \
np.log(1. - np.exp(-energies))
S_v *= units.kB
if return_list:
return S_v
else:
return np.sum(S_v)
def _vprint(self, text):
"""Print output if verbose flag True."""
if self.verbose:
sys.stdout.write(text + os.linesep)
[docs]
def get_ideal_entropy(
self,
temperature: float,
translation: bool = False,
vibration: bool = False,
rotation: bool = False,
geometry: Optional[_GEOMETRY_OPTIONS] = None,
electronic: bool = False,
pressure: Optional[float] = None,
symmetrynumber: Optional[int] = None) -> _FLOATWITHDICT:
"""Returns the entropy, in eV/K and a dict of the contributions.
Parameters
----------
temperature : float
The temperature in Kelvin.
translation : bool
Include translational entropy.
vibration : bool
Include vibrational entropy.
rotation : bool
Include rotational entropy.
geometry : str
The geometry of the molecule. Options are 'nonlinear',
'linear', and 'monatomic'.
electronic : bool
Include electronic entropy.
pressure : float
The pressure in Pa. Only needed for the translational entropy.
symmetrynumber : int
The symmetry number of the molecule. Only needed for linear and
nonlinear molecules.
Returns
-------
Tuple of one float and one dict
The float is the total entropy in eV/K.
The dict contains the contributions to the entropy.
"""
if (geometry in ['linear', 'nonlinear']) and (symmetrynumber is None):
raise ValueError(
'Symmetry number required for linear and nonlinear molecules.')
if not hasattr(self, 'atoms'):
raise ValueError(
'Atoms object required for ideal entropy calculation.')
if electronic and (self.spin is None):
raise ValueError(
'Spin value required for electronic entropy calculation.')
S: float = 0.0
ret = {}
if translation:
S_t = self.get_ideal_trans_entropy(self.atoms, temperature)
ret['S_t'] = S_t
S += S_t
if pressure:
# Pressure correction to translational entropy.
S_p = - units.kB * np.log(pressure / self.referencepressure)
S += S_p
ret['S_p'] = S_p
if rotation:
# Rotational entropy (term inside the log is in SI units).
if geometry == 'monatomic':
S_r = 0.0
elif geometry == 'nonlinear':
inertias = (self.atoms.get_moments_of_inertia() * units._amu /
1e10**2) # kg m^2
S_r = np.sqrt(np.pi * np.prod(inertias)) / symmetrynumber
S_r *= (8.0 * np.pi**2 * units._k * temperature /
units._hplanck**2)**(3.0 / 2.0)
S_r = units.kB * (np.log(S_r) + 3.0 / 2.0)
elif geometry == 'linear':
inertias = (self.atoms.get_moments_of_inertia() * units._amu /
(10.0**10)**2) # kg m^2
inertia = max(inertias) # should be two identical and one zero
S_r = (8 * np.pi**2 * inertia * units._k * temperature /
symmetrynumber / units._hplanck**2)
S_r = units.kB * (np.log(S_r) + 1.)
else:
raise RuntimeError(f"Invalid geometry: {geometry}")
S += S_r
ret['S_r'] = S_r
# Electronic entropy.
if electronic:
assert self.spin is not None # for mypy, error is raised above
S_e = units.kB * np.log(2 * self.spin + 1)
S += S_e
ret['S_e'] = S_e
# Vibrational entropy
if vibration:
S_v = self.get_vib_entropy_contribution(temperature,
return_list=False)
assert isinstance(S_v, float) # make mypy happy
S += S_v
ret['S_v'] = S_v
return S, ret
[docs]
class HarmonicThermo(BaseThermoChem):
"""Class for calculating thermodynamic properties in the approximation
that all degrees of freedom are treated harmonically. Often used for
adsorbates.
Note: This class does not include the translational and rotational
contributions to the entropy by default. Use the get_ideal_entropy method
for that and add them manually.
Inputs:
vib_energies : list
a list of the harmonic energies of the adsorbate (e.g., from
ase.vibrations.Vibrations.get_energies). The number of
energies should match the number of degrees of freedom of the
adsorbate; i.e., 3*n, where n is the number of atoms. Note that
this class does not check that the user has supplied the correct
number of energies. Units of energies are eV. Note, that if 'modes'
are given, these energies are not used.
potentialenergy : float
the potential energy in eV (e.g., from atoms.get_potential_energy)
(if potentialenergy is unspecified, then the methods of this
class can be interpreted as the energy corrections)
ignore_imag_modes : bool
If 'True', any imaginary frequencies will be removed in the
calculation of the thermochemical properties.
If 'False' (default), an error will be raised if any imaginary
frequencies are present.
modes : list of AbstractMode
A list of mode objects. If not provided, :class:`HarmonicMode` objects
will be created from the vib_energies. This is useful if you want to
replace individual modes with non-harmonic modes.
"""
def __init__(self, vib_energies: Sequence[complex],
potentialenergy: float = 0.,
ignore_imag_modes: bool = False,
modes: Optional[Sequence[AbstractMode]] = None) -> None:
# Check for imaginary frequencies.
vib_energies, n_imag = _clean_vib_energies(
vib_energies, ignore_imag_modes)
if modes is None:
modes = [HarmonicMode(energy) for energy in vib_energies]
super().__init__(modes=modes)
self.n_imag = n_imag
self.potentialenergy = potentialenergy
[docs]
def get_internal_energy(self, temperature: float,
verbose: bool = True) -> float:
"""Returns the internal energy, in eV, in the harmonic approximation
at a specified temperature (K)."""
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.3f eV'
vprint('Internal energy components at T = %.2f K:' % temperature)
vprint('=' * 31)
vprint(fmt % ('E_pot', self.potentialenergy))
U, contribs = zip(*[mode.get_internal_energy(
temperature, contributions=True) for mode in self.modes])
U = np.sum(U)
U += self.potentialenergy
self.print_contributions(self.combine_contributions(contribs), verbose)
vprint('-' * 31)
vprint(fmt % ('U', U))
vprint('=' * 31)
return U
[docs]
def get_entropy(self, temperature: float,
pressure: float = units.bar,
verbose: bool = True) -> float:
"""Returns the entropy, in eV/K at a specified temperature (K).
Note: This does not include the translational and rotational
contributions to the entropy. Use the get_ideal_entropy method
for that.
Parameters
----------
temperature : float
The temperature in Kelvin.
pressure : float
Not used, but kept for compatibility with other classes.
verbose : bool
If True, print the contributions to the entropy.
Returns
-------
float
"""
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.7f eV/K%13.3f eV'
vprint('Entropy components at T = %.2f K:' % temperature)
vprint('=' * 49)
vprint('%15s%13s %13s' % ('', 'S', 'T*S'))
S, contribs = zip(*[mode.get_entropy(
temperature, contributions=True) for mode in self.modes])
S = np.sum(S)
self.print_contributions(self.combine_contributions(contribs), verbose)
vprint('-' * 49)
vprint(fmt % ('S', S, S * temperature))
vprint('=' * 49)
return S
[docs]
def get_helmholtz_energy(self, temperature: float,
verbose: bool = True) -> float:
"""Returns the Helmholtz free energy, in eV, in the harmonic
approximation at a specified temperature (K)."""
self.verbose = True
vprint = self._vprint
U = self.get_internal_energy(temperature, verbose=verbose)
vprint('')
S = self.get_entropy(temperature, verbose=verbose)
F = U - temperature * S
vprint('')
vprint('Free energy components at T = %.2f K:' % temperature)
vprint('=' * 23)
fmt = '%5s%15.3f eV'
vprint(fmt % ('U', U))
vprint(fmt % ('-T*S', -temperature * S))
vprint('-' * 23)
vprint(fmt % ('F', F))
vprint('=' * 23)
return F
[docs]
class QuasiHarmonicThermo(HarmonicThermo):
"""Subclass of :class:`HarmonicThermo`, including the quasi-harmonic
approximation of Cramer, Truhlar and coworkers :doi:`10.1021/jp205508z`.
Note: This class not include the translational and rotational
contributions to the entropy by default. Use the get_ideal_entropy method
for that and add them manually.
Inputs:
vib_energies : list
a list of the energies of the vibrations (e.g., from
ase.vibrations.Vibrations.get_energies). The number of
energies should match the number of degrees of freedom of the
adsorbate; i.e., 3*n, where n is the number of atoms. Note that
this class does not check that the user has supplied the correct
number of energies. Units of energies are eV.
potentialenergy : float
the potential energy in eV (e.g., from atoms.get_potential_energy)
(if potentialenergy is unspecified, then the methods of this
class can be interpreted as the energy corrections)
ignore_imag_modes : bool
If 'True', any imaginary frequencies will be removed in the
calculation of the thermochemical properties.
If 'False' (default), an error will be raised if any imaginary
frequencies are present.
modes : list of AbstractMode
A list of mode objects. If not provided, :class:`HarmonicMode` objects
will be created from the raised vib_energies. This is useful if you want
to replace individual modes with non-harmonic modes.
raise_to : float
The value to which all frequencies smaller than this value will be
raised. Unit is eV. Defaults to :math:`100cm^{-1} = 0.012398 eV`.
"""
@staticmethod
def _raise(input: Sequence[float], raise_to: float) -> Sequence[float]:
return [raise_to if x < raise_to else x for x in input]
def __init__(self, vib_energies: Sequence[complex],
potentialenergy: float = 0.,
ignore_imag_modes: bool = False,
modes: Optional[Sequence[AbstractMode]] = None,
raise_to: float = 100 * units.invcm) -> None:
# Check for imaginary frequencies.
vib_energies, _ = _clean_vib_energies(
vib_energies, ignore_imag_modes)
# raise the low frequencies to a certain value
vib_energies = self._raise(vib_energies, raise_to)
if modes is None:
modes = [HarmonicMode(energy) for energy in vib_energies]
super().__init__(vib_energies,
potentialenergy=potentialenergy,
ignore_imag_modes=ignore_imag_modes,
modes=modes)
[docs]
class MSRRHOThermo(QuasiHarmonicThermo):
"""Subclass of :class:`QuasiHarmonicThermo`,
including Grimme's scaling method based on
:doi:`10.1002/chem.201200497` and :doi:`10.1039/D1SC00621E`.
Note: This class not include the translational and rotational
contributions to the entropy by default. Use the get_ideal_entropy method
for that and add them manually.
We enforce treating imaginary modes as Grimme suggests (converting
them to real by multiplying them with :math:`-i`). So make sure to check
your input energies.
Inputs:
vib_energies : list
a list of the energies of the vibrations (e.g., from
ase.vibrations.Vibrations.get_energies). The number of
energies should match the number of degrees of freedom of the
atoms object; i.e., 3*n, where n is the number of atoms. Note that
this class does not check that the user has supplied the correct
number of energies. Units of energies are eV. Note, that if 'modes'
are given, these energies are not used.
atoms: an ASE atoms object
used to calculate rotational moments of inertia and molecular mass
tau : float
the vibrational energy threshold in :math:`cm^{-1}`, namcomplexed
:math:`\\tau` in :doi:`10.1039/D1SC00621E`.
Values close or equal to 0 will result in the standard harmonic
approximation. Defaults to :math:`35cm^{-1}`.
nu_scal : float
Linear scaling factor for the vibrational frequencies. Named
:math:`\\nu_{scal}` in :doi:`10.1039/D1SC00621E`.
Defaults to 1.0, check the `Truhlar group database
<https://comp.chem.umn.edu/freqscale/index.html>`_
for values corresponding to your level of theory.
Note that for `\\nu_{scal}=1.0` this method is equivalent to
the quasi-RRHO method in :doi:`10.1002/chem.201200497`.
treat_int_energy : bool
Extend the msRRHO treatement to the internal energy. If False, only
the entropy contribution as in Grimmmes paper is considered.
If true, the approach of Otlyotov and Minenkov
:doi:`10.1002/jcc.27129` is used.
modes : list of AbstractMode
A list of mode objects. If not provided, :class:`RRHOMode` objects will
be created from the raised vib_energies. This is useful if you want to
replace individual modes with non-harmonic modes.
"""
def __init__(self, vib_energies: Sequence[complex], atoms: Atoms,
potentialenergy: float = 0.,
tau: float = 35., nu_scal: float = 1.0,
treat_int_energy: bool = False,
modes: Optional[Sequence[AbstractMode]] = None) -> None:
inertia = np.mean(atoms.get_moments_of_inertia())
self.atoms = atoms
# clean the energies by converting imaginary to real
# i.e. multiply with -i
vib_e: Sequence[float] # make mypy happy
vib_e = [np.imag(v) if np.iscomplex(v)
else np.real(v) for v in vib_energies]
self.nu_scal = nu_scal
# scale the frequencies (i.e. energies) before passing them on
vib_e_scaled = [float(v) for v in np.multiply(vib_e, nu_scal)]
if modes is None:
modes = [RRHOMode(energy, inertia,
tau=tau,
treat_int_energy=treat_int_energy
) for energy in vib_e_scaled]
super().__init__(vib_e_scaled,
potentialenergy=potentialenergy,
ignore_imag_modes=False,
modes=modes,
raise_to=0.0)
self.treat_int_energy = treat_int_energy
[docs]
class HinderedThermo(BaseThermoChem):
"""Class for calculating thermodynamic properties in the hindered
translator and hindered rotor model where all but three degrees of
freedom are treated as harmonic vibrations, two are treated as
hindered translations, and one is treated as a hindered rotation.
Inputs:
vib_energies : list
a list of all the vibrational energies of the adsorbate (e.g., from
ase.vibrations.Vibrations.get_energies). If atoms is not provided,
then the number of energies must match the number of degrees of freedom
of the adsorbate; i.e., 3*n, where n is the number of atoms. Note
that this class does not check that the user has supplied the
correct number of energies.
Units of energies are eV.
trans_barrier_energy : float
the translational energy barrier in eV. This is the barrier for an
adsorbate to diffuse on the surface.
rot_barrier_energy : float
the rotational energy barrier in eV. This is the barrier for an
adsorbate to rotate about an axis perpendicular to the surface.
sitedensity : float
density of surface sites in cm^-2
rotationalminima : integer
the number of equivalent minima for an adsorbate's full rotation.
For example, 6 for an adsorbate on an fcc(111) top site
potentialenergy : float
the potential energy in eV (e.g., from atoms.get_potential_energy)
(if potentialenergy is unspecified, then the methods of this class
can be interpreted as the energy corrections)
mass : float
the mass of the adsorbate in amu (if mass is unspecified, then it will
be calculated from the atoms class)
inertia : float
the reduced moment of inertia of the adsorbate in amu*Ang^-2
(if inertia is unspecified, then it will be calculated from the
atoms class)
atoms : an ASE atoms object
used to calculate rotational moments of inertia and molecular mass
symmetrynumber : integer
symmetry number of the adsorbate. This is the number of symmetric arms
of the adsorbate and depends upon how it is bound to the surface.
For example, propane bound through its end carbon has a symmetry
number of 1 but propane bound through its middle carbon has a symmetry
number of 2. (if symmetrynumber is unspecified, then the default is 1)
ignore_imag_modes : bool
If 'True', any imaginary frequencies present after the 3N-3 cut will
be removed in the calculation of the thermochemical properties.
If 'False' (default), an error will be raised if imaginary frequencies
are present after the 3N-3 cut.
"""
def __init__(self, vib_energies, trans_barrier_energy, rot_barrier_energy,
sitedensity, rotationalminima, potentialenergy=0.,
mass=None, inertia=None, atoms=None, symmetrynumber=1,
ignore_imag_modes: bool = False) -> None:
self.trans_barrier_energy = trans_barrier_energy * units._e
self.rot_barrier_energy = rot_barrier_energy * units._e
self.area = 1. / sitedensity / 100.0**2
self.rotationalminima = rotationalminima
self.potentialenergy = potentialenergy
self.atoms = atoms
self.symmetry = symmetrynumber
# Sort the vibrations
vib_energies = list(vib_energies)
vib_energies.sort(key=np.abs)
# Keep only the relevant vibrational energies (3N-3)
if atoms:
vib_energies = vib_energies[-(3 * len(atoms) - 3):]
else:
vib_energies = vib_energies[-(len(vib_energies) - 3):]
# Check for imaginary frequencies.
vib_energies, n_imag = _clean_vib_energies(
vib_energies, ignore_imag_modes)
super().__init__(vib_energies=vib_energies)
self.n_imag = n_imag
if (mass or atoms) and (inertia or atoms):
if mass:
self.mass = mass * units._amu
elif atoms:
self.mass = np.sum(atoms.get_masses()) * units._amu
if inertia:
self.inertia = inertia * units._amu / units.m**2
elif atoms:
self.inertia = (atoms.get_moments_of_inertia()[2] *
units._amu / units.m**2)
else:
raise RuntimeError('Either mass and inertia of the '
'adsorbate must be specified or '
'atoms must be specified.')
# Calculate hindered translational and rotational frequencies
self.freq_t = np.sqrt(self.trans_barrier_energy /
(2 * self.mass * self.area))
self.freq_r = 1. / (2 * np.pi) * np.sqrt(self.rotationalminima**2 *
self.rot_barrier_energy /
(2 * self.inertia))
[docs]
def get_internal_energy(self, temperature, verbose=True):
"""Returns the internal energy (including the zero point energy),
in eV, in the hindered translator and hindered rotor model at a
specified temperature (K)."""
from scipy.special import iv
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.3f eV'
vprint('Internal energy components at T = %.2f K:' % temperature)
vprint('=' * 31)
U = 0.
vprint(fmt % ('E_pot', self.potentialenergy))
U += self.potentialenergy
# Translational Energy
T_t = units._k * temperature / (units._hplanck * self.freq_t)
R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)
dU_t = 2 * (-1. / 2 - 1. / T_t / (2 + 16 * R_t) + R_t / 2 / T_t -
R_t / 2 / T_t *
iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +
1. / T_t / (np.exp(1. / T_t) - 1))
dU_t *= units.kB * temperature
vprint(fmt % ('E_trans', dU_t))
U += dU_t
# Rotational Energy
T_r = units._k * temperature / (units._hplanck * self.freq_r)
R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)
dU_r = (-1. / 2 - 1. / T_r / (2 + 16 * R_r) + R_r / 2 / T_r -
R_r / 2 / T_r *
iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +
1. / T_r / (np.exp(1. / T_r) - 1))
dU_r *= units.kB * temperature
vprint(fmt % ('E_rot', dU_r))
U += dU_r
# Vibrational Energy
dU_v = self.get_vib_energy_contribution(temperature)
vprint(fmt % ('E_vib', dU_v))
U += dU_v
# Zero Point Energy
dU_zpe = self.get_zero_point_energy()
vprint(fmt % ('E_ZPE', dU_zpe))
U += dU_zpe
vprint('-' * 31)
vprint(fmt % ('U', U))
vprint('=' * 31)
return U
[docs]
def get_zero_point_energy(self, verbose=True):
"""Returns the zero point energy, in eV, in the hindered
translator and hindered rotor model"""
zpe_t = 2 * (1. / 2 * self.freq_t * units._hplanck / units._e)
zpe_r = 1. / 2 * self.freq_r * units._hplanck / units._e
zpe_v = self.get_ZPE_correction()
zpe = zpe_t + zpe_r + zpe_v
return zpe
[docs]
def get_entropy(self, temperature,
pressure=units.bar,
verbose=True):
"""Returns the entropy, in eV/K, in the hindered translator
and hindered rotor model at a specified temperature (K)."""
from scipy.special import iv
self.verbose = verbose
vrpint = self._vprint
fmt = '%-15s%13.7f eV/K%13.3f eV'
vrpint('Entropy components at T = %.2f K:' % temperature)
vrpint('=' * 49)
vrpint('%15s%13s %13s' % ('', 'S', 'T*S'))
S = 0.
# Translational Entropy
T_t = units._k * temperature / (units._hplanck * self.freq_t)
R_t = self.trans_barrier_energy / (units._hplanck * self.freq_t)
S_t = 2 * (-1. / 2 + 1. / 2 * np.log(np.pi * R_t / T_t) -
R_t / 2 / T_t *
iv(1, R_t / 2 / T_t) / iv(0, R_t / 2 / T_t) +
np.log(iv(0, R_t / 2 / T_t)) +
1. / T_t / (np.exp(1. / T_t) - 1) -
np.log(1 - np.exp(-1. / T_t)))
S_t *= units.kB
vrpint(fmt % ('S_trans', S_t, S_t * temperature))
S += S_t
# Rotational Entropy
T_r = units._k * temperature / (units._hplanck * self.freq_r)
R_r = self.rot_barrier_energy / (units._hplanck * self.freq_r)
S_r = (-1. / 2 + 1. / 2 * np.log(np.pi * R_r / T_r) -
np.log(self.symmetry) -
R_r / 2 / T_r * iv(1, R_r / 2 / T_r) / iv(0, R_r / 2 / T_r) +
np.log(iv(0, R_r / 2 / T_r)) +
1. / T_r / (np.exp(1. / T_r) - 1) -
np.log(1 - np.exp(-1. / T_r)))
S_r *= units.kB
vrpint(fmt % ('S_rot', S_r, S_r * temperature))
S += S_r
# Vibrational Entropy
S_v = self.get_vib_entropy_contribution(temperature)
vrpint(fmt % ('S_vib', S_v, S_v * temperature))
S += S_v
# Concentration Related Entropy
N_over_A = np.exp(1. / 3) * (10.0**5 /
(units._k * temperature))**(2. / 3)
S_c = 1 - np.log(N_over_A) - np.log(self.area)
S_c *= units.kB
vrpint(fmt % ('S_con', S_c, S_c * temperature))
S += S_c
vrpint('-' * 49)
vrpint(fmt % ('S', S, S * temperature))
vrpint('=' * 49)
return S
[docs]
def get_helmholtz_energy(self, temperature, verbose=True):
"""Returns the Helmholtz free energy, in eV, in the hindered
translator and hindered rotor model at a specified temperature
(K)."""
self.verbose = True
vprint = self._vprint
U = self.get_internal_energy(temperature, verbose=verbose)
vprint('')
S = self.get_entropy(temperature, verbose=verbose)
F = U - temperature * S
vprint('')
vprint('Free energy components at T = %.2f K:' % temperature)
vprint('=' * 23)
fmt = '%5s%15.3f eV'
vprint(fmt % ('U', U))
vprint(fmt % ('-T*S', -temperature * S))
vprint('-' * 23)
vprint(fmt % ('F', F))
vprint('=' * 23)
return F
[docs]
class IdealGasThermo(BaseThermoChem):
"""Class for calculating thermodynamic properties of a molecule
based on statistical mechanical treatments in the ideal gas
approximation.
Inputs for enthalpy calculations:
vib_energies : list
a list of the vibrational energies of the molecule (e.g., from
ase.vibrations.Vibrations.get_energies). The number of expected
vibrations is calculated by the geometry and the number of atoms.
By default, the lowest energies are cut to enforce the expected
length. By setting vib_selection, this behaviour can be changed.
If either atoms and natoms is unspecified, the full list is used.
Units are eV.
geometry : 'monatomic', 'linear', or 'nonlinear'
geometry of the molecule
potentialenergy : float
the potential energy in eV (e.g., from atoms.get_potential_energy)
(if potentialenergy is unspecified, then the methods of this
class can be interpreted as the energy corrections)
natoms : integer
the number of atoms, used along with 'geometry' to determine how
many vibrations to use. (Not needed if an atoms object is supplied
in 'atoms' or if the user desires the entire list of vibrations
to be used.)
vib_selection : 'exact', 'highest' (default), 'abs_highest', 'all'
selection of input vibrational energies considered to be true
vibrations (excluding translations and rotations) implied by the
geometry and number of atoms. Only applied if number of atoms
is provided, either with natoms or atoms.
- 'exact': no selection; the number of input vibrational energies must
be equal to the number of true vibrations
- 'highest': select vibrational energies whose square is the highest
, i.e. large real energies followed by small real energies, small
imaginary energies and large imaginary energies. Will not catch
unrealistically large imaginary frequencies, resulting from e.g.
unrelaxed molecules. This is the default.
- 'abs_highest': select vibrational energies whose absolute number are
highest, real or imaginary. Will fail intentionally for
unrealistically large imaginary frequencies, but also unintentionally
for small real energies, corresponding e.g. to frustrated rotations
in a larger molecule.
- 'all': Use all input vibrational energies without checking the
number of them.
ignore_imag_modes : bool
If 'True', any imaginary frequencies present after the 3N-3 cut will
be removed in the calculation of the thermochemical properties.
If 'False' (default), an error will be raised if imaginary frequencies
are present after the 3N-3 cut.
Extra inputs needed for entropy / free energy calculations:
atoms : an ASE atoms object
used to calculate rotational moments of inertia and molecular mass
symmetrynumber : integer
symmetry number of the molecule. See, for example, Table 10.1 and
Appendix B of C. Cramer "Essentials of Computational Chemistry",
2nd Ed.
spin : float
the total electronic spin. (0 for molecules in which all electrons
are paired, 0.5 for a free radical with a single unpaired electron,
1.0 for a triplet with two unpaired electrons, such as O_2.)
"""
def __init__(self, vib_energies: Sequence[complex],
geometry: _GEOMETRY_OPTIONS,
potentialenergy: float = 0.,
atoms: Optional[Atoms] = None,
symmetrynumber: Optional[int] = None,
spin: Optional[float] = None,
natoms: Optional[int] = None,
vib_selection: Optional[_VIB_SELECTION_OPTIONS] = 'highest',
ignore_imag_modes: bool = False,
modes: Optional[Sequence[AbstractMode]] = None) -> None:
self.potentialenergy = potentialenergy
self.geometry = geometry
self.sigma = symmetrynumber
if natoms is None and atoms:
natoms = len(atoms)
# Preliminary sorting vibrations by square
vib_energies = list(vib_energies)
vib_energies.sort(key=lambda f: (f ** 2).real)
if natoms and vib_selection != 'all':
# Determine number of true vibrations from geometry
if geometry == 'nonlinear':
num_vibs = 3 * natoms - 6
elif geometry == 'linear':
num_vibs = 3 * natoms - 5
elif geometry == 'monatomic':
num_vibs = 0
else:
raise ValueError(f"Unsupported geometry: {geometry}")
if num_vibs < 0:
raise ValueError(
f"Number of atoms ({natoms}) too small for " +
f"geometry '{geometry}'."
)
if vib_selection == 'exact':
# Demand the expected number of true vibrations
if len(vib_energies) != num_vibs:
raise ValueError(
f"{num_vibs} vibrational energies expected, " +
f"{len(vib_energies)} received.\n" +
"To select true vibrational energies automatically," +
" set vib_selection."
)
elif vib_selection in ('highest', 'abs_highest'):
# Cut vibrations to get the expected number of true vibrations
if num_vibs == 0:
vib_energies = []
else:
if vib_selection == 'abs_highest':
vib_energies.sort(key=np.abs)
vib_energies = vib_energies[-num_vibs:]
if len(vib_energies) != num_vibs:
raise ValueError(
f"Too few vibration modes ({len(vib_energies)}) " +
"after selection. {num_vibs} expected"
)
else:
raise ValueError(
f"(Unsupported vib_selection: {vib_selection})"
)
elif not natoms and not atoms and vib_selection != 'all':
raise ValueError(
"Either natoms or atoms must be specified to select " +
"true vibrational energies. Else, set vib_selection to 'all'."
)
# Check for imaginary frequencies.
vib_energies, n_imag = _clean_vib_energies(
vib_energies, ignore_imag_modes)
super().__init__(vib_energies=vib_energies,
atoms=atoms,
spin=spin)
self.n_imag = n_imag
[docs]
def get_internal_energy(self, temperature: float,
verbose: bool = True) -> float:
"""Returns the internal energy, in eV, in the ideal gas approximation
at a specified temperature (K)."""
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.3f eV'
vprint('Enthalpy components at T = %.2f K:' % temperature)
vprint('=' * 31)
U = 0.
vprint(fmt % ('E_pot', self.potentialenergy))
U += self.potentialenergy
zpe = self.get_ZPE_correction()
vprint(fmt % ('E_ZPE', zpe))
U += zpe
Cv_tT = self.get_ideal_translational_energy(temperature)
vprint(fmt % ('Cv_trans (0->T)', Cv_tT))
U += Cv_tT
Cv_rT = self.get_ideal_rotational_energy(self.geometry, temperature)
vprint(fmt % ('Cv_rot (0->T)', Cv_rT))
U += Cv_rT
dU_v = self.get_vib_energy_contribution(temperature)
vprint(fmt % ('Cv_vib (0->T)', dU_v))
U += dU_v
vprint('-' * 31)
vprint(fmt % ('U', U))
vprint('=' * 31)
return U
[docs]
def get_enthalpy(self, temperature: float,
verbose: bool = True) -> float:
"""Returns the enthalpy, in eV, in the ideal gas approximation
at a specified temperature (K)."""
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.3f eV'
vprint('Enthalpy components at T = %.2f K:' % temperature)
vprint('=' * 31)
H = 0.
H += self.get_internal_energy(temperature, verbose=verbose)
Cp_corr = units.kB * temperature
vprint(fmt % ('(C_v -> C_p)', Cp_corr))
H += Cp_corr
vprint('-' * 31)
vprint(fmt % ('H', H))
vprint('=' * 31)
return H
[docs]
def get_entropy(self, temperature: float,
pressure: float = units.bar,
verbose: bool = True) -> float:
"""Returns the entropy, in eV/K, in the ideal gas approximation
at a specified temperature (K) and pressure (Pa)."""
if self.atoms is None or self.sigma is None or self.spin is None:
raise RuntimeError('atoms, symmetrynumber, and spin must be '
'specified for entropy and free energy '
'calculations.')
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.7f eV/K%13.3f eV'
vprint(f'Entropy components at T = {temperature:.2f} K and'
f' P = {pressure:.1f} Pa:')
vprint('=' * 49)
vprint('{"":15s}{"S":13s} {"T*S:13s}')
S, S_dict = self.get_ideal_entropy(temperature,
translation=True,
vibration=True,
rotation=True,
geometry=self.geometry,
electronic=True,
pressure=pressure,
symmetrynumber=self.sigma)
vprint(
fmt %
('S_trans (1 bar)',
S_dict['S_t'],
S_dict['S_t'] *
temperature))
vprint(fmt % ('S_rot', S_dict['S_r'], S_dict['S_r'] * temperature))
vprint(fmt % ('S_elec', S_dict['S_e'], S_dict['S_e'] * temperature))
vprint(fmt % ('S_vib', S_dict['S_v'], S_dict['S_v'] * temperature))
vprint(
fmt %
('S (1 bar -> P)',
S_dict['S_p'],
S_dict['S_p'] * temperature))
vprint('-' * 49)
vprint(fmt % ('S', S, S * temperature))
vprint('=' * 49)
return S
[docs]
def get_gibbs_energy(self, temperature: float,
pressure: float,
verbose: bool = True) -> float:
"""Returns the Gibbs free energy, in eV, in the ideal gas
approximation at a specified temperature (K) and pressure (Pa)."""
self.verbose = verbose
vprint = self._vprint
H = self.get_enthalpy(temperature, verbose=verbose)
vprint('')
S = self.get_entropy(temperature, pressure=pressure, verbose=verbose)
G = H - temperature * S
vprint('')
vprint('Free energy components at T = %.2f K and P = %.1f Pa:' %
(temperature, pressure))
vprint('=' * 23)
fmt = '%5s%15.3f eV'
vprint(fmt % ('H', H))
vprint(fmt % ('-T*S', -temperature * S))
vprint('-' * 23)
vprint(fmt % ('G', G))
vprint('=' * 23)
return G
[docs]
class CrystalThermo(BaseThermoChem):
"""Class for calculating thermodynamic properties of a crystalline
solid in the approximation that a lattice of N atoms behaves as a
system of 3N independent harmonic oscillators.
Inputs:
phonon_DOS : list
a list of the phonon density of states,
where each value represents the phonon DOS at the vibrational energy
value of the corresponding index in phonon_energies.
phonon_energies : list
a list of the range of vibrational energies (hbar*omega) over which
the phonon density of states has been evaluated. This list should be
the same length as phonon_DOS and integrating phonon_DOS over
phonon_energies should yield approximately 3N, where N is the number
of atoms per unit cell. If the first element of this list is
zero-valued it will be deleted along with the first element of
phonon_DOS. Units of vibrational energies are eV.
potentialenergy : float
the potential energy in eV (e.g., from atoms.get_potential_energy)
(if potentialenergy is unspecified, then the methods of this
class can be interpreted as the energy corrections)
formula_units : int
the number of formula units per unit cell. If unspecified, the
thermodynamic quantities calculated will be listed on a
per-unit-cell basis.
"""
[docs]
@classmethod
def from_transition_state(cls, *args, **kwargs) -> "CrystalThermo":
"""
Not yet implemented for CrystalThermo but needed to overwrite
BaseThermoChem method.
"""
raise NotImplementedError(
"from_transition_state is not implemented for CrystalThermo."
)
def __init__(self, phonon_DOS, phonon_energies,
formula_units=None, potentialenergy=0.):
self.phonon_energies = phonon_energies
self.phonon_DOS = phonon_DOS
if formula_units:
self.formula_units = formula_units
self.potentialenergy = potentialenergy / formula_units
else:
self.formula_units = 0
self.potentialenergy = potentialenergy
[docs]
def get_internal_energy(self, temperature, verbose=True):
"""Returns the internal energy, in eV, of crystalline solid
at a specified temperature (K)."""
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.4f eV'
if self.formula_units == 0:
vprint('Internal energy components at '
'T = %.2f K,\non a per-unit-cell basis:' % temperature)
else:
vprint('Internal energy components at '
'T = %.2f K,\non a per-formula-unit basis:' % temperature)
vprint('=' * 31)
U = 0.
omega_e = self.phonon_energies
dos_e = self.phonon_DOS
if omega_e[0] == 0.:
omega_e = np.delete(omega_e, 0)
dos_e = np.delete(dos_e, 0)
vprint(fmt % ('E_pot', self.potentialenergy))
U += self.potentialenergy
zpe_list = omega_e / 2.
if self.formula_units == 0:
zpe = trapezoid(zpe_list * dos_e, omega_e)
else:
zpe = trapezoid(zpe_list * dos_e, omega_e) / self.formula_units
vprint(fmt % ('E_ZPE', zpe))
U += zpe
B = 1. / (units.kB * temperature)
E_vib = omega_e / (np.exp(omega_e * B) - 1.)
if self.formula_units == 0:
E_phonon = trapezoid(E_vib * dos_e, omega_e)
else:
E_phonon = trapezoid(E_vib * dos_e, omega_e) / self.formula_units
vprint(fmt % ('E_phonon', E_phonon))
U += E_phonon
vprint('-' * 31)
vprint(fmt % ('U', U))
vprint('=' * 31)
return U
[docs]
def get_entropy(self, temperature, verbose=True):
"""Returns the entropy, in eV/K, of crystalline solid
at a specified temperature (K)."""
self.verbose = verbose
vprint = self._vprint
fmt = '%-15s%13.7f eV/K%13.4f eV'
if self.formula_units == 0:
vprint('Entropy components at '
'T = %.2f K,\non a per-unit-cell basis:' % temperature)
else:
vprint('Entropy components at '
'T = %.2f K,\non a per-formula-unit basis:' % temperature)
vprint('=' * 49)
vprint('%15s%13s %13s' % ('', 'S', 'T*S'))
omega_e = self.phonon_energies
dos_e = self.phonon_DOS
if omega_e[0] == 0.:
omega_e = np.delete(omega_e, 0)
dos_e = np.delete(dos_e, 0)
B = 1. / (units.kB * temperature)
S_vib = (omega_e / (temperature * (np.exp(omega_e * B) - 1.)) -
units.kB * np.log(1. - np.exp(-omega_e * B)))
if self.formula_units == 0:
S = trapezoid(S_vib * dos_e, omega_e)
else:
S = trapezoid(S_vib * dos_e, omega_e) / self.formula_units
vprint('-' * 49)
vprint(fmt % ('S', S, S * temperature))
vprint('=' * 49)
return S
[docs]
def get_helmholtz_energy(self, temperature, verbose=True):
"""Returns the Helmholtz free energy, in eV, of crystalline solid
at a specified temperature (K)."""
self.verbose = True
vprint = self._vprint
U = self.get_internal_energy(temperature, verbose=verbose)
vprint('')
S = self.get_entropy(temperature, verbose=verbose)
F = U - temperature * S
vprint('')
if self.formula_units == 0:
vprint('Helmholtz free energy components at '
'T = %.2f K,\non a per-unit-cell basis:' % temperature)
else:
vprint('Helmholtz free energy components at '
'T = %.2f K,\non a per-formula-unit basis:' % temperature)
vprint('=' * 23)
fmt = '%5s%15.4f eV'
vprint(fmt % ('U', U))
vprint(fmt % ('-T*S', -temperature * S))
vprint('-' * 23)
vprint(fmt % ('F', F))
vprint('=' * 23)
return F
def _clean_vib_energies(vib_energies: Sequence[complex],
ignore_imag_modes: bool = False
) -> Tuple[Sequence[float], int]:
"""Checks and deal with the presence of imaginary vibrational modes
Also removes +0.j from real vibrational energies.
Inputs:
vib_energies : list
a list of the vibrational energies
ignore_imag_modes : bool
If 'True', any imaginary frequencies will be removed.
If 'False' (default), an error will be raised if imaginary
frequencies are present.
Outputs:
vib_energies : list
a list of the real vibrational energies.
n_imag : int
the number of imaginary frequencies treated.
"""
# raise to a value, accept int and float in contrast to documentation
if ignore_imag_modes:
n_vib_energies = len(vib_energies)
vib_energies = [np.real(v) for v in vib_energies if np.real(v) > 0]
n_imag = n_vib_energies - len(vib_energies)
if n_imag > 0:
warn(f"{n_imag} imag modes removed", UserWarning)
else:
if sum(np.iscomplex(vib_energies)):
raise ValueError('Imaginary vibrational energies are present.')
n_imag = 0
ret = [np.real(v) for v in vib_energies] # clear +0.j
return ret, n_imag
