Hello,
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('enthalpy')) print(get_variable('thermal energy correction')) print(get_variable('thermal energy')) print(get_variable('gibbs free energy correction')) print(get_variable('gibbs free energy')) print(get_variable('zpve'))
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.core.clean()
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)
rot_const=np.asarray(ch4_thermo_mol.rotational_constants())
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!
John