patonlab / GoodVibes

Calculate quasi-harmonic free energies from Gaussian output files with temperature and other corrections
http://www.patonlab.colostate.edu
MIT License
129 stars 50 forks source link
compchem dispersion gaussian quasi-harmonic rrho temperature thermochemistry vibrational-entropies

GoodVibes

Build Status PyPI version Anaconda-Server Badge Documentation Status DOI Anaconda-Server Badge

GoodVibes is a Python program to compute thermochemical data from one or a series of electronic structure calculations. It has been used since 2015 by several groups, primarily to correct the poor description of low frequency vibrations by the rigid-rotor harmonic oscillator treatment. The current version includes thermochemistry at variable temperature/concentration, various quasi-harmonic entropy and enthalpy schemes, automated detection of frequency scaling factors, D3-dispersion corrections calculations, Boltzmann averaging, duplicate conformer detection, automated tabulation and plotting of energy profiles, and error checking. Developed by Robert Paton, Ignacio Funes-Ardoiz, and members of the Paton Research Group, Colorado State: Guilian Luchini, Juan V. Alegre-Requena, and Yanfei Guan . Integration with Travis CI testing by Jaime Rodríguez-Guerra with additions from Guilian Luchini.

All (electronic, translational, rotational and vibrational) partition functions are recomputed and will be adjusted to any temperature or concentration. These default to 298.15 Kelvin and 1 atmosphere.

The program will attempt to parse the level of theory and basis set used in the calculations and then try to apply the appropriate vibrational (zpe) scaling factor. Scaling factors are taken from the Truhlar group database.

Documentation

GoodVibes documentation can be found on our read-the-docs page.

Quasi-Harmonic Approximation

Two types of quasi-harmonic approximation are readily applied. The first is vibrational entropy: below a given cut-off value vibrational normal modes are not well described by the rigid-rotor-harmonic-oscillator (RRHO) approximation and an alternative expression is instead used to compute the associated entropy. The quasi-harmonic vibrational entropy is always less than or equal to the standard (RRHO) value obtained using Gaussian. Two literature approaches have been implemented. In the simplest approach, from Cramer and Truhlar,1 all frequencies below the cut-off are uniformly shifted up to the cut-off value before entropy calculation in the RRHO approximation. Alternatively, as proposed by Grimme,2 entropic terms for frequencies below the cut-off are obtained from the free-rotor approximation; for those above the RRHO expression is retained. A damping function is used to interpolate between these two expressions close to the cut-off frequency.

The second type of quasi-harmonic approximation available is applied to the vibrational energy used in enthalpy calculations. Similar to the entropy corrections, the enthalpy correction implements a quasi-harmonic correction to the RRHO vibrational energy computed in DFT methods. The quasi-harmonic enthalpy value as specified by Head-Gordon3 will be less than or equal to the uncorrected value using the RRHO approach, as the quasi-RRHO value of the vibrational energy used to compute the enthalpy is damped to approach a value of 0.5RT, opposed to the RRHO value of RT. Because of this, the quasi-harmonic enthalpy correction is appropriate for use in systems and reactions resulting in a loss of a rotational or translational degree of freedom.

Installation

Citing GoodVibes

Luchini, G.; Alegre-Requena, J. V.; Funes-Ardoiz, I.; Paton, R. S. GoodVibes: Automated Thermochemistry for Heterogeneous Computational Chemistry Data. F1000Research, 2020, 9, 291 DOI: 10.12688/f1000research.22758.1

Using GoodVibes

python -m goodvibes [-q] [--qs grimme/truhlar] [--qh] [-f cutoff_freq] [--fs S_cutoff_freq] [--fh H_cutoff_freq]
[--check] [-t temperature] [-c concentration] [--ti 't_initial, t_final, step'] [--ee] [--bav "global" or "conf"]
[--cosmo cosmo_filename] [--cosmoint cosmo_filename,initial_temp,final_temp] [-v frequency_scale_factor]
[--vmm mm_freq_scale_factor][--ssymm] [--spc link/filename] [--boltz] [--dup][--pes pes_yaml] [--nogconf]
[--graph graph_yaml] [--cpu] [--imag] [--invertifreq] [--freespace solvent_name] [--output output_name]
[--media solvent_name] [--xyz] [--csv] [--custom_ext file_extension] <output_file(s)>

Example 1: Grimme-type quasi-harmonic correction with a (Grimme type) cut-off of 150 cm-1

python -m goodvibes examples/methylaniline.out -f 150

   Structure                    E        ZPE             H        T.S     T.qh-S          G(T)       qh-G(T)
   *********************************************************************************************************
o  methylaniline      -326.664901   0.142118   -326.514489   0.039668   0.039465   -326.554157   -326.553954
   *********************************************************************************************************

The output shows both standard harmonic and quasi-harmonic corrected thermochemical data (in Hartree). The corrected enthalpy and entropy values are always less than or equal to the harmonic value.

Example 2: Quasi-harmonic thermochemistry with a larger basis set single point energy correction link job

python -m goodvibes examples/ethane_spc.out --spc link

   Structure                E_SPC             E        ZPE         H_SPC        T.S     T.qh-S      G(T)_SPC   qh-G(T)_SPC
   ***********************************************************************************************************************
o  ethane_spc          -79.858399    -79.830421   0.073508    -79.780448   0.027569   0.027570    -79.808017    -79.808019
   ***********************************************************************************************************************

This calculation contains a multi-step job: an optimization and frequency calculation with a small basis set followed by (--Link1--) a larger basis set single point energy. Note the use of the --spc link option. The standard harmonic and quasi-harmonic corrected thermochemical data are obtained from the small basis set partition function combined with the larger basis set single point electronic energy. In this example, GoodVibes automatically recognizes the level of theory used in the frequency calculation, B3LYP/6-31G(d), and applies the appropriate scaling factor of 0.977 (this can be suppressed to apply no scaling with -v 1.0)

Alternatively, if a single point energy calculation has been performed separately, provided both file names share a common root e.g. ethane.out and ethane_TZ.out then use of the --spc TZ option is appropriate. This will give identical results as above.

python -m goodvibes examples/ethane.out --spc TZ

   Structure                E_SPC             E        ZPE         H_SPC        T.S     T.qh-S      G(T)_SPC   qh-G(T)_SPC
   ***********************************************************************************************************************
o  ethane              -79.858399    -79.830421   0.073508    -79.780448   0.027569   0.027570    -79.808017    -79.808019
   ***********************************************************************************************************************

Example 3: Changing the temperature (from standard 298.15 K to 1000 K) and concentration (from standard state in gas phase, 1 atm, to standard state in solution, 1 mol/l)

python -m goodvibes examples/methylaniline.out –t 1000 –c 1.0

   Structure                    E        ZPE             H        T.S     T.qh-S          G(T)       qh-G(T)
   *********************************************************************************************************
o  methylaniline      -326.664901   0.142118   -326.452307   0.218212   0.216559   -326.670519   -326.668866
   *********************************************************************************************************

This correction from 1 atm to 1 mol/l is responsible for the addition 1.89 kcal/mol to the Gibbs energy of each species (at 298K). It affects the translational entropy, which is the only component of the molecular partition function to show concentration dependence. In the example above the correction is larger due to the increase in temperature.

Example 4: Analyzing the Gibbs energy across an interval of temperatures 300-1000 K with a stepsize of 100 K, applying a (Truhlar type) cut-off of 100 cm-1

python -m goodvibes examples/methylaniline.out --ti '300,1000,100' --qs truhlar -f 120

   Structure               Temp/K                        H        T.S     T.qh-S          G(T)       qh-G(T)
   ******************************************************************************************************
o  methylaniline            300.0              -326.514399   0.040005   0.039842   -326.554404   -326.554241
o  methylaniline            400.0              -326.508735   0.059816   0.059596   -326.568551   -326.568331
o  methylaniline            500.0              -326.501670   0.082625   0.082349   -326.584296   -326.584020
o  methylaniline            600.0              -326.493429   0.108148   0.107816   -326.601577   -326.601245
o  methylaniline            700.0              -326.484222   0.136095   0.135707   -326.620317   -326.619930
o  methylaniline            800.0              -326.474218   0.166216   0.165772   -326.640434   -326.639990
o  methylaniline            900.0              -326.463545   0.198300   0.197800   -326.661845   -326.661346
o  methylaniline           1000.0              -326.452307   0.232169   0.231614   -326.684476   -326.683921
   ******************************************************************************************************

Note that the energy and ZPE are not printed in this instance since they are temperature-independent. The Truhlar-type quasi-harmonic correction sets all frequencies below than 120 cm-1 to a value of 100. Constant pressure is assumed, so that the concentration is recomputed at each temperature.

Example 5: Analyzing the Gibbs Energy using scaled vibrational frequencies

python -m goodvibes examples/methylaniline.out -v 0.95

   Structure                    E        ZPE             H        T.S     T.qh-S          G(T)       qh-G(T)
   *********************************************************************************************************
o  methylaniline      -326.664901   0.135012   -326.521265   0.040238   0.040091   -326.561503   -326.561356
   *********************************************************************************************************

The frequencies are scaled by a factor of 0.95 before they are used in the computation of the vibrational energies (including ZPE) and entropies.

Example 6: Writing Cartesian coordinates

python -m goodvibes examples/HCN*.out --xyz

Optimized cartesian-coordinates found in files HCN_singlet.out and HCN_triplet.out are written to Goodvibes_output.xyz

Example 7: Analyzing multiple files at once

python -m goodvibes examples/*.out --cpu

   Structure                    E        ZPE             H        T.S     T.qh-S          G(T)       qh-G(T)
   *********************************************************************************************************
o  Al_298K            -242.328708   0.000000   -242.326347   0.017670   0.017670   -242.344018   -242.344018
o  Al_400K            -242.328708   0.000000   -242.326347   0.017670   0.017670   -242.344018   -242.344018
o  H2O                 -76.368128   0.020772    -76.343577   0.021458   0.021458    -76.365035    -76.365035
o  HCN_singlet         -93.358851   0.015978    -93.339373   0.022896   0.022896    -93.362269    -93.362269
o  HCN_triplet         -93.153787   0.012567    -93.137780   0.024070   0.024070    -93.161850    -93.161850
o  allene             -116.569605   0.053913   -116.510916   0.027618   0.027621   -116.538534   -116.538537
o  benzene            -232.227201   0.101377   -232.120521   0.032742   0.032745   -232.153263   -232.153265
o  ethane              -79.830421   0.075238    -79.750770   0.027523   0.027525    -79.778293    -79.778295
o  isobutane          -158.458811   0.132380   -158.319804   0.034241   0.034252   -158.354046   -158.354056
o  methylaniline      -326.664901   0.142118   -326.514489   0.039668   0.039535   -326.554157   -326.554024
o  neopentane         -197.772980   0.160311   -197.604824   0.036952   0.036966   -197.641776   -197.641791
   *********************************************************************************************************
TOTAL CPU      0 days  2 hrs 29 mins 28 secs

The program will detect several different levels of theory and give a warning that any vibrational scaling factor other than 1 would be inappropriate in this case. Wildcard characters (*) can be used to represent any character or string of characters.

Example 8: Entropic Symmetry Correction

python -m goodvibes examples/allene.out examples/benzene.out examples/ethane.out examples/isobutane.out examples/neopentane.out --ssymm

   Structure                    E        ZPE             H        T.S     T.qh-S          G(T)       qh-G(T)  Point Group
   **********************************************************************************************************************
o  allene             -116.569605   0.053913   -116.510916   0.026309   0.026312   -116.537225   -116.537228          D2d
o  benzene            -232.227201   0.101377   -232.120521   0.030396   0.030399   -232.150917   -232.150919          D6h
o  ethane              -79.830421   0.075238    -79.750770   0.025831   0.025833    -79.776601    -79.776603          D3d
o  isobutane          -158.458811   0.132380   -158.319804   0.033204   0.033214   -158.353008   -158.353019          C3v
o  neopentane         -197.772980   0.160311   -197.604824   0.034606   0.034620   -197.639430   -197.639444           Td
   *********************************************************************************************************************************************

GoodVibes will apply a symmetry correction described above to the entropy term of each molecule after determining the symmetry number. It is always a good idea to double-check that the point group GoodVibes returns is the correct group.

Example 9: Potential Energy Surface (PES) Comparison with Accessible Conformer Correction

python -m goodvibes examples/gconf_ee_boltz/*.log --pes examples/gconf_ee_boltz/gconf_aminox_cat.yaml

   Structure                       E        ZPE             H        T.S     T.qh-S          G(T)       qh-G(T)
   ************************************************************************************************************
o  Aminoxylation_TS1_R   -879.405138   0.295352   -879.091374   0.063746   0.061481   -879.155120   -879.152855
o  Aminoxylation_TS2_S   -879.404445   0.295301   -879.090562   0.064366   0.061891   -879.154928   -879.152453
o  aminox_cat_conf212_S  -517.875165   0.200338   -517.662195   0.051817   0.049814   -517.714012   -517.712009
o  aminox_cat_conf280_R  -517.877308   0.200869   -517.664171   0.049996   0.048777   -517.714167   -517.712948
o  aminox_cat_conf65_S   -517.877161   0.200789   -517.664159   0.049790   0.048656   -517.713949   -517.712815
o  aminox_subs_conf713   -361.535757   0.095336   -361.433167   0.037824   0.037696   -361.470991   -361.470863
   ************************************************************************************************************

   Gconf correction requested to be applied to below relative values using quasi-harmonic Boltzmann factors

   RXN: Reaction (kcal/mol)       DE       DZPE            DH       T.DS    T.qh-DS         DG(T)      qh-DG(T)
   ************************************************************************************************************
o  Cat+Subs                     0.00       0.00          0.00       0.00       0.00          0.00          0.00
o  TS                           4.72      -0.46          3.53     -15.85     -16.37         19.39         19.90
   ************************************************************************************************************

A .yaml file is given to the --pes argument which specifies the reaction: Catalyst + Substrate -> TS. Because multiple conformers for the catalysts and transition states have been provided, GoodVibes will calculate a correction to the free energy based on the number of accessible conformations based on the Boltzmann-weighted energies of the conformers. To turn this correction off, --nogconf should be specified. An example .yaml file is shown below to show how these files should be formatted.

Example 10: Stereoselectivity and Boltzmann populations

python -m goodvibes examples/gconf_ee_boltz/Aminoxylation_TS1_R.log examples/gconf_ee_boltz/Aminoxylation_TS2_S.log --boltz --ee "*_R*:*_S*"

   Structure                       E        ZPE             H        T.S     T.qh-S          G(T)       qh-G(T)  Boltz
   *******************************************************************************************************************
o  Aminoxylation_TS1_R   -879.405138   0.295352   -879.091374   0.063746   0.061481   -879.155120   -879.152855  0.605
o  Aminoxylation_TS2_S   -879.404445   0.295301   -879.090562   0.064366   0.061891   -879.154928   -879.152453  0.395
   *******************************************************************************************************************

   Selectivity            Excess (%)     Ratio (%)         Ratio     Major Iso           ddG
   *****************************************************************************************
o                              20.98         60:40         1.5:1             R          0.25
   *****************************************************************************************

The --boltz option will provide Boltzmann probabilities to the right of energy results under the boltz tab. With the --ee option, %ee, er and a reduced ratio are shown along with the dominant isomer and a calculated transition state energy value, ddG or ΔG‡.

Checks

A computational workflow can become less effective without consistency throughout the process. By using the --check option, GoodVibes will enforce a number of pass/fail checks on the input files given to make sure uniform options were used. Checks employed are:

Gaussian Output Checks
Single Point Calculation Checks

Symmetry

GoodVibes is able to detect a probable symmetry point group for each species and apply a symmetry correction to the entropy (Ssym) by finding a molecule's internal symmetry number using atom connectivity, and external symmetry with the help of the external open source C program, "Brute Force Symmetry Analyzer" developed by S. Patchkovskii. These numbers are combined to give a symmetry number, n, and Ssym is then defined as -Rln(n), which is applied to the GoodVibes calculated entropy. Note: this option may not function properly on some versions of Windows.

File Naming Conventions

Some options (--pes, --graph, --spc, --ee, --media) require the calculation output files to be named in a certain way for GoodVibes to recognize them and perform extra calculations properly.

.yaml File Formatting

When using the --pes or --graph options in GoodVibes, a .yaml file must be provided to the program to specify qualities like reaction pathways, provided conformers, and other formatting options. The same .yaml may be used for both --pes and --graphing options. An example .yaml file from an external Zenodo repository is shown below:

--- # PES
    Ph: [Ph-Int1 + EtOH, Ph-TS1 + EtOH, Ph-Int2 + EtOH, Ph-TS2 + EtOH, Ph-Int3 + EtOH]
    Py: [Py-Int1 + EtOH, Py-TS1 + EtOH, Py-Int2 + EtOH, Py-TS2 + EtOH, Py-Int3 + EtOH]

--- # SPECIES
    Ph-Int1     : Int_I_Ph*
    Ph-TS1      : TS_1_Ph*
    Ph-Int2     : Int_II_Ph*
    Ph-TS2      : TS_II_Ph*
    Ph-Int3     : Int_III_Ph*
    Py-Int1     : Int_I_Py*
    Py-TS1      : TS_1_Py*
    Py-Int2     : Int_II_Py*
    Py-TS2      : TS_II_Py*
    Py-Int3     : Int_III_Py*
    EtOH        : ethanol*

--- # FORMAT
    dec :  2 
    legend : False
    color : black,#26a6a4
    pointlabel : False 
    gridlines: True
    show_conformers: True
    show_gconf: False
    dpi : 400
    title: Potential Energy Surface

options in the # FORMAT block are optional, but allow for stylistic choices to be employed, especially when graphing. All current options that can be specified for either --pes or --graph options are:

    dec : 1 or 2 (decimal points in output)
    legend : True or False (puts legend on graph)
    ylim : y_min,y_max (y axis limits on graph)
    color : Color (color of line for a reaction pathway, multiple pathways can have different colors i.e. color1,color2,color3 etc., this follows rules for matplotlib standard colors)
    pointlabel : True or False (labels relative energy on graph at point)
    xlabel : True or False (displays structure labels at pathway points on x axis)
    title : Title (title displayed on graph)
    gridlines: True or False (displays gridlines on graph)
    dpi : number (specify dpi (dots per inch) of an image, will automatically save output image at specified dpi)
    show_conformers : True or False (displays a point for each conformer of a certain compound at its relative energy on graph)
    show_gconf : True or False (displays the effect of multiple accessible conformers correction if applied)

Tips and Troubleshooting

Papers from other research groups citing GoodVibes

References for the underlying theory

  1. Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2011, 115, 14556-14562 DOI: 10.1021/jp205508z
  2. Grimme, S. Chem. Eur. J. 2012, 18, 9955–9964 DOI: 10.1021/jp509921r
  3. Li, Y.; Gomes, J.; Sharada, S. M.; Bell, A. T.; Head-Gordon, M. J. Phys. Chem. C 2015, 119, 1840-1850 DOI: 10.1002/chem.201200497
  4. Alecu, I. M.; Zheng, J.; Zhao, Y.; Truhlar, D. G.; J. Chem. Theory Comput. 2010, 6, 2872-2887 DOI: 10.1021/ct100326h
  5. Mammen, M.; Shakhnovich, E. I.; Deutch, J. M.; Whitesides, G. M. J. Org. Chem. 1998, 63, 3821-3830 DOI: 10.1021/jo970944f

License:

GoodVibes is freely available under an MIT License