API reference¶

High-level pipeline¶

ThermoScreening.thermo.api.dftbplus_thermo(atoms, temperature=298.15, pressure=101325, charge=0.0, directory=None, parameter_set='3ob', spin=None, spin_constants=None, solvent=None, solvation_param_file=None, dispersion=None, quasi_rrho=False, **kwargs)[source]¶

Run the thermo calculation with geometry optimization. Returns thermo calculation object.

Parameters:

atoms (ase.Atoms) – Atoms object. The initial geometry.
temperature (float) – The temperature in K. Default is 298.15.
pressure (float) – The pressure in Pa. Default is 101325.
charge (float) – The system charge. Default is 0.0.
directory (str, optional) – Working directory to run the DFTB+ steps in (created if needed). The geometry/Hessian/modes files are written here, so separate jobs can run without clobbering each other. Defaults to the current directory.
parameter_set (str) – Slater-Koster set to locate when slako_dir is not supplied. Defaults to "3ob".
spin (float, optional) – Spin quantum number S. Defaults to the minimum-spin electron-count guess (even -> 0, odd -> 0.5). When S > 0 the DFTB+ steps run colinear spin-polarised (so radicals are treated open-shell automatically); S = 0 keeps the restricted closed-shell calculation.
spin_constants (dict, optional) – Element -> spin-constant mapping matching the Slater-Koster set in use. Defaults to the 3ob constants; pass SPIN_CONSTANTS_MIO (or the value from resolve_parameter_set()) when running the mio set.
solvent (str, optional) – Solvent name (e.g. "water") for GBSA/ALPB implicit solvation. When set, the optimisation, energy, and Hessian all run in solution. Requires the solvent’s parameter file (thermo setup-dftb --solvent <name>). Defaults to a gas-phase calculation.
solvation_param_file (str, optional) – Explicit path to a GBSA parameter file, overriding solvent (use a method-consistent set instead of the default GFN-fit one).
dispersion (str, optional) – Dispersion correction to add to the Hamiltonian. "d3-bj" uses Grimme D3 with Becke-Johnson damping (3ob-recommended parameters). Defaults to no dispersion. Requires a DFTB+ build with D3 support.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO treatment for the vibrational entropy, which tames the entropy of low-frequency modes. Default False.
**kwargs (dict) – Additional keyword arguments to pass to the Geoopt and Hessian classes.

Returns:

The thermo calculation object.

Return type:

Thermo

Notes

Reproducibility: identical inputs on the same machine give bit-identical results (DFTB+ is thread-invariant). Across machines the optimiser stopping point can vary by ~3e-5 A, giving ~2e-9 Hartree / ~3e-4 cal/mol/K scatter – physically negligible (~1e-6 of the method error).

Solvation accuracy: with 3ob/mio the GBSA parameters are GFN-xTB-fit, so solvation free energies are only qualitative and for small neutral solutes can be positive/non-monotonic in the solvent dielectric.

ThermoScreening.thermo.api.xtb_thermo(atoms, temperature=298.15, pressure=101325, charge=0.0, spin=None, method='GFN2-xTB', directory=None, quasi_rrho=False, fmax=0.01)[source]¶

Run the thermo pipeline with GFN-xTB (tblite): optimise, compute the Hessian/frequencies, and evaluate the thermochemistry.

Needs no Slater-Koster files or spin constants – GFN-xTB parameters are built in for H-Rn and open-shell systems are handled natively via the number of unpaired electrons.

Parameters:

atoms (ase.Atoms) – Atoms object. The initial geometry.
temperature (float) – The temperature in K. Default is 298.15.
pressure (float) – The pressure in Pa. Default is 101325.
charge (float) – The system charge. Default is 0.0.
spin (float, optional) – Spin quantum number S. Defaults to the minimum-spin electron-count guess (even -> 0, odd -> 0.5). round(2*S) unpaired electrons are passed to xTB, so radicals run open-shell automatically.
method (str) – GFN-xTB parametrisation, "GFN2-xTB" (default) or "GFN1-xTB".
directory (str, optional) – Working directory to run in (created if needed). Defaults to the current directory.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO vibrational entropy. Default False.
fmax (float) – Force convergence threshold for the optimisation (eV/A). Default 0.01.

Returns:

The thermo calculation object.

Return type:

Thermo

Notes

Reproducibility: with OMP_NUM_THREADS=1 identical inputs give bit-identical energies on the same machine; the entropy is thread-invariant regardless. The optimiser stopping point is geometry/tolerance sensitive, so tighten fmax (e.g. 1e-3) for more environment-stable geometries.

ThermoScreening.thermo.api.xtb_cli_thermo(atoms, temperature=298.15, pressure=101325, charge=0.0, spin=None, method='GFN2-xTB', solvent=None, directory=None, quasi_rrho=False)[source]¶

Run the thermo pipeline with the native xtb program (GFN-xTB).

Unlike xtb_thermo() (in-process tblite, gas-phase), this uses the xtb binary, which natively supports open-shell (--uhf), charge (--chrg) and implicit solvation (--alpb) – so charged radicals in solution are handled with method-consistent parameters. A single xtb --ohess optimises and computes the Hessian.

Parameters:

atoms (ase.Atoms) – The initial geometry.
temperature (float) – Temperature in K. Default 298.15.
pressure (float) – Pressure in Pa. Default 101325.
charge (float) – System charge. Default 0.0.
spin (float, optional) – Spin quantum number S. Defaults to the minimum-spin electron-count guess (even -> 0, odd -> 0.5); round(2*S) unpaired electrons are passed to xtb, so radicals run open-shell automatically.
method (str) – GFN parametrisation: "GFN2-xTB" (default), "GFN1-xTB" or "GFN0-xTB".
solvent (str, optional) – ALPB implicit-solvation solvent (e.g. "water"). Native to xtb, so the parameters are method-consistent. Defaults to gas phase.
directory (str, optional) – Working directory (created if needed). Defaults to the current directory.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO vibrational entropy. Default False.

Returns:

The thermo calculation object.

Return type:

Thermo

ThermoScreening.thermo.api.xtb_fukui_indices(atoms, charge=0.0, spin=None, method='GFN2-xTB', solvent=None, directory=None)[source]¶

Per-atom Fukui reactivity indices via the native xtb program.

A single-point analysis (no geometry optimisation) – pass an already optimised geometry, e.g. from xtb_cli_thermo(). Each index estimates how much electron density a site would gain or lose if the whole molecule were reduced or oxidised by one electron, which is a direct, cheap way to localise where a computed reduction/oxidation is likely to happen (or where an unexpected side reaction might), rather than relying on whole- molecule quantities like total_EeGtot() or frontier orbital energies.

Parameters:

atoms (ase.Atoms) – The (already optimised) geometry to analyse.
charge (float) – System charge. Default 0.0.
spin (float, optional) – Spin quantum number S. Defaults to the minimum-spin electron-count guess (even -> 0, odd -> 0.5), matching xtb_cli_thermo().
method (str) – GFN parametrisation: "GFN2-xTB" (default), "GFN1-xTB" or "GFN0-xTB".
solvent (str, optional) – ALPB implicit-solvation solvent (e.g. "water"). Defaults to gas phase.
directory (str, optional) – Working directory (created if needed). Defaults to the current directory.

Returns:

(symbol, f_plus, f_minus, f_zero) per atom, in atoms order. f_plus is the susceptibility to nucleophilic attack (electron gain, i.e. reduction) at that site, f_minus to electrophilic attack (electron loss, oxidation), and f_zero (their average) to radical attack.

Return type:

list of tuple(str, float, float, float)

Raises:

ValueError – If method is unknown, or xtb’s output has no Fukui table.
RuntimeError – If the xtb run fails.

ThermoScreening.thermo.api.orca_thermo(hess_file, energy=None, temperature=298.15, pressure=101325, charge=0.0, spin=None, quasi_rrho=False, transition_state=False)[source]¶

Run the thermochemistry from an ORCA .hess file (DFT-quality data).

Reads the geometry, vibrational frequencies and energy from the ORCA frequency output and evaluates the RRHO thermochemistry on them, so the absolute energetics (e.g. reaction/redox free energies) are DFT-accurate rather than semiempirical.

Parameters:

hess_file (str) – Path to an ORCA .hess file.
energy (float, optional) – Electronic energy in Hartree. Defaults to the companion <basename>.property.txt file ORCA writes alongside the .hess (see ThermoScreening.calculator.orca.read_orca_hess()); pass this to override it, or when that file is unavailable.
temperature (float) – Temperature in K. Default 298.15.
pressure (float) – Pressure in Pa. Default 101325.
charge (float) – System charge. Default 0.0.
spin (float, optional) – Spin quantum number S. Defaults to the minimum-spin electron-count guess.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO vibrational entropy. Default False.
transition_state (bool) – If True, treat the file as a first-order saddle point (e.g. an ORCA OptTS + freq run): its one imaginary frequency is required and excluded from the vibrational thermochemistry instead of raising. See Thermo.imaginary_mode_wavenumber. Default False (a minimum is expected).

Returns:

The thermo calculation object.

Return type:

Thermo

Raises:

TSValueError – If no energy is available (no companion .property.txt/usable $act_energy and no energy).

ThermoScreening.thermo.api.cclib_thermo(output_file, energy=None, temperature=298.15, pressure=101325, charge=0.0, spin=None, quasi_rrho=False, transition_state=False)[source]¶

Run the thermochemistry from a QM output file via cclib.

Consumes a frequency calculation from any cclib-supported program (Gaussian, Turbomole, ORCA, Psi4, NWChem, …) and evaluates the RRHO thermochemistry on its DFT-quality geometry, frequencies and energy. Requires the optional qm extra (pip install thermoscreening[qm]).

Parameters:

output_file (str or list of str) – Path to a QM frequency-calculation output file. Turbomole splits a job’s output across many small files instead of one logfile (control, coord, aoforce.out, …); pass a list of every relevant file’s path for it (cclib’s multi-file mode) rather than a single path.
energy (float, optional) – Electronic energy in Hartree. Defaults to the best energy cclib parses from the file(s); pass this to override it (e.g. a higher-level single point) or when cclib finds no energy.
temperature (float) – Temperature in K. Default 298.15.
pressure (float) – Pressure in Pa. Default 101325.
charge (float) – System charge. Default 0.0.
spin (float, optional) – Spin quantum number S. Defaults to the minimum-spin electron-count guess.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO vibrational entropy. Default False.
transition_state (bool) – If True, treat the file as a first-order saddle point (e.g. a Gaussian TS optimization + freq run): its one imaginary frequency is required and excluded from the vibrational thermochemistry instead of raising. See Thermo.imaginary_mode_wavenumber. Default False (a minimum is expected).

Returns:

The thermo calculation object.

Return type:

Thermo

Raises:

TSValueError – If cclib is missing, the file cannot be parsed, or no energy is available.

ThermoScreening.thermo.api.pyscf_thermo(mean_field, frequencies=None, hessian=None, energy=None, temperature=298.15, pressure=101325, charge=0.0, spin=None, quasi_rrho=False, transition_state=False)[source]¶

Run the thermochemistry from an in-memory PySCF calculation.

PySCF is a Python library, so its results live in memory rather than a file. Pass a converged mean-field object; its geometry (Bohr -> Angstrom) and energy (mean_field.e_tot) are used, and the vibrational frequencies are taken from frequencies or computed from hessian via pyscf.hessian.thermo.

Parameters:

mean_field (object) – A converged PySCF mean-field object (has .mol and .e_tot).
frequencies (array_like, optional) – Vibrational frequencies in cm^-1. If omitted, hessian is used.
hessian (array_like, optional) – A PySCF Hessian (from mf.Hessian().kernel()); frequencies are derived from it. Requires the pyscf extra. Ignored if frequencies is given.
energy (float, optional) – Electronic energy in Hartree. Defaults to mean_field.e_tot.
temperature (float) – Temperature in K. Default 298.15.
pressure (float) – Pressure in Pa. Default 101325.
charge (float) – System charge. Default 0.0.
spin (float, optional) – Spin quantum number S. Defaults to the minimum-spin electron-count guess.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO vibrational entropy. Default False.
transition_state (bool) – If True, treat the geometry as a first-order saddle point: its one imaginary frequency is required and excluded from the vibrational thermochemistry instead of raising. See Thermo.imaginary_mode_wavenumber. Default False (a minimum is expected).

Returns:

The thermo calculation object.

Return type:

Thermo

Raises:

TSValueError – If neither frequencies nor hessian is given (or pyscf is missing when a hessian is used).

ThermoScreening.thermo.api.run_thermo(vibrational_frequencies, coord_file='geo_opt.xyz', temperature=298.15, pressure=101325, energy=0.0, engine='dftb+', charge=0.0, atoms=None, spin=None, quasi_rrho=False, transition_state=False)[source]¶

Run the thermo calculation. Returns thermo calculation object.

Parameters:

vibrational_frequencies (np.ndarray) – Vibrational frequencies in cm^-1.
coord_file (str) – The name of the coordinate file. Default is ‘geo_opt.xyz’. Ignored when atoms is provided.
temperature (float) – The temperature in K. Default is 298.15.
pressure (float) – The pressure in Pa. Default is 101325.
energy (float) – The energy in Hartree. Default is 0.0.
engine (str) – The name of the engine. Default is ‘dftb+’.
charge (float) – The system charge.
atoms (ase.Atoms, optional) – An optimized geometry to use directly instead of reading coord_file.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO treatment for the vibrational entropy instead of the pure harmonic oscillator. Default False.
transition_state (bool) – If True, treat the geometry as a first-order saddle point: exactly one imaginary vibrational frequency is required and excluded from the vibrational thermochemistry (see Thermo.imaginary_mode_wavenumber), instead of raising. Default False (a minimum is expected).
format. (The coordinate file should be in xyz)

Returns:

The thermo calculation object.

Return type:

Thermo

Raises:

TSValueError – If the number of vibrational frequencies does not match with the degree of freedom, or if the imaginary-mode count doesn’t match transition_state (see Thermo).

ThermoScreening.thermo.api.execute(input_file)[source]¶

Execute the thermo calculation. Returns thermo calculation object. Input file should be in the following format:

` coord_file=geo_opt.xyz vibrational_file=vibrational.out temperature=298.15 pressure=101325 energy=0.0 engine=dftb+ `

Parameters:: input_file (str) – The name of the input file.
Returns:: The thermo calculation object.
Return type:: Thermo

Screening¶

ThermoScreening.thermo.screening.screen(source, out='results', charge=0.0, temperature=298.15, pressure=101325, directory='screening', parameters=None, spin=None, parameter_set='3ob', solvent=None, dispersion=None, quasi_rrho=False, engine='dftb+', method='GFN2-xTB', resume=False, jobs=1, shard_index=None, shard_count=None)[source]¶

Run a thermochemistry screen over a set of molecules.

Parameters:

source (str) – A directory of .xyz/.gen structures (all run at charge), or a .csv manifest with path and optional name/charge/spin columns.
out (str) – Output stem; results are written to <out>.csv and <out>.json.
charge (float) – Charge used for directory input and for manifest rows without a charge.
temperature (float) – Temperature in K.
pressure (float) – Pressure in Pa.
directory (str) – Root working directory; each molecule runs in <directory>/<name>.
parameters (dict, optional) – DFTB+ Hamiltonian parameters. Defaults to those of parameter_set.
spin (float, optional) – Spin quantum number S applied to molecules without a manifest spin; defaults to an electron-count guess per molecule.
parameter_set (str) – Slater-Koster parameter set to use ("3ob" or "mio"). Selects both the default Hamiltonian parameters and the matching spin constants.
solvent (str, optional) – Solvent name for GBSA/ALPB implicit solvation applied to every molecule (e.g. "water"). Defaults to gas phase. The solvent parameter file must be installed (thermo setup-dftb --solvent <name>). With the DFTB+ engine the GBSA parameters are GFN-xTB-fit, so solvation free energies are only qualitative (and can be non-monotonic in the dielectric for small neutral solutes).
dispersion (str, optional) – Dispersion correction for the DFTB+ engine. "d3-bj" adds Grimme D3(BJ) with the 3ob-recommended parameters. Defaults to none.
quasi_rrho (bool) – If True, use Grimme’s quasi-RRHO treatment for the vibrational entropy (recommended for flexible molecules with low-frequency modes). Default False (pure harmonic oscillator).
engine (str) – Calculation engine: "dftb+" (default), "xtb" (GFN-xTB via tblite, gas-phase) or "xtb-cli" (native xtb binary, which additionally supports implicit solvation of charged radicals). parameter_set applies only to DFTB+; solvent applies to DFTB+ and xtb-cli.
method (str) – GFN-xTB parametrisation for the xtb engines ("GFN2-xTB" default).
resume (bool) – If True, reuse the successful (status="ok") records from a prior <out>.json and skip those molecules; molecules that previously failed (or were never run) are (re-)run. Results are written incrementally after every molecule regardless, so an interrupted screen can be resumed.
jobs (int) – Number of molecules to run concurrently. Default 1 (serial). With jobs > 1 the molecules run in a process pool (process-based because each job changes the working directory). Results and output ordering are unchanged; only the wall-clock time differs.
shard_index (int, optional) – Select one deterministic zero-based shard of the input for execution on a cluster. Both values are required together. Shards write to isolated output and calculation directories and can be combined with collect_screen_shards().
shard_count (int, optional) – Select one deterministic zero-based shard of the input for execution on a cluster. Both values are required together. Shards write to isolated output and calculation directories and can be combined with collect_screen_shards().

Returns:

One result record per molecule, including failed ones (status=”error”).

Return type:

list of dict

ThermoScreening.thermo.screening.collect_shards(shard_directory, out='results')[source]¶: Detect the shard workflow and dispatch to its strict collector.

ThermoScreening.thermo.screening.collect_screen_shards(shard_directory, out='results')[source]¶: Validate and combine all shards from a distributed screen.

ThermoScreening.thermo.screening.redox_screen(source, out='redox-results', charge=None, temperature=298.15, pressure=101325, directory='redox-screening', parameters=None, spin=None, parameter_set=None, solvent=None, dispersion=None, quasi_rrho=False, engine='dftb+', method=None, resume=False, jobs=1, reference=None, reference_e1=None, reference_e2=None, reference_charge=None, potential_scale=None, max_conformers=20, shard_index=None, shard_count=None)[source]¶

Run a two-step molecular reduction screen from one starting geometry.

The workflow prepares one lowest-MMFF starting conformer from SMILES when needed, computes the oxidized, one-electron-reduced and two-electron-reduced states, and reports the first, second and overall two-electron reduction potentials. Each electronic state is optimized independently, but this is not a charge-state conformer search or ensemble-free-energy calculation. A reference molecule and measured E1/E2 may calibrate both steps.

shard_index and shard_count select a deterministic candidate subset for scheduler arrays. Use collect_redox_shards() to validate and combine every shard.

ThermoScreening.thermo.screening.collect_redox_shards(shard_directory, out='redox-results')[source]¶: Validate and combine all shards from a distributed redox screen.

ThermoScreening.thermo.screening.rank_by_gibbs(results)[source]¶

Sort comparable successful records by Gibbs free energy.

Parameters:: results (list of dict) – Records as returned by screen().
Returns:: Successful records sorted by total Gibbs free energy, or an empty list when the screen mixes molecular formulas or charges. Absolute molecular energies are only comparable for equal-composition species calculated with the same settings.
Return type:: list of dict

class ThermoScreening.thermo.screening.ScreeningJob(name, path, charge, spin=None)[source]¶

A single molecule to screen.

Parameters:

name (str)
path (Path)
charge (float)
spin (float | None)

Conformer generation¶

ThermoScreening.thermo.conformers.generate(smiles, max_conformers=10, optimize=True, prune_rms_thresh=0.5, energy_window=None, random_seed=42)[source]¶

Generate conformers for a molecule from its SMILES.

Parameters:

smiles (str) – The molecule as a SMILES string.
max_conformers (int) – Maximum number of conformers to embed (ETKDG). Default 10.
optimize (bool) – If True, optimise each conformer with the MMFF force field and sort the results by MMFF energy (lowest first). Default True.
prune_rms_thresh (float) – RMSD threshold (A) for pruning duplicate embedded conformers. Default 0.5.
energy_window (float, optional) – If given (and optimize), keep only conformers within this many kcal/mol of the lowest-energy conformer.
random_seed (int) – Random seed for reproducible embedding. Default 42.

Returns:

The generated conformers (energy-sorted when optimised).

Return type:

list of ase.Atoms

Raises:

ImportError – If RDKit is not installed.
ValueError – If the SMILES cannot be parsed or no conformer could be embedded.

ThermoScreening.thermo.conformers.generate_thermo_ensemble(smiles, thermo_fn, charge=0.0, max_conformers=10, max_attempts=5, prune_rms_thresh=0.5, energy_window=None, random_seed=42, **thermo_kwargs)[source]¶

Generate a conformer ensemble and compute Thermo for each, retrying past saddle points.

A conformer that fails to converge to a true minimum (TSValueError, e.g. an imaginary frequency) or whose engine process fails outright (RuntimeError, e.g. an xtb CLI non-zero exit) is skipped; if fewer than max_conformers succeed, generate() is re-run with a new random seed, up to max_attempts seeds total. This is the loop a charged-species ensemble (radical anion, dianion, …) otherwise requires by hand before it can be passed to EnsembleThermo.

Parameters:

smiles (str) – The molecule as a SMILES string.
thermo_fn (callable) – A Thermo-computing engine, e.g. xtb_cli_thermo, called as thermo_fn(atoms, charge=charge, **thermo_kwargs).
charge (float) – System charge, forwarded to thermo_fn. Default 0.0.
max_conformers (int) – Target number of successfully-computed conformers. Default 10.
max_attempts (int) – Number of distinct random seeds (each embedding up to max_conformers conformers) to try before giving up. Default 5.
prune_rms_thresh (see generate().)
energy_window (see generate().)
random_seed (int) – Random seed for the first attempt; later attempts use random_seed + 1, random_seed + 2, etc. Default 42.
**thermo_kwargs – Forwarded to thermo_fn (e.g. solvent="water").

Returns:

The successfully computed conformers (fewer than max_conformers if max_attempts is exhausted first).

Return type:

list of Thermo

Raises:

ValueError – If no conformer converges to a true minimum within max_attempts.

ThermoScreening.thermo.conformers.write_conformers(conformers, directory, prefix='conformer')[source]¶

Write conformers to directory as <prefix>_<i>.xyz files.

Returns the list of written paths, so the directory can be passed straight to ThermoScreening.thermo.screening.screen().

Reactions and redox¶

ThermoScreening.thermo.reactions.reaction_free_energy(reactants, products, unit='kcal')[source]¶

Reaction free energy dG = sum(products) - sum(reactants).

Parameters:

reactants (iterable) – Each entry is a Thermo (stoichiometric coefficient 1) or a (coefficient, Thermo) tuple. All species should be computed with the same engine and conditions for the difference to be meaningful.
products (iterable) – Each entry is a Thermo (stoichiometric coefficient 1) or a (coefficient, Thermo) tuple. All species should be computed with the same engine and conditions for the difference to be meaningful.
unit (str) – Output unit: "kcal" (kcal/mol, default), "kJ" (kJ/mol), "eV" (per particle), or "H" (Hartree per particle).

Returns:

The reaction free energy in the requested unit.

Return type:

float

Raises:

ValueError – If unit is not supported.

ThermoScreening.thermo.reactions.reduction_potential(oxidized, reduced, n_electrons=1, reference_potential=4.44)[source]¶

Reduction potential (V) for Ox + n e- -> Red.

Uses E = -dG_red / (n F) - E_reference with dG_red = G(Red) - G(Ox) and the free energy of the electron taken as zero. Because one electron-volt per electron is one volt, the absolute potential simplifies to -dG_red[eV] / n, which is then referenced to reference_potential.

Parameters:

oxidized (Thermo) – The oxidised and reduced species, computed consistently (same engine, same solvent, and open-shell/charged as appropriate).
reduced (Thermo) – The oxidised and reduced species, computed consistently (same engine, same solvent, and open-shell/charged as appropriate).
n_electrons (int) – Positive number of electrons transferred. Default 1.
reference_potential (float) – Absolute potential (V) of the reference electrode to report against. Defaults to the SHE (SHE_ABSOLUTE_POTENTIAL); pass 0.0 for the absolute reduction potential.

Returns:

The reduction potential in volts (versus reference_potential).

Return type:

float

Raises:

ValueError – If n_electrons is not positive.

ThermoScreening.thermo.reactions.calibrate_reduction_reference(oxidized, reduced, experimental_potential, n_electrons=1)[source]¶

Calibrate a reduction-potential reference against one known redox pair.

The returned value can be passed as reference_potential to reduction_potential() for other compounds computed with the same method and conditions. This assumes that the systematic offset between the computed absolute energy scale and the experimental reference electrode is transferable from the reference pair to the target compounds.

Parameters:

oxidized (Thermo) – The reference redox pair, computed consistently with the target compounds.
reduced (Thermo) – The reference redox pair, computed consistently with the target compounds.
experimental_potential (float) – Experimental reduction potential of the reference pair, in volts.
n_electrons (int) – Positive number of electrons transferred. Default 1.

Returns:

The calibrated value for reference_potential, in volts.

Return type:

float

Raises:

ValueError – If n_electrons is not positive.

Acid dissociation (pKa)¶

ThermoScreening.thermo.pka.pKa(acid, base, temperature=298.15, reference_free_energy=-270.28)[source]¶

pKa for HA(soln) -> A-(soln) + H+(soln) (the “direct method”).

pKa = dG / (R T ln 10), with dG = G(base) - G(acid) + reference_free_energy the deprotonation free energy including the proton’s aqueous reference free energy (the proton has no electronic structure, so no engine can compute it directly).

Parameters:

acid (Thermo) – The acid (HA) and its conjugate base (A-), computed consistently (same engine, same conditions – ideally in solution, e.g. solvent="water" – and open-shell/charged as appropriate).
base (Thermo) – The acid (HA) and its conjugate base (A-), computed consistently (same engine, same conditions – ideally in solution, e.g. solvent="water" – and open-shell/charged as appropriate).
temperature (float) – Temperature in K. Should match the temperature acid/base were computed at. Default 298.15.
reference_free_energy (float) – The proton’s absolute aqueous free energy in kcal/mol. Defaults to PROTON_AQUEOUS_FREE_ENERGY_KCAL; pass a value calibrated to your method (see calibrate_proton_reference()) for quantitative accuracy.

Returns:

The (dimensionless) pKa.

Return type:

float

Notes

The raw “direct method” with a literature proton reference has several-pKa-unit systematic error even at DFT+continuum-solvent levels (e.g. Ho & Coote, Theor. Chem. Acc. 2010, 125, 3). With a semiempirical engine (GFN-xTB, DFTB) it is far worse than that – not usable at all uncalibrated: verified with real xtb-cli (GFN2-xTB + ALPB water) calculations, phenol’s raw pKa comes out at approximately -100 (vs. the experimental 9.99), an error of about 110 pKa units. This is not a solvation-model inaccuracy; semiempirical tight-binding methods do not preserve an absolute ab-initio/experimental energy scale, so their electronic energies cannot be combined directly with a literature constant derived assuming DFT/ab-initio-quality absolute energies.

calibrate_proton_reference against one reference compound absorbs this offset entirely (it is an additive constant, not molecule-specific) and is required, not merely recommended, for a semiempirical engine: in the same verification, hydroquinone’s pKa calibrated against phenol (pKa_exp = 9.99) came out at 11.00 vs. the experimental 10.35 – a good result. See the “Acid dissociation (pKa)” section of the usage docs for the full worked example.

ThermoScreening.thermo.pka.calibrate_proton_reference(acid, base, experimental_pKa, temperature=298.15)[source]¶

The reference_free_energy that reproduces experimental_pKa for a reference acid/base pair.

Solves pKa() for reference_free_energy instead of pKa. Use the result as pKa(..., reference_free_energy=...) for other acids computed the same way (same engine/conditions), ideally structurally similar to the calibration pair. With a semiempirical engine (GFN-xTB, DFTB) this is required, not just recommended for extra accuracy: the raw literature proton reference alone gives an unusable result (errors of ~100 pKa units, not just a few – see pKa()’s Notes).

Parameters:

acid (Thermo) – A reference acid/base pair with a known experimental pKa.
base (Thermo) – A reference acid/base pair with a known experimental pKa.
experimental_pKa (float) – The reference acid’s experimental pKa.
temperature (float) – Temperature in K, matching acid/base. Default 298.15.

Returns:

The calibrated reference_free_energy (kcal/mol).

Return type:

float

Transition states and kinetics¶

ThermoScreening.thermo.kinetics.eyring_rate_constant(reactants, ts, temperature=298.15, kappa=1.0)[source]¶

Eyring transition-state-theory rate constant.

k = kappa * (kB T / h) * exp(-dG-double-dagger / (R T)), with the activation free energy dG-double-dagger = G(ts) - sum(G(reactants)).

Parameters:

reactants (iterable) – The reactant species, in the same convention as ThermoScreening.thermo.reactions.reaction_free_energy(): each entry is a Thermo (stoichiometric coefficient 1) or a (coefficient, Thermo) tuple. Pass a single-entry list for a unimolecular reaction.
ts (Thermo) – The transition state, computed with transition_state=True (see ThermoScreening.thermo.thermo.Thermo).
temperature (float) – Temperature in K, matching the temperature reactants and ts were computed at. Default 298.15.
kappa (float) – Transmission coefficient (a tunneling/recrossing correction). Default 1.0 (no correction); see wigner_tunneling_correction() for a simple estimate from the transition state’s imaginary frequency.

Returns:

The TST rate constant. For a single reactant this is a first-order rate constant in s^-1. For multiple reactants (a bimolecular or higher reaction) this is the pseudo-first-order rate assuming each reactant is at its computed standard state; converting to a true bimolecular+ rate constant (e.g. L mol^-1 s^-1) requires an additional standard-state correction that is not applied here.

Return type:

float

ThermoScreening.thermo.kinetics.wigner_tunneling_correction(imaginary_wavenumber, temperature=298.15)[source]¶

Wigner’s first-order tunneling correction to a TST rate constant.

kappa = 1 + (1/24) * (h c |nu_imag| / (kB T))^2 (Wigner, Z. Phys. Chem. B 19, 203 (1932)), a small quantum correction for tunneling along the reaction coordinate. Valid only for small corrections; for deep tunneling (large kappa) use a more complete treatment (e.g. Eckart).

Parameters:

imaginary_wavenumber (float) – The transition state’s imaginary-mode wavenumber in cm^-1 (as returned by Thermo.imaginary_mode_wavenumber(); the sign is ignored).
temperature (float) – Temperature in K. Default 298.15.

Returns:

The dimensionless transmission coefficient kappa (>= 1), for use as eyring_rate_constant(..., kappa=...).

Return type:

float

Conformer ensembles¶

ThermoScreening.thermo.ensemble.boltzmann_weights(thermos, temperature=298.15)[source]¶

Boltzmann populations of a conformer ensemble at temperature (K).

Parameters:

thermos (iterable of Thermo) – Conformers of the same molecule, computed with matching settings.
temperature (float) – Temperature in kelvin. Default 298.15.

Returns:

Normalised populations (summing to 1) in the order of thermos.

Return type:

list of float

Raises:

ValueError – If thermos is empty or temperature is not positive.

ThermoScreening.thermo.ensemble.ensemble_free_energy(thermos, temperature=298.15, unit='H')[source]¶

Ensemble (Boltzmann) free energy of a conformer set.

G_ensemble = -k_B T ln( sum_i exp(-G_i / k_B T) ), which lies at or below the lowest conformer free energy by the mixing (conformational) entropy.

Parameters:

thermos (iterable of Thermo) – Conformers of the same molecule, computed with matching settings.
temperature (float) – Temperature in kelvin. Default 298.15.
unit (str) – Output unit: "H" (Hartree per particle, default), "eV" (per particle), "kcal" (kcal/mol), or "kJ" (kJ/mol).

Returns:

The ensemble free energy in the requested unit.

Return type:

float

Raises:

ValueError – If thermos is empty, temperature is not positive, or unit is not supported.

ThermoScreening.thermo.ensemble.lowest_gibbs(thermos)[source]¶

Return the conformer with the lowest absolute Gibbs free energy.

Parameters:: thermos (iterable of Thermo) – Conformers of the same molecule, computed with matching settings.
Returns:: The lowest-free-energy conformer.
Return type:: Thermo
Raises:: ValueError – If thermos is empty.

class ThermoScreening.thermo.ensemble.EnsembleThermo(thermos, temperature=298.15)[source]¶

A Thermo-like adapter exposing a conformer ensemble’s Boltzmann free energy through total_EeGtot().

pKa, calibrate_proton_reference, reduction_potential, and reaction_free_energy only ever call .total_EeGtot() on the species they’re given, so wrapping a conformer ensemble in this class lets it be passed anywhere a single Thermo is expected today – no changes needed to those helpers.

Parameters:

thermos (iterable of Thermo) – Conformers of the same species, computed with matching settings.
temperature (float) – Temperature in kelvin, matching how thermos were computed and any downstream helper it’s passed to. Default 298.15.

Raises:

ValueError – If thermos is empty or temperature is not positive.

Thermochemistry core¶

class ThermoScreening.thermo.system.System(atoms=None, charge=None, periodicity=None, pbc=None, cell=None, solvation=None, solvent=None, electronic_energy=None, vibrational_frequencies=None, spin=None)[source]¶

A class for storing the system information of a chemical system.

Parameters:

atoms (list[Atom] | None)
charge (float | None)
periodicity (bool | None)
pbc (list[bool] | bool | None)
cell (Cell | ndarray | None)
solvation (bool | None)
solvent (str | None)
electronic_energy (float | None)
vibrational_frequencies (ndarray | None)
spin (float | None)

atoms¶

A list of Atom objects.

Type:: List[Atom]

charge¶

The charge of the system.

Type:: float

periodicity¶

Whether the system has periodic boundary conditions.

Type:: bool

pbc¶

The periodic boundary conditions of the system.

Type:: List[bool]

cell¶

The cell of the system.

Type:: Cell

solvation¶

Whether the system is solvated.

Type:: bool

solvent¶

The solvent of the system.

Type:: str

electronic_energy¶

The electronic energy of the system.

Type:: float

vibrational_frequencies¶

The vibrational frequencies of the system.

Type:: np.ndarray

atom_names()[source]¶: The atom names of the system.

atomic_masses()[source]¶: A list of the atomic masses of the system.

coord()[source]¶: The atom positions of the system.

x()[source]¶: The x coordinates of the system.

y()[source]¶: The y coordinates of the system.

z()[source]¶: The z coordinates of the system.

Properties()¶

----------

number_of_atoms : int: The number of atoms in the system.

dim : int: The dimension of the system.

charge : float: The charge of the system.

dof : int: The degree of freedom of the system.

spin : float: The spin of the system.

rotational_symmetry_number : int: The symmetry number of the system.

rotational_group : str: The rotational group of the system.

spacegroup_number : int: The spacegroup of the system.

spacegroup : str: The spacegroup of the system.

center_of_mass : np.ndarray: The center of mass of the system.

periodicity : bool: The periodicity of the system.

pbc : List[bool]: The periodic boundary conditions of the system.

cell : np.ndarray: The cell of the system.

solvation : str: The solvation of the system.

solvent : str: The solvent of the system.

electronic_energy : float: The electronic energy of the system.

vibrational_frequencies : np.ndarray: The vibrational frequencies of the system.

imaginary_frequencies : np.ndarray: The imaginary vibrational frequencies of the system.

has_imaginary_frequencies : bool: Whether the system has imaginary vibrational frequencies.

mass : float: The mass of the system.

real_vibrational_frequencies : np.ndarray: The real frequency of the system.

check_frequency_length : bool: Whether the system has the correct number of vibrational frequencies.

Examples

>>> from ThermoScreening.thermo.system import System
>>> from ThermoScreening.thermo.atoms import Atom
>>> from ThermoScreening.thermo.cell import Cell
>>> import numpy as np
>>> atom1 = Atom("H", [0.0, 0.0, 0.0], 1.0)
>>> atom2 = Atom("H", [0.0, 0.0, 1.0], 1.0)
>>> atoms = [atom1, atom2]
>>> cell = Cell([10.0, 10.0, 10.0])
>>> system = System(atoms=atoms, cell=cell, vibrational_frequencies=np.array([1.0, 2.0, 3.0]))
>>> system.atoms
[Atom(symbol='H', position=[0.0, 0.0, 0.0], mass=1.0), Atom(symbol='H', position=[0.0, 0.0, 1.0], mass=1.0)]

class ThermoScreening.thermo.thermo.Thermo(temperature, pressure, system, engine, quasi_rrho=False, transition_state=False)[source]¶

A class for calculating the thermochemistral properties e.g. Gibbs free energy, enthalpy, entropy, heat capacity, etc. of a system. The class needs as input the temperature, pressure, system information (coordinates, degree of freedom, electronic spin, symmetry number, electronic energy and vibrational frequencies).

Parameters:

temperature (float)
pressure (float)
system (System)
engine (str)
quasi_rrho (bool)
transition_state (bool)

temperature¶

The temperature of the system in Kelvin.

Type:: float

pressure¶

The pressure of the system in Pa.

Type:: float

system¶

The system information of the system.

Type:: System

engine¶

The engine used for the calculation to compute the thermochemical properties with correct units.

Type:: str

run()[source]¶: Calculates the thermochemical properties of the system.

electronic_energy()[source]¶

Returns the electronic energy of the system.

Returns:: The electronic energy of the system in Hartree per particle.
Return type:: float

imaginary_mode_wavenumber()[source]¶

The transition state’s imaginary-mode wavenumber (cm^-1, negative).

Set by run() when this object was constructed with transition_state=True; useful for a tunneling correction (see ThermoScreening.thermo.kinetics.wigner_tunneling_correction()).

Returns:: The imaginary mode’s wavenumber, or None for a non-transition- state Thermo (the default) or before run() has been called.
Return type:: float or None

temperature_scan(temperatures)[source]¶

Recompute the thermochemistry at each of temperatures (in Kelvin).

The electronic energy, geometry, and vibrational frequencies are temperature-independent, so this reuses the same System (and this object’s pressure, engine, quasi-RRHO and transition-state setting) and only the thermal terms are recomputed – i.e. a temperature scan from a single Hessian. This object is left unchanged.

Parameters:: temperatures (iterable of float) – Temperatures in Kelvin.
Returns:: A Thermo computed (run() already called) for each temperature, in the given order.
Return type:: list of Thermo

total_EeGtot()[source]¶

Returns the sum of the electronic energy and the total Gibbs free energy.

Returns:: The sum of the electronic energy and the total Gibbs free energy in Hartree per particle.
Return type:: float

total_energy(unit='H')[source]¶

Returns the total energy of the system based on the unit.

Parameters:: unit (str) – The unit of the total energy ("H", "kcal" or "cal"). Defaults to "H" (Hartree).
Raises:: TSValueError – If the unit is not supported.
Returns:: The total energy of the system.
Return type:: float

total_enthalpy(unit='H')[source]¶

Returns the total enthalpy of the system based on the unit.

Parameters:: unit (str) – The unit of the total enthalpy ("H", "kcal" or "cal"). Defaults to "H" (Hartree).
Raises:: TSValueError – If the unit is not supported.
Returns:: The total enthalpy of the system.
Return type:: float

total_entropy(unit='cal/(mol*K)')[source]¶

Returns the total entropy of the system based on the unit.

Parameters:: unit (str) – The unit of the total entropy ("cal/(mol*K)" or "H/T"). Defaults to "cal/(mol*K)".
Raises:: TSValueError – If the unit is not supported.
Returns:: The total entropy of the system.
Return type:: float

total_gibbs_free_energy(unit='H')[source]¶

Returns the total Gibbs free energy of the system based on the unit.

Parameters:: unit (str) – The unit of the total Gibbs free energy ("H", "kcal" or "cal"). Defaults to "H" (Hartree).
Raises:: TSValueError – If the unit is not supported.
Returns:: The total Gibbs free energy of the system.
Return type:: float

total_heat_capacity(unit='cal/(mol*K)')[source]¶

Returns the constant-volume (ideal-gas) heat capacity Cv of the system.

Cv = translational (3/2 R) + rotational (n_rot/2 R) + vibrational + electronic (0). For the constant-pressure heat capacity of an ideal gas add R: Cp = Cv + R. (Cross-checked against native xTB, whose reported heat capacity is Cp: tool Cv == xtb Cp - R.)

Parameters:: unit (str) – The unit of the total heat capacity ("cal/(mol*K)" or "H/T"). Defaults to "cal/(mol*K)".
Raises:: TSValueError – If the unit is not supported.
Returns:: The total constant-volume heat capacity Cv of the system.
Return type:: float

Coordinate and frequency readers¶

ThermoScreening.thermo.api.read_xyz(coord_file)[source]¶

Reads xyz-file and returns number of atoms, chemical symbols of atoms, coordinates of atoms, and cell parameters.

Parameter:¶

coord_filestr: The name of the coordinate file.

Returns:¶

List: number of atoms, chemical symbols of atoms, coordinates of atoms, cell parameters

Parameters:: coord_file (str)

ThermoScreening.thermo.api.read_gen(coord_file)[source]¶

Reads gen-file and returns number of atoms, chemical symbols of atoms, coordinates of atoms, and cell vector.

Parameter:¶

coord_filestr: The name of the coordinate file.

Returns:¶

List: number of atoms, chemical symbols of atoms, coordinates of atoms, cell vector

Parameters:: coord_file (str)

ThermoScreening.thermo.api.read_coord(coord_file, engine)[source]¶

Reads the coordinate file and returns the coordinates as a numpy array.

Parameters:

coord_file (str) – The name of the coordinate file.
engine (str) – The name of the engine.

Returns:

number of atoms, chemical symbols of atoms, coordinates of atoms

Return type:

List

Raises:

TSNotImplementedError – If the engine is not supported. If the input file is not supported. If the gen file is not tested yet.

ThermoScreening.thermo.api.read_vibrational(vibrational_file, engine)[source]¶

Reads the vibrational file and returns the vibrational frequencies as a numpy array.

Parameters:

vibrational_file (str) – The name of the vibrational file.
engine (str)

Returns:

The vibrational frequencies as a numpy array.

Return type:

np.ndarray

Raises:

TSNotImplementedError – If the engine is not supported.

ThermoScreening.calculator.orca.read_orca_hess(path)[source]¶

Read geometry, vibrational frequencies and energy from an ORCA .hess file.

Parameters:

path (str) – Path to an ORCA .hess file (from an ORCA frequency calculation).

Returns:

atoms (ase.Atoms) – The geometry, converted from Bohr to Angstrom.
frequencies (np.ndarray) – The vibrational frequencies in cm^-1, including the near-zero translational/rotational modes as ORCA writes them.
energy (float or None) – The electronic energy in Hartree, preferably from the companion <basename>.property.txt file ORCA writes alongside the .hess (its $Single_Point_Data / &FinalEnergy entry); falls back to a nonzero $act_energy in the .hess itself (populated for relaxed surface scans); None if neither is available.

Notes

$act_energy (along with $act_atom/$act_coord) is an ORCA relaxed-surface-scan field for the active point on the scan; for an ordinary optimization + frequency job it is 0.0 and does not carry the electronic energy at all, so it is not read unless it is genuinely nonzero.

The per-atom mass column is not used; ASE’s standard isotope masses are used for the rotational/translational terms (as with the other engines).

Raises:: TSValueError – If the required $atoms or $vibrational_frequencies block is missing or malformed.

ThermoScreening.calculator.qm.read_cclib(path)[source]¶

Read geometry, vibrational frequencies and energy from a QM output file.

Uses cclib to parse any supported program’s output (Gaussian, ORCA, Psi4, NWChem, …). Most programs write one logfile per job; pass its path.

Turbomole splits a job’s output across many small files instead of one logfile (control, coord, aoforce.out, …); pass a list of every relevant file’s path (cclib’s own multi-file mode) rather than a single path.

Parameters:

path (str or list of str) – Path to a QM frequency-calculation output file, or (for a multi-file program such as Turbomole) a list of paths to every relevant file from the same job.

Returns:

atoms (ase.Atoms) – The final geometry (cclib reports coordinates in Angstrom).
frequencies (np.ndarray) – The vibrational frequencies in cm^-1 (cclib gives the real vibrational modes; imaginary modes appear as negatives).
energy (float or None) – The electronic energy in Hartree (converted from cclib’s eV), or None if the output has no parseable energy.

Notes

cclib returns only the real vibrational modes (3N-6, or 3N-5 for a linear molecule). The thermochemistry keeps the top dof modes for the geometry’s own linearity classification, so this matches for clearly (non)linear molecules; a near-linear geometry where the classification disagrees with the QM program’s mode count could drop or miss a mode.

Raises:: TSValueError – If cclib is not installed, the file(s) cannot be parsed (including a parser crash on an unusual output format, wrapped from cclib), or there are no vibrational frequencies.

Backend setup¶

ThermoScreening.cli.dftb_setup.install_slakos(install_root=None, url=None, force=False, parameter_set='3ob-3-1')[source]¶

Download and extract a Slater-Koster parameter set.

Parameters:

install_root (str or Path, optional) – Directory where parameter sets are installed. Defaults to the user-local share directory.
url (str, optional) – Archive URL. Defaults to the release URL for parameter_set.
force (bool) – Re-download even when the set already exists.
parameter_set (str) – Set to install ("3ob"/"3ob-3-1" or "mio"/"mio-1-1").

Return type:

Path

ThermoScreening.cli.dftb_setup.install_gbsa_param(solvent, install_root=None, url=None, force=False)[source]¶

Download the GBSA implicit-solvation parameter file for solvent.

Parameters:

solvent (str) – Solvent name (e.g. "water"); see GBSA_SOLVENTS.
install_root (str or Path, optional) – Directory root for downloaded parameters. Defaults to the user-local share directory.
url (str, optional) – Explicit file URL override. Defaults to the grimme-lab release file.
force (bool) – Re-download even when the file already exists.

Return type:

Path

ThermoScreening.cli.dftb_setup.check_dftb_setup(env=None, parameter_set='3ob')[source]¶

Check whether the calculation backends are available.

Reports the DFTB+ toolchain (dftb+, modes, DFTB_PREFIX + a Slater-Koster file) as required, and the xTB toolchain (xtb binary, tblite) as optional.

Parameters:

env (dict[str, str] | None)
parameter_set (str)

Return type:

list[Diagnostic]