Repository: duartegroup/autodE (python)

Issue:
A new feature for the calculation of thermochemical contributions
Hi, Tom I recommend a new feature for the calculation of thermochemical contributions.
qRRHO treatment proposed by Grimme only applies to entropy.  Recently,  another paper (J. Comput. Chem., 44, 1807 (2023))  proved that using Grimme's scheme to calculate the enthalpy contribution will ensure the internal consistency and improve the quality of Gibbs free energies (ΔGr).

This approach can be easily realized by incorporating a few lines of code. For example:
```
U_free_rotor = 0.5 * SIConstants.k_b * temp
U_qRRHO = w * U_RROH + (1 - w) * U_free_rotor
 ```
where w is calculation by(line 378 in autode.thermochemistry.igm.py):
`w = 1.0 / (1.0 + (omega_0 / freq) ** alpha)`

Pull request:
Minenkov thermochem
<!--
Please replace _#ISSUE_ with just e.g. #12 if this PR resolves issue 12.
-->
## Resolves #289 

Adds the suggested method and tidies up the `igm` module. There's much more that could be better but 🤷🏼, at least it's better than it was.

***

## Checklist

* [x] The changes include an associated explanation of how/why
* [x] Test pass
* [x] Documentation has been updated
* [x] Changelog has been updated


Patch:
diff --git a/autode/config.py b/autode/config.py
index 1fd883a..0fee60c 100644
--- a/autode/config.py
+++ b/autode/config.py
@@ -148,11 +148,13 @@ class _ConfigClass:
     #
     # Method to treat low frequency modes (LFMs). Either standard RRHO ('igm'),
     # Truhlar's method where all frequencies below a threshold are scaled to
-    # a shifted value (see J. Phys. Chem. B, 2011, 115, 14556) or Grimme's
+    # a shifted value (see J. Phys. Chem. B, 2011, 115, 14556), Grimme's
     # method of interpolating between HO and RR (i.e. qRRHO, see
-    # Chem. Eur. J. 2012, 18, 9955)
+    # Chem. Eur. J. 2012, 18, 9955), or 'minenkov' where free rotor/vibrational
+    # interpolation is useed for U and S (i.e. mRRHO, see
+    # J. Comput. Chem., 2023 44, 1807)
     #
-    # One of: 'igm', 'truhlar', 'grimme'
+    # One of: 'igm', 'truhlar', 'grimme', 'minenkov'
     lfm_method = "grimme"
     #
     # Parameters for Grimme's method (only used when lfm_method='grimme'),
diff --git a/autode/thermochemistry/igm.py b/autode/thermochemistry/igm.py
index 3d3d273..783dda4 100644
--- a/autode/thermochemistry/igm.py
+++ b/autode/thermochemistry/igm.py
@@ -9,10 +9,15 @@ All calculations performed in SI units, for simplicity.
 """
 import numpy as np
 from enum import Enum
+from typing import TYPE_CHECKING, Any
 from autode.log import logger
 from autode.config import Config
 from autode.constants import Constants
-from autode.values import FreeEnergyCont, EnthalpyCont
+from autode.values import FreeEnergyCont, EnthalpyCont, Frequency
+
+
+if TYPE_CHECKING:
+    from autode.species.species import Species
 
 
 class SIConstants:
@@ -34,9 +39,27 @@ class LFMethod(Enum):
     igm = 0
     truhlar = 1
     grimme = 2
+    minenkov = 3
+
+
+class _ThermoParams:
+    def __init__(self, default_sigma_r: int = 1, **kwargs: Any) -> None:
+
+        self.method = kwargs.get("lfm_method", LFMethod[Config.lfm_method])
+        if isinstance(self.method, str):
+            self.method = LFMethod[self.method.lower()]
+
+        self.temp = kwargs["temp"]
+        self.ss = kwargs.get("ss", Config.standard_state)
+        self.shift = Frequency(kwargs.get("freq_shift", Config.vib_freq_shift))
+        self.w0 = Frequency(kwargs.get("w0", Config.grimme_w0))
+        self.alpha = int(kwargs.get("alpha", Config.grimme_alpha))
+        self.sigma_r = kwargs.get("sn", default_sigma_r)
 
 
-def calculate_thermo_cont(species, temp=298.15, **kwargs):
+def calculate_thermo_cont(
+    species: "Species", temp: float = 298.15, **kwargs: Any
+):
     """
     Calculate and set the thermochemical contributions (Enthalpic, Free Energy)
     of a species using a variant of the ideal gas model (RRHO). Appends
@@ -47,6 +70,8 @@ def calculate_thermo_cont(species, temp=298.15, **kwargs):
 
     [2] J. Phys. Chem. B, 2011, 115, 14556
 
+    [3] J. Comput. Chem., 2023, 44, 1807
+
     ---------------------------------------------------------------------------
     Arguments:
         species (autode.species.Species):
@@ -56,7 +81,7 @@ def calculate_thermo_cont(species, temp=298.15, **kwargs):
 
         lfm_method (str | LFMethod): Method used to calculate the molecular
                                     entropy by treating the low frequency modes.
-                                    One of: {'igm', 'truhlar', 'grimme'}.
+                                    One of: {'igm', 'truhlar', 'grimme', 'minenkov'}.
                                     Default: Config.lfm_method
 
         ss (str): Standard state at which the molecular entropy is calculated.
@@ -80,6 +105,8 @@ def calculate_thermo_cont(species, temp=298.15, **kwargs):
     Raises:
         (KeyError | ValueError): If frequencies are not defined
     """
+    params = _ThermoParams(default_sigma_r=species.sn, temp=temp, **kwargs)
+
     if species.n_atoms == 0:
         logger.warning(
             "Species had no atoms. Cannot calculate thermochemical "
@@ -93,22 +120,12 @@ def calculate_thermo_cont(species, temp=298.15, **kwargs):
             f" no frequencies present for {species.name}."
         )
 
-    if "lfm_method" in kwargs and type(kwargs["lfm_method"]) is str:
-        kwargs["lfm_method"] = LFMethod[kwargs["lfm_method"].lower()]
-
-    S_cont = _entropy(
-        species,
-        method=kwargs.get("lfm_method", LFMethod[Config.lfm_method]),
-        temp=temp,
-        ss=kwargs.get("ss", Config.standard_state),
-        shift=kwargs.get("freq_shift", Config.vib_freq_shift),
-        w0=kwargs.get("w0", Config.grimme_w0),
-        alpha=kwargs.get("alpha", Config.grimme_alpha),
-        sigma_r=kwargs.get("sn", species.sn),
+    logger.info(
+        f"Calculating themochemistry with {params.method} at {params.temp} K"
     )
 
-    U_cont = _internal_energy(species, temp=temp)
-
+    S_cont = _entropy(species, params)
+    U_cont = _internal_energy(species, params)
     H_cont = EnthalpyCont(U_cont + SIConstants.k_b * temp, units="J").to("Ha")
 
     # Add a method string for how this enthalpic contribution was calculated
@@ -124,7 +141,7 @@ def calculate_thermo_cont(species, temp=298.15, **kwargs):
     return None
 
 
-def _thermo_method_str(species, **kwargs):
+def _thermo_method_str(species: "Species", **kwargs: Any) -> str:
     """
     Brief summary of the important methods used in evaluating the free energy
     using entropy and enthalpy methods
@@ -152,7 +169,7 @@ def _thermo_method_str(species, **kwargs):
     return string
 
 
-def _q_trans_igm(species, ss, temp):
+def _q_trans_igm(species: "Species", ss: str, temp: float) -> float:
     """
     Calculate the translational partition function using the PIB model,
     coupled with an effective volume
@@ -193,7 +210,7 @@ def _q_trans_igm(species, ss, temp):
     return q_trans
 
 
-def _q_rot_igm(species, temp, sigma_r):
+def _q_rot_igm(species: "Species", temp: float, sigma_r: int) -> float:
     """
     Calculate the rotational partition function using the IGM method. Uses the
     rotational symmetry number, default = 1
@@ -211,6 +228,10 @@ def _q_rot_igm(species, temp, sigma_r):
     if species.n_atoms == 1:
         return 1
 
+    assert species.com is not None, "Must have a COM"
+    assert species.atoms is not None, "Must have atoms"
+    assert species.moi is not None, "Must have a moment of inertia"
+
     if species.is_linear():
         com = species.com.to("m")
         i_val = sum(
@@ -236,7 +257,7 @@ def _q_rot_igm(species, temp, sigma_r):
     return temp**1.5 / sigma_r * np.sqrt(np.pi / np.prod(omega_diag))
 
 
-def _s_trans_pib(species, ss, temp):
+def _s_trans_pib(species: "Species", ss: str, temp: float) -> float:
     """
     Calculate the translational entropy using a particle in a box model
 
@@ -255,7 +276,7 @@ def _s_trans_pib(species, ss, temp):
     return SIConstants.k_b * (np.log(q_trans) + 1.0 + 1.5)
 
 
-def _s_rot_rr(species, temp, sigma_r):
+def _s_rot_rr(species: "Species", temp: float, sigma_r: int) -> float:
     """
     Calculate the rigid rotor (RR) entropy
 
@@ -280,7 +301,7 @@ def _s_rot_rr(species, temp, sigma_r):
         return SIConstants.k_b * (np.log(q_rot) + 1.5)
 
 
-def _igm_s_vib(species, temp):
+def _igm_s_vib(species: "Species", temp: float) -> float:
     """
     Calculate the entropy of a molecule according to the Ideal Gas Model (IGM)
     RRHO method
@@ -293,6 +314,7 @@ def _igm_s_vib(species, temp):
     Returns:
         (float): S_vib
     """
+    assert species.vib_frequencies is not None, "Must have frequenecies"
     s = 0.0
 
     for freq in species.vib_frequencies:
@@ -304,7 +326,9 @@ def _igm_s_vib(species, temp):
     return float(s)
 
 
-def _truhlar_s_vib(species, temp, shift_freq):
+def _truhlar_s_vib(
+    species: "Species", temp: float, shift_freq: Frequency
+) -> float:
     """
     Calculate the entropy of a molecule according to the Truhlar's method of
     shifting low frequency modes
@@ -318,17 +342,17 @@ def _truhlar_s_vib(species, temp, shift_freq):
     Returns:
         (float): S_vib in J K-1 mol-1
     """
+    assert species.vib_frequencies is not None, "Must have frequenecies"
     s = 0
 
-    if hasattr(shift_freq, "to"):
-        shift_freq = float(shift_freq.to("cm-1"))
+    shift_cm = float(shift_freq.to("cm-1"))
 
     for freq in species.vib_frequencies:
 
         # Threshold lower bound of the frequency
-        freq = max(float(freq.real), shift_freq)
+        freq_cm = max(float(freq.to("cm-1").real), shift_cm)
 
-        x = freq * Constants.c_in_cm * SIConstants.h / SIConstants.k_b
+        x = freq_cm * Constants.c_in_cm * SIConstants.h / SIConstants.k_b
         s += SIConstants.k_b * (
             ((x / temp) / (np.exp(x / temp) - 1.0))
             - np.log(1.0 - np.exp(-x / temp))
@@ -337,7 +361,16 @@ def _truhlar_s_vib(species, temp, shift_freq):
     return float(s)
 
 
-def _grimme_s_vib(species, temp, omega_0, alpha):
+def _grimme_w(omega_0: float, freq: float, alpha: int) -> float:
+    assert (
+        abs(freq) < 6000 and abs(omega_0) < 6000
+    ), "Units may be wrong - expecing cm-1"
+    return 1.0 / (1.0 + (omega_0 / freq) ** alpha)
+
+
+def _grimme_s_vib(
+    species: "Species", temp: float, omega_0: Frequency, alpha: int
+) -> float:
     """
     Calculate the entropy according to Grimme's qRRHO method of RR-HO
     interpolation in Chem. Eur. J. 2012, 18, 9955
@@ -352,15 +385,18 @@ def _grimme_s_vib(species, temp, omega_0, alpha):
     Returns:
         (float): S_vib
     """
+    assert species.vib_frequencies is not None, "Must have frequenecies"
+    assert species.moi is not None, "Must have a moment of inertia"
+
     s = 0.0
-    omega_0 = float(omega_0.to("cm-1")) if hasattr(omega_0, "to") else omega_0
+    w0 = float(omega_0.to("cm-1")) if hasattr(omega_0, "to") else omega_0
 
     # Average I = (I_xx + I_yy + I_zz) / 3.0
     b_avg = np.trace(species.moi.to("kg m^2")) / 3.0
 
     for freq in species.vib_frequencies:
 
-        omega = freq.real.to("hz")
+        omega = float(freq.real.to("hz"))
 
         mu = SIConstants.h / (8.0 * np.pi**2 * omega)
         mu_prime = (mu * b_avg) / (mu + b_avg)
@@ -375,26 +411,21 @@ def _grimme_s_vib(species, temp, omega_0, alpha):
         ) / SIConstants.h**2
         s_r = SIConstants.k_b * (0.5 + np.log(np.sqrt(factor)))
 
-        w = 1.0 / (1.0 + (omega_0 / freq) ** alpha)
+        w = _grimme_w(omega_0=w0, freq=freq, alpha=alpha)
 
         s += w * s_v + (1.0 - w) * s_r
 
     return float(s)
 
 
-def _entropy(species, method, temp, ss, shift, w0, alpha, sigma_r):
+def _entropy(species: "Species", params: _ThermoParams) -> float:
     """
     Calculate the entropy
 
     ---------------------------------------------------------------------------
     Arguments:
         species (autode.species.Species):
-        method (str):
-        temp (float):
-        ss (str):
-        shift (float)
-        w0 (float):
-        alpha (float | int):
+        params (autode.thermochemistry.igm._ThermoParams):
 
     Returns:
         (float): S in SI units
@@ -402,30 +433,35 @@ def _entropy(species, method, temp, ss, shift, w0, alpha, sigma_r):
     Raises:
         (NotImplementedError):
     """
-    logger.info(f"Calculating molecular entropy. σ_R = {sigma_r}")
+    logger.info(f"Calculating molecular entropy. σ_R = {params.sigma_r}")
+    temp = params.temp
 
     # Translational entropy component
-    s_trans = _s_trans_pib(species, ss=ss, temp=temp)
+    s_trans = _s_trans_pib(species, ss=params.ss, temp=params.temp)
 
     if species.n_atoms < 2:
         # A molecule one or no atoms has no rotational/vibrational DOF
         return s_trans
 
     # Rotational entropy component
-    s_rot = _s_rot_rr(species, temp=temp, sigma_r=sigma_r)
+    s_rot = _s_rot_rr(species, temp=temp, sigma_r=params.sigma_r)
 
     # Vibrational entropy component
-    if method == LFMethod.igm:
+    if params.method == LFMethod.igm:
         s_vib = _igm_s_vib(species, temp)
 
-    elif method == LFMethod.truhlar:
-        s_vib = _truhlar_s_vib(species, temp, shift_freq=shift)
+    elif params.method == LFMethod.truhlar:
+        s_vib = _truhlar_s_vib(species, temp, shift_freq=params.shift)
 
-    elif method == LFMethod.grimme:
-        s_vib = _grimme_s_vib(species, temp, omega_0=w0, alpha=alpha)
+    elif (
+        params.method == LFMethod.grimme or params.method == LFMethod.minenkov
+    ):
+        s_vib = _grimme_s_vib(
+            species, temp, omega_0=params.w0, alpha=params.alpha
+        )
 
     else:
-        raise NotImplementedError(f"Unrecognised method: {method}")
+        raise NotImplementedError(f"Unrecognised method: {params.method}")
 
     logger.info(
         f"S_trans = {s_trans*Constants.n_a:.3f} J K-1 mol-1\n"
@@ -437,7 +473,7 @@ def _entropy(species, method, temp, ss, shift, w0, alpha, sigma_r):
     return s_trans + s_rot + s_vib
 
 
-def _zpe(species):
+def _zpe(species: "Species"):
     """
     Calculate the zero point energy of a molecule, contributed to by the real
     (positive) frequencies
@@ -453,52 +489,63 @@ def _zpe(species):
     if species.n_atoms < 2:
         return 0.0
 
+    assert species.vib_frequencies is not None, "Must have frequenecies"
+
     zpe = 0.0
     for freq in species.vib_frequencies:
-        zpe += 0.5 * SIConstants.h * freq.real.to("hz")
+        zpe += 0.5 * SIConstants.h * float(freq.real.to("hz"))
 
     return float(zpe)
 
 
-def _internal_vib_energy(species, temp):
+def _internal_vib_energy(species: "Species", params: _ThermoParams) -> float:
     """
     Calculate the internal energy from vibrational motion within the IGM
 
     ---------------------------------------------------------------------------
     Arguments:
         species (autode.species.Species):
-        temp (float): Temperature in K
+        params (_ThermoParams):
 
     Returns:
         (float): U_vib in SI units
     """
-    e_vib = 0.0
+    assert species.vib_frequencies is not None, "Must have frequenecies"
+
+    temp = params.temp
+    w0_cm = float(params.w0.to("cm-1"))
+
+    u = 0.0  # Total internal vibrational energy
+    u_r = 0.5 * SIConstants.k_b * temp  # Free rotor internal energy
 
     # Final 6 vibrational frequencies are translational/rotational
     for freq in species.vib_frequencies:
-        x = (
-            float(freq.real.to("cm-1"))
-            * Constants.c_in_cm
-            * SIConstants.h
-            / SIConstants.k_b
-        )
-        e_vib += SIConstants.k_b * x * (1.0 / (np.exp(x / temp) - 1.0))
+        freq_cm = float(freq.real.to("cm-1"))
+        x = freq_cm * Constants.c_in_cm * SIConstants.h / SIConstants.k_b
 
-    return float(e_vib)
+        u_v = SIConstants.k_b * x * (1.0 / (np.exp(x / temp) - 1.0))
+        if params.method == LFMethod.minenkov:
+            w = _grimme_w(omega_0=w0_cm, freq=freq_cm, alpha=params.alpha)
+            u += w * u_v + (1 - w) * u_r
+        else:
+            u += u_v
 
+    return float(u)
 
-def _internal_energy(species, temp):
+
+def _internal_energy(species: "Species", params: _ThermoParams) -> float:
     """
     Calculate the internal energy of a molecule
 
     ---------------------------------------------------------------------------
     Arguments:
         species (autode.species.Species):
-        temp (float): Temperature in K
+        params (_ThermoParams):
 
     Returns:
         (float): U_cont in SI units
     """
+    temp = params.temp
     e_trns = 1.5 * SIConstants.k_b * temp
 
     if species.n_atoms < 2:
@@ -514,7 +561,7 @@ def _internal_energy(species, temp):
         e_rot = 1.5 * SIConstants.k_b * temp
 
     zpe = _zpe(species)
-    e_vib = _internal_vib_energy(species, temp=temp)
+    e_vib = _internal_vib_energy(species, params)
 
     logger.info(
         f"ZPE =     {zpe*Constants.n_a/1E3:.3f}    kJ mol-1\n"
diff --git a/doc/changelog.rst b/doc/changelog.rst
index 4e1c085..3275c06 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -6,6 +6,10 @@ Changelog
 ------
 -------
 
+Functionality improvements
+**************************
+- Adds the thermochemistry method from A. Otlyotov, Y. Minenkov in https://doi.org/10.1002/jcc.27129
+
 
 Bug Fixes
 *********
diff --git a/doc/install.rst b/doc/install.rst
index ffe063d..7d55491 100644
--- a/doc/install.rst
+++ b/doc/install.rst
@@ -5,7 +5,7 @@ Dependencies
 ------------
 **autodE** is a Python package that relies on external electronic structure codes and requires:
 
-- `Python <https://www.python.org/>`_ > v. 3.6
+- `Python <https://www.python.org/>`_ > v. 3.7
 
 - One of:
 


Test patch:
diff --git a/tests/test_thermochem.py b/tests/test_thermochem.py
index a057b0a..464952c 100644
--- a/tests/test_thermochem.py
+++ b/tests/test_thermochem.py
@@ -1,6 +1,7 @@
 import os
 import pytest
 import numpy as np
+
 from autode import Molecule, Atom, Calculation, HessianKeywords
 from autode.transition_states import TSguess, TransitionState
 from autode.thermochemistry import calculate_thermo_cont
@@ -9,7 +10,7 @@ from autode.values import Energy
 from autode.species import Species
 from autode.methods import ORCA, G09
 from . import testutils
-from autode.thermochemistry.igm import _q_rot_igm, _s_rot_rr, _entropy, _zpe
+from autode.thermochemistry.igm import _q_rot_igm, _s_rot_rr, _zpe
 
 here = os.path.dirname(os.path.abspath(__file__))
 
@@ -250,10 +251,10 @@ def test_unknown_entropy_method():
 
     h2 = Molecule(atoms=[Atom("H"), Atom("H", x=0.7)])
 
-    with pytest.raises(NotImplementedError):
-        _ = _entropy(
+    with pytest.raises(KeyError):
+        _ = calculate_thermo_cont(
             species=h2,
-            method="an_unkown_method",
+            lfm_method="an_unkown_method",
             temp=298,
             ss="1M",
             shift=100,
@@ -299,3 +300,22 @@ def test_calc_thermo_with_calc():
 
     assert os.path.exists("h2_calc_hess_orca.inp")
     assert "B3LYP " in open("h2_calc_hess_orca.inp", "r").readline()
+
+
+@testutils.work_in_zipped_dir(os.path.join(here, "data", "thermochem.zip"))
+def test_themochem_minenkov_is_close_to_grimme():
+
+    mol = Molecule("alkane_hess_orca.xyz")
+    calc = Calculation(
+        name="tmp", molecule=mol, method=orca, keywords=orca.keywords.hess
+    )
+    calc.set_output_filename("alkane_hess_orca.hess")
+
+    mol.calc_thermo(calc=calc, ss="1atm", sn=1, lfm_method="grimme")
+    G_grimme = mol.g_cont
+
+    mol.calc_thermo(calc=calc, ss="1atm", sn=1, lfm_method="minenkov")
+    G_minenkov = mol.g_cont
+
+    mad = abs(float(G_grimme.to("kcal mol-1") - G_minenkov.to("kcal mol-1")))
+    assert 1e-5 < mad < 1  # Absolute differnce should be small but non zero