Thermodynamic Analysis


I have been running thermodynamic analysis for hydrogen bonded complexes. I would like to know how to retrieve detailed results of thermodynamic analysis as a dictionary or using get_variable() when running the freq() and thermo() functions. In this case, ‘detailed’ means I would like to know how to access the electronic, translational, rotational, and vibration components of the calculated thermodynamic variables directly in my code.

I started out in the Harmonic Vibrational Analysis and Visualization of Normal Modes documentation and set up the following input file. Frequency returns a tuple containing total energy and a wavefunction object. Thermodynamic analysis is also run whenever frequency is called, but the freq() function doesn’t have an option to include a dictionary of thermodynamic results in the returned tuple. It is fairly straight forward to access the selected thermodynamic variables using get_variable().

molecule water {
O           -0.055763247595     0.035150988973     0.000000000000
H            0.846077608192     0.361016510624     0.000000000000
H            0.038925597540    -0.918888214072     0.000000000000

set basis aug-cc-pVTZ
set T 273.15
set P 101325

water_e, water_wfn = freq('mp2',molecule=water,return_wfn=True)

print(get_variable('enthalpy correction')) 
print(get_variable('thermal energy correction'))
print(get_variable('thermal energy')) 
print(get_variable('gibbs free energy correction'))
print(get_variable('gibbs free energy'))

The output file contains a detailed thermodynamic analysis table including details like the electronic, translational, rotational, and vibrational components. I included a portion of this summary table below for reference.

  Thermal Energy, E (includes ZPE)
    Electronic E          0.000 [kcal/mol]       0.000 [kJ/mol]      0.00000000 [Eh]
    Translational E       1.192 [kcal/mol]       4.989 [kJ/mol]      0.00190009 [Eh]
    Rotational E          1.192 [kcal/mol]       4.989 [kJ/mol]      0.00190009 [Eh]
    Vibrational E        13.510 [kcal/mol]      56.524 [kJ/mol]      0.02152895 [Eh]
  Correction E           15.894 [kcal/mol]      66.502 [kJ/mol]      0.02532913 [Eh]
  Total E, Electronic energy at  400.00 [K]                        -76.31876251 [Eh]

  Enthalpy, H_trans = E_trans + k_B * T
    Electronic H          0.000 [kcal/mol]       0.000 [kJ/mol]      0.00000000 [Eh]
    Translational H       1.987 [kcal/mol]       8.314 [kJ/mol]      0.00316681 [Eh]
    Rotational H          1.192 [kcal/mol]       4.989 [kJ/mol]      0.00190009 [Eh]
    Vibrational H        13.510 [kcal/mol]      56.524 [kJ/mol]      0.02152895 [Eh]
  Correction H           16.689 [kcal/mol]      69.827 [kJ/mol]      0.02659585 [Eh]
  Total H, Enthalpy at  400.00 [K]                                 -76.31749579 [Eh]

  Gibbs free energy, G = H - T * S
    Electronic G          0.000 [kcal/mol]       0.000 [kJ/mol]      0.00000000 [Eh]
    Translational G     -12.440 [kcal/mol]     -52.049 [kJ/mol]     -0.01982455 [Eh]
    Rotational G         -3.892 [kcal/mol]     -16.286 [kJ/mol]     -0.00620294 [Eh]
    Vibrational G        13.494 [kcal/mol]      56.458 [kJ/mol]      0.02150372 [Eh]
  Correction G           -2.839 [kcal/mol]     -11.877 [kJ/mol]     -0.00452377 [Eh]
  Total G, Free enthalpy at  400.00 [K]                            -76.34861541 [Eh]

I would really like to access the electronic, translational, rotational, and vibration components of the calculated thermodynamic variables in my code. Unfortunately, the same approach I took before didn’t work, some failed examples include get_variable(‘vibrational enthalpy’) and get_variable(‘vibrational entropy’).

The Vibrational and Thermodynamic Analysis documentation indicates that thermo() should return a dictionary that contains detailed thermodynamic variables.

But additionally, every valid combination of {S, Cv, Cp, ZPE, E, H, G} with {elec, trans, rot, vib, corr, tot} (e.g., vibrational entropy, S_vib, and enthalpy correction, H_corr) is returned by dictionary from the thermo function. See python/vibanalysis (near the end) for an example.

I called thermo() using the syntax outlined in the Harmonic Vibrational Analysis documentation (see link above) and the function seemed to work just fine. The output file contains the same detailed thermodynamic analysis table as freq(). However, the return was a float and not a tuple or dictionary.

therm = thermo(water_wfn,water_wfn.frequencies())

The Vibrational and Thermodynamic Analysis documentation links to the vibanalysis python example and to the freq-isotope2 example for additional information.

Unfortunately, the syntax for using thermo() shown in vibanalysis looks like it is for use in Python and is quite different from what I have seen when using psi4 as an executable. I am not really sure what the required inputs are to thermo() as it was used in vibanalysis (maybe there is a documentation page??). I tried adapting it in my input file and didn’t get very far (error messages like " AttributeError: module ‘psi4.driver.qcdb’ has no attribute ‘vib’ " were common).

vibanalysis snippet:

psi4.optimize(‘hf/3-21g’, molecule=ch4_thermo_mol)
E, thermo_wfn = psi4.freq(‘hf/3-21g’, return_wfn=True, molecule=ch4_thermo_mol, dertype=1)
print(psi4.get_variable(‘CURRENT ENERGY’))

thermo_vibinfo = psi4.driver.vibanal_wfn(thermo_wfn)

therminfo, text = qcdb.vib.thermo(thermo_vibinfo, T=298.15, P=101325.0, multiplicity=1, molecular_mass=16.0313, E0=-39.97687755, sigma=12, rot_const=rot_const)

The syntax for using thermo() shown in freq-isotope2 looks similar to what I am used to when using psi4 as an executable. I tried adapting this and also didn’t get very far since it looks like vibanal_wfn() isn’t actually a callable function (error was “NameError: name ‘vibanal_wfn’ is not defined”).

Right now, thermo() triggers the thermodynamic analysis which is then written to the output file, but thermo() doesn’t seem to return any useful objects. Is it possible to modify the thermo() function so that it returns a dictionary containing all the details of the thermodynamic analysis?


Thanks for the detailed post. I agree the interface needs some work for general (not python-happy) use. In the meantime, the below changes to freq-isotope2 show how to grab the info.

>>> git diff tests/freq-isotope2/input.dat 
diff --git a/tests/freq-isotope2/input.dat b/tests/freq-isotope2/input.dat
index b43ddae..2df606b 100644
--- a/tests/freq-isotope2/input.dat
+++ b/tests/freq-isotope2/input.dat
@@ -89,6 +89,7 @@ compare_values(46.450, entropy, 2, 'HDO S')  # molpro default #TEST
+set hessian_write on
 vibanal_wfn(wfn, molecule=dto)
 compare_values(-74.96590119, get_variable('current energy'), 6, 'DTO E0')  #TEST
 compare_values(0.016317, get_variable('ZPVE'), 4, 'DTO ZPVE')  #TEST
@@ -97,3 +98,19 @@ compare_values(0.020098, get_variable('enthalpy correction'), 4, 'DTO dH')  #TES
 entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_g
 compare_values(49.603, entropy, 2, 'DTO S')  # molpro sym=cs #TEST
+# A. print all the vib-related info from original `e, wfn = freq()` call
+for k, v in wfn.frequency_analysis.items():
+    print(v)
+# B. `set hessian_write on` above wrote all the vib- & thermo-related info from last analysis.
+#    now, print all the vib- & thermo-related info from that file
+with open(core.get_writer_file_prefix( + ".vibrec") as handle:
+    import json
+    jsondata = json.load(handle)
+for k, v in jsondata.items():
+    print(k, v)
+compare_values(0.016317, jsondata['ZPE_corr']['data'], 4, 'DTO ZPVE from file')  #TEST
+compare_values(0.020098, jsondata['H_corr']['data'], 4, 'DTO dH from file')  #TEST
+compare_values(298.15, jsondata['T']['data'], 2, 'temp') #TEST