feature/quasi-RRHO

.gitignore

@@ -7,5 +7,6 @@ firecode.egg-info/
 
 firecode/calculators/uma-s-1p1.pt
 firecode/calculators/uma-s-1p2.pt
+firecode/calculators/esen_sm_conserving_all.pt
 coverage.xml
 .coverage

CHANGELOG.md

@@ -6,6 +6,14 @@
 <!-- FINALSPLEVEL keyword? -->
  <!-- add documentation: SCRAMBLECHECK kw, fsm>, neb>, rdkit_search>, standalone optimizer, non-inline constraints-->
 
+
+## FIRECODE 1.7.0 🔥 (WIP)
+- Restructured and expanded solvent module.
+- Implemented vibrational analysis via ASE including quasi-RRHO thermochemistry (as in GoodVibes) to get free energies.
+- Added FREQ keyword, that computes free energies at the end of the geometry optimization cycles.
+- Added T (temperature, in K), T_C (temperature, in °C), P (pressure, in atm), and C (concentration, in mol/L) keywords.
+<!-- add number of active constraints printout for parallel multithread functions -->
+
 ## FIRECODE 1.6.0 🔥 (March 11 2026)
 - Refreshed constraints handling in operators.
 - Added non-inline constraints specification after a molecule line. Syntax: "B/A/D i1 i2 [i3] [i4] [auto/value]"

docs/operators_keywords.rst

@@ -23,26 +23,26 @@ Here is a list of the currently available operators:
    and will constrain the specified distances (both "fixed" and not). Generates a new ``mol_crest_confs.xyz``
    file with the crest-optimized conformers. The default level is GFN2//GFN-FF (see CREST docs).
    Charge and multiplicity are passed from the molecule or keyword line.
-   
+
    ::
-   
+
        mtd> molecule.xyz 4A 8A charge=-1
 
 -  ``rdkit_search>`` - Performs a conformational search with the `ETKDGv3 algorithm <https://pubs.acs.org/doi/10.1021/acs.jcim.0c00025>`__
    via RDKit. Generates a new ``mol_rdkit_confs.xyz`` file with the generated conformers.
 
 -  ``refine>`` - Acts on a multimolecular input file, pruning similar structures and optimizing it
-   at the chosen theory level. Soft, non-fixed constraints ("interactions") are enforced during a 
+   at the chosen theory level. Soft, non-fixed constraints ("interactions") are enforced during a
    first loose optimization, and subsequently relaxed.
 
 -  ``neb>`` - Runs a NEB (Nudged Elastic Band) procedure on the provided structures.
-   The implementation is a climbing image nudged elastic band (CI-NEB) TS search.  
+   The implementation is a climbing image nudged elastic band (CI-NEB) TS search.
    The operator should be used on a single file that contains two or more structures.
    The first and last structures will be used as start and endpoints. If three are provided, the intermediate
    structure will be used as a TS guess in the interpolation. If more are provided, they will all be used
    in place of (or to aid) interpolation to reach the desired number of images
    (set with the ``IMAGES=n`` or NEB(IMAGES=n) keywords). Default is 7.
-   
+
 -  ``scan>`` - Runs a scan of the desired distance or dihedral angle, based on the number of indices provided.
    For distance scans, the atomic indices will be approached if they are not bound, and separated if they are.
    This could be chained with, for example, a NEB calculation, where the scan trajectory will be used as a starting
@@ -53,8 +53,8 @@ Here is a list of the currently available operators:
        neb> scan> mol.xyz 3 4
 
    For dihedral scans, the dihedral is rotated the full 360° in both clockwise and counterclockwise directions.
-   Maxima above a certain threshold (set with ``KCAL=n``) are scanned again in smaller steps. All maxima for the second mode accurate 
-   scan are saved to firecode_maxima_mol.xyz. The optimized geometries for the initial scans are saved with the relative 
+   Maxima above a certain threshold (set with ``KCAL=n``) are scanned again in smaller steps. All maxima for the second mode accurate
+   scan are saved to firecode_maxima_mol.xyz. The optimized geometries for the initial scans are saved with the relative
    names. This version of the ``scan>`` operator is terminating and cannot be chained.
 
 Deprecated (legacy):
@@ -80,16 +80,16 @@ General keywords
 
 Keywords are case-insensitive and are divided by at least one blank space.
 Some of them are self-sufficient (*i.e.* ``NCI``), while some others require an
-additional input (*i.e.* ``STEPS=10`` or ``DIST(a=1.8, b=2, c=1.34)``). 
+additional input (*i.e.* ``STEPS=10`` or ``DIST(a=1.8, b=2, c=1.34)``).
 Commas are used to divide keyword arguments where more than
 one is accepted, like in ``DIST``.
 
--  **CALC** - Overrides default calculator in ``settings.py``. 
+-  **CALC** - Overrides default calculator in ``settings.py``.
    Available calculators are ORCA, XTB, TBLITE, AIMNET2 and UMA. Syntax: ``CALC=ORCA``
 
 -  **CHARGE** - Specify the charge to be used in optimizations.
 
--  **CONFS** - Override the maximum number of conformers to be considered 
+-  **CONFS** - Override the maximum number of conformers to be considered
    for the refinement or embedding of each molecule (default is 1000). Syntax: ``CONFS=10000``
 
 -  **CRESTNCI** - crest_search> runs: passes the "--nci" argument to CREST, running
@@ -114,7 +114,7 @@ one is accepted, like in ``DIST``.
    Value can only be ``XTB`` or None. Syntax: ``FFCALC=XTB``
 
 -  **FFLEVEL** - Manually set the theory level to be used for force field
-   calculations. Default is GFN-FF for XTB. Support for other FF calculators has been discontinued.  
+   calculations. Default is GFN-FF for XTB. Support for other FF calculators has been discontinued.
    Standard values can be modified by running the module with the -s flag
    (>>> python -m firecode -s).
 
@@ -132,7 +132,7 @@ one is accepted, like in ``DIST``.
 
 -  **LEVEL** - Manually set the theory level to be used.
    White spaces, if needed, can be expressed in input files with underscores. Syntax (ORCA):
-   ``LEVEL(B3LYP_def2-TZVP)``. Defaults can be found in settings.py, and can be modified by running 
+   ``LEVEL(B3LYP_def2-TZVP)``. Defaults can be found in settings.py, and can be modified by running
    the module with the -s flag (``>>> firecode -s``).
 
 -  **NEB** - Set NEB options. Syntax and default values: ``NEB(images=7, preopt=true, ci=true)``.
@@ -141,6 +141,8 @@ one is accepted, like in ``DIST``.
    refine bonding distances. Set by default with the ``SHRINK`` keyword
    and for monomolecular TSs.
 
+-  **P** - Pressure, in atm (default 1.0). Syntax: ``P=100.0``
+
 -  **PKA** - Specify the reference pKa for a compound in multimolecular
    pKa calculation runs. Syntax: ``PKA(mol.xyz)=11``
 
@@ -149,21 +151,25 @@ one is accepted, like in ``DIST``.
    ``settings.py``. Syntax: ``PROCS=32``
 
 -  **RMSD** - RMSD threshold (Angstroms) for structure pruning.
-   The smaller, the more retained structures (default is 0.25 A). 
+   The smaller, the more retained structures (default is 0.25 A).
    Syntax: ``THRESH=n``
 
+-  **T** - Temperature, in Kelvin (default 298.15). Syntax: ``T=300.0``
+
+-  **T_C** - Temperature, in Celsius (default 25.0). Syntax: ``T_C=25.0``
+
 
 Embedding-specific keywords (legacy)
 ++++++++++++++++++++++++++++++++++++
 
 -  **BYPASS** - Debug keyword. Used to skip all pruning steps and
    directly output all the embedded geometries.
-   
+
 -  **CLASHES** - Manually specify the max number of clashes and/or
    the distance threshold at which two atoms are considered clashing.
    The more forgiving (higher number, smaller dist), the more structures will reach the geometry
    optimization step. Default values are num=0 and dist=1.5 (A). Syntax: ``CLASHES(num=3,dist=1.2)``
-   
+
 -  **DEEP** - Performs a deeper search, retaining more starting
    points for calculations and smaller turning angles. Equivalent to
    ``THRESH=0.3 STEPS=72 CLASHES=(num=1,dist=1.3)``. **Use with care!**
@@ -172,18 +178,18 @@ Embedding-specific keywords (legacy)
    (C=C and C=N) during the embed. Likely to be useful only for
    monomolecular embeds, where molecular distortion is often important, and
    undesired isomerization processes can occur.
- 
+
 -  **NOOPT** - Skip the optimization steps, directly writing
    structures to file after compenetration and similarity pruning.
    Dihedral embeds: performs rigid scans instead of relaxed ones.
-   
+
 -  **RIGID** - Only applies to "cyclical"/"chelotropic" embeds.
    Avoid bending structures to better build TSs.
 
 -  **ROTRANGE** - Only applies to "cyclical"/"chelotropic" embeds.
    Manually specify the rotation range to be explored around the
    structure pivot. Default is 45. Syntax: ``ROTRANGE=90``
-   
+
 -  **SHRINK** - Exaggerate orbital dimensions during embed, scaling
    them by a specified factor. If used as a single keyword (``SHRINK``),
    orbital dimensions are scaled by a factor of one and a half. A syntax
@@ -206,9 +212,9 @@ Embedding-specific keywords (legacy)
    default ``STEPS=24`` will perform 15° turns. For "cyclical" and
    "chelotropic" embeds, the rotation range to be explored is
    +-\ ``ROTRANGE`` degrees. Therefore the default values, equivalent to
-   ``ROTRANGE=45 STEPS=5``, will sample five equally spaced positions between 
+   ``ROTRANGE=45 STEPS=5``, will sample five equally spaced positions between
    +45 and -45 degrees (going through zero).
 
 -  **SUPRAFAC** - Only retain suprafacial orbital configurations in
    cyclical TSs. Thought for Diels-Alder and other cycloaddition
-   reactions.
+   reactions.

firecode/__main__.py

@@ -14,7 +14,7 @@
 
 https://github.com/ntampellini/firecode
 
-Nicolo' Tampellini - nicolo.tampellini@yale.edu
+Nicolo' Tampellini - ntamp@mit.edu
 
 """
 
@@ -29,7 +29,6 @@
 
 
 def main() -> None:
-
     # Redirect stdout and stderr to handle encoding errors
     sys.stdout = TextIOWrapper(sys.stdout.buffer, encoding="utf-8", errors="replace")
 

firecode/ase_manipulations.py

@@ -26,7 +26,6 @@
 import time
 from dataclasses import dataclass
 from functools import wraps
-from shutil import rmtree
 from subprocess import getoutput
 from typing import TYPE_CHECKING, Any, Callable, Iterable, ParamSpec, Sequence, TypeVar, cast
 
@@ -37,8 +36,6 @@
 from ase.constraints import FixInternals
 from ase.mep import DyNEB
 from ase.optimize import FIRE, LBFGS
-from ase.thermochemistry import IdealGasThermo
-from ase.vibrations import Vibrations
 from prism_pruner.algebra import dihedral, normalize
 from prism_pruner.graph_manipulations import d_min_bond, find_paths
 from prism_pruner.rmsd import rmsd_and_max
@@ -48,13 +45,12 @@
 from firecode.calculators._xtb import xtb_gsolv
 from firecode.settings import DEFAULT_LEVELS
 from firecode.typing_ import Array1D_float, Array1D_str, Array2D_float, Array3D_float, MaybeNone
-from firecode.units import EH_TO_EV, EH_TO_KCAL, EV_TO_KCAL
+from firecode.units import EV_TO_KCAL
 from firecode.utils import (
     HiddenPrints,
     NewFolderContext,
     cartesian_product,
     read_xyz,
-    suppress_stdout_stderr,
     write_xyz,
 )
 
@@ -590,7 +586,8 @@ def ase_neb(
                 opt.run(fmax=0.05, steps=200 + opt.nsteps)
 
             energies = [image.get_total_energy() * EV_TO_KCAL for image in images]  # type: ignore[no-untyped-call]
-            print(f"--> Updated temporary MEP at {title}_MEP_temp_pre_CI.xyz")
+            if verbose_print:
+                print(f"--> Updated temporary MEP at {title}_MEP_temp_pre_CI.xyz")
             ase_dump(f"{title}_MEP_temp_pre_CI.xyz", atoms, neb.images, energies)
 
             if climbing_image:
@@ -867,7 +864,6 @@ def ase_popt(
         ase_calc = Opt_func_dispatcher(calculator).get_ase_calc(method, solvent)
 
     ase_atoms.calc = ase_calc
-    ase_calc.verbosity = 0  # type: ignore[union-attr]
 
     ase_atoms = set_charge_and_mult_on_ase_atoms(ase_atoms, charge, mult)
     maxiter = maxiter or 750
@@ -929,7 +925,7 @@ def ase_popt(
     with NewFolderContext(title, delete_after=(not debug)):
         t_start_opt = time.perf_counter()
 
-        with suppress_stdout_stderr():
+        with HiddenPrints():
             with LBFGS(ase_atoms, maxstep=0.2, trajectory=traj) as opt:
                 opt.run(fmax=fmax, steps=maxiter)
                 iterations = opt.nsteps
@@ -1034,91 +1030,6 @@ def is_linear(atoms: Atoms, tol: float = 1e-3) -> bool:
     return cast("bool", s[1] < tol)
 
 
-def ase_vib(
-    atoms: Array1D_str,
-    coords: Array2D_float,
-    ase_calc: ASECalculator,
-    charge: int,
-    mult: int,
-    temp: float = 298.15,
-    title: str = "temp",
-    return_gcorr: bool = True,
-    write_log: bool = True,
-) -> float:
-    """returns: G(corr) or Free energy, in kcal/mol."""
-    ase_atoms = Atoms(atoms, positions=coords)
-    ase_atoms = set_charge_and_mult_on_ase_atoms(ase_atoms, charge=charge, mult=mult)
-
-    ase_calc.verbosity = 0  # type: ignore[attr-defined]
-    ase_atoms.calc = ase_calc
-
-    vib = Vibrations(ase_atoms, delta=0.005)  # type: ignore[no-untyped-call]
-
-    with open(f"vib_{title}.out", "w", encoding="utf-8") as f:
-        with HiddenPrints():
-            f.write("--> FIRECODE ASE Frequency calculation report\n\n")
-
-            # tighten convergence
-            f.write("--> Tightening geom. opt. convergence to fmax=1E-4\n")
-            opt = LBFGS(ase_atoms, maxstep=0.01)
-            opt.run(fmax=1e-4)  # type: ignore[no-untyped-call]
-            energy = ase_atoms.get_potential_energy()  # type: ignore[no-untyped-call]
-
-            # remove cache folder
-            if "vib" in os.listdir():
-                rmtree(os.path.join(os.getcwd(), "vib"))
-
-            # run vibrational analysis
-            vib.run()  # type: ignore[no-untyped-call]
-
-            # get energies (frequencies) and print summary
-            vib_energies = vib.get_energies()  # type: ignore[no-untyped-call]
-            vib.summary(log=f)  # type: ignore[no-untyped-call]
-
-            if len(atoms) == 1:
-                geometry = "monatomic"
-            elif is_linear(ase_atoms):
-                geometry = "linear"
-            else:
-                geometry = "nonlinear"
-
-            # compute thermo
-            thermo = IdealGasThermo(  # type: ignore[no-untyped-call]
-                vib_energies=vib_energies,
-                potentialenergy=energy,
-                atoms=ase_atoms,
-                spin=(mult - 1) / 2,
-                geometry=geometry,
-                symmetrynumber=1,
-                ignore_imag_modes=True,
-            )
-
-            EE = energy / EH_TO_EV  # was eV, now Eh
-            G = thermo.get_gibbs_energy(temperature=temp, pressure=101325.0) / EH_TO_EV  # type: ignore[no-untyped-call]
-            H = thermo.get_enthalpy(temperature=temp) / EH_TO_EV  # type: ignore[no-untyped-call]
-            S = thermo.get_entropy(temperature=temp, pressure=101325.0) / EH_TO_EV  # type: ignore[no-untyped-call]
-            gcorr = G - EE
-
-            f.write("\n--> What follows mocks an ORCA output for scraping purposes:\n\n")
-            f.write(f"T: {temp:.2f} K ({temp - 273.15:.2f} °C)\n")
-            f.write(f"FINAL SINGLE POINT ENERGY {EE:.8f} Eh\n")
-            f.write(f"FINAL GIBBS FREE ENERGY {G:.8f} Eh\n")
-            f.write(f"G-E(el) ... {gcorr:.8f} Eh     {gcorr * EH_TO_KCAL:.2f} kcal/mol\n")
-            f.write(f"Total enthalpy ... {H:.8f} Eh\n")
-            f.write(f"Final entropy term ... {S:.8f} Eh/K\n")
-
-    del vib
-    if os.path.isdir("vib"):
-        rmtree("vib")
-
-    if not write_log:
-        os.remove(f"vib_{title}.out")
-
-    if return_gcorr:
-        return cast("float", gcorr * EH_TO_KCAL)
-    return cast("float", G * EH_TO_KCAL)
-
-
 def fsm_operator(embedder: Embedder, optimize_endpoints: bool = True) -> str:
     """Connect two structures using the Freezing String Method and return the resulting list of coordinates."""
     assert len(embedder.objects) == 1, "FSM requires a single input structure."