Source code for pysisyphus.thermo

import h5py
import jinja2
import numpy as np

try:
    from thermoanalysis.QCData import QCData
    from thermoanalysis.thermo import thermochemistry

    can_thermoanalysis = True
except ModuleNotFoundError:
    can_thermoanalysis = False

from pysisyphus.config import p_DEFAULT, T_DEFAULT
from pysisyphus.constants import AU2KJPERMOL, AU2KCALPERMOL, BOHR2ANG
from pysisyphus.helpers_pure import highlight_text
from pysisyphus.Geometry import Geometry


[docs] def get_thermoanalysis_from_hess_h5( h5_fn, T=T_DEFAULT, p=p_DEFAULT, point_group="c1", return_geom=False, kind="qrrho", scale_factor=1.0, ): with h5py.File(h5_fn, "r") as handle: masses = handle["masses"][:] vibfreqs = handle["vibfreqs"][:] coords3d = handle["coords3d"][:] energy = handle.attrs["energy"] mult = handle.attrs["mult"] atoms = handle.attrs["atoms"] thermo_dict = { "masses": masses, "wavenumbers": vibfreqs, # thermoanalysis expects coordinates in Angstrom "coords3d": coords3d * BOHR2ANG, "scf_energy": energy, "mult": mult, } qcd = QCData(thermo_dict, point_group=point_group) thermo = thermochemistry( qcd, temperature=T, pressure=p, kind=kind, scale_factor=scale_factor, ) if return_geom: geom = Geometry(atoms=atoms, coords=coords3d) return thermo, geom else: return thermo
THERMO_TPL = jinja2.Template( """ {% if geom -%} Geometry : {{ geom }}, {{ geom.atoms|length }} atoms {%- endif %} Temperature : {{ "%0.2f" % thermo.T }} K Pressure : {{ thermo.p }} Pa Total Mass : {{ "%0.4f" % thermo.M }} amu ! Symmetry is currently not supported in pysisyphus. ! ! If not given otherwise, c1 and σ = 1 are assumed. ! Point Group : {{ thermo.point_group }} Symmetry Number σ : {{ thermo.sym_num }} Linear : {{ thermo.linear }} Low frequency ansatz : {{ thermo.kind }} + | Normal Mode Wavenumbers (1-based indices) {{ sep }} {% for nu in org_nus -%} {{ "\t%04d" % loop.index }}: {{ fmt_nu(nu) }} {% endfor -%} {% if is_ts %}This should be a TS.{% endif %} Expected {{ nexpected }} real wavenumbers, got {{ nused }}. {%- if nus_below|length > 0 %} Imaginary modes present: {%- for nu_below in nus_below %} {{ "\t%04d" % loop.index }}: {{ fmt_nu(nu_below) }} {%- endfor %} {%- endif %} {{ sep }} + | Partition functions {{ sep }} {{ fmt_pf("Electronic", Q_el) }} {{ fmt_pf("Translation", Q_trans) }} {{ fmt_pf("Rotation", Q_rot) }} {{ fmt_pf("Vibration, (Bot)", Q_vib) }} {{ fmt_pf("Vibration, (V=0)", Q_vib0) }} {{ sep }} {{ fmt_pf("Tot (Bot)", Q_tot) }} {{ fmt_pf("Tot (Bot)", Q_tot0) }} + | Inner energy U = U_el + U_vib + U_rot + U_trans {{ sep }} {{ fmt("U_el", thermo.U_el) }} {{ fmt("ZPE", thermo.ZPE) }} {{ fmt("U_vib (incl. ZPE)", thermo.U_vib) }} {{ fmt("U_rot", thermo.U_rot) }} {{ fmt("U_trans", thermo.U_trans) }} {{ sep }} {{ fmt("U", thermo.U_tot) }} + | Enthalpy H = U + kB*T {{ sep }} {{ fmt("U", thermo.U_tot) }} {{ fmt("kB*T", thermo.kBT) }} {{ sep }} {{ fmt("H", thermo.H) }} + | Entropy correction T*S = T*(S_el + S_vib + S_rot + S_trans) {{ sep }} {{ fmt("T*S_el", thermo.TS_el) }} {{ fmt("T*S_vib", thermo.TS_vib) }} {{ fmt("T*S_rot", thermo.TS_rot) }} {{ fmt("T*S_trans", thermo.TS_trans) }} {{ sep }} {{ fmt("T*S", thermo.TS_tot) }} + | Gibbs free energy G = H - T*S {{ sep }} {{ fmt("H", thermo.H) }} {{ fmt("T*S", thermo.TS_tot) }} {{ sep }} {{ fmt("G", thermo.G) }} {{ fmt("dG", thermo.dG) }} """.strip() )