Add transition-state thermochemistry and Eyring rate constants - #105
Merged
Conversation
Thermo.run() previously raised on any imaginary vibrational frequency, so a transition state (a first-order saddle point) could not be evaluated at all. Add Thermo(..., transition_state=True): requires exactly one imaginary mode among the kept frequencies, excludes it from the vibrational partition function/entropy/energy/heat-capacity sums (the reaction coordinate contributes no thermal vibrational term), and exposes its wavenumber via Thermo.imaginary_mode_wavenumber(). Non-TS behavior (the default) is unchanged; temperature_scan propagates the flag. Thread transition_state=False through run_thermo, orca_thermo, cclib_thermo and pyscf_thermo -- the "import an externally-computed structure" engines, since a DFTB+/xtb geometry optimization is a minimizer and cannot land on a saddle point. Add ThermoScreening.thermo.kinetics: - eyring_rate_constant(reactants, ts, temperature, kappa=1.0): TST rate constant from the activation free energy, using the same Thermo/(coefficient, Thermo) stoichiometry convention as reaction_free_energy. - wigner_tunneling_correction(imaginary_wavenumber, temperature): a simple first-order tunneling estimate from the TS's imaginary frequency. Closes #104
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #105 +/- ##
==========================================
+ Coverage 97.22% 97.25% +0.03%
==========================================
Files 30 31 +1
Lines 1872 1895 +23
==========================================
+ Hits 1820 1843 +23
Misses 52 52
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Thermo.run()previously raised on any imaginary vibrational frequency, so a transition state (a first-order saddle point) could not be evaluated at all — blocking barrier/rate work.Core:
Thermo(..., transition_state=True)Thermo.imaginary_mode_wavenumber()(cm⁻¹, negative).temperature_scanpropagates the flag.Plumbing
transition_state=Falsethreaded throughrun_thermo,orca_thermo,cclib_thermo,pyscf_thermo— the "import an externally-computed structure" engines, since a DFTB+/xtb geometry optimization is a minimizer and can't itself land on a saddle point.New:
ThermoScreening.thermo.kineticseyring_rate_constant(reactants, ts, temperature, kappa=1.0)— TST rate constantk = kappa * (kB T/h) * exp(-ΔG‡/RT), using the sameThermo/(coefficient, Thermo)stoichiometry convention asreaction_free_energy.wigner_tunneling_correction(imaginary_wavenumber, temperature)— a simple first-order tunneling estimate from the TS's imaginary frequency.Verification
k(ΔG‡=0) = kB*T/h(~6.212×10¹² s⁻¹ at 298.15 K, the exact textbook TST prefactor) and cross-checked a 20 kcal/mol barrier against the well-known "ΔG‡≈17.5 kcal/mol → k≈1 s⁻¹" rule of thumb — both match exactly.temperature_scanpropagation),kinetics.py(closed-form cross-checks, stoichiometry, kappa scaling), and one transition-state forwarding test per engine (run_thermo,orca_thermo,cclib_thermo,pyscf_thermo).Closes #104