7. Performing a typical computational study

There are of course many different types of computational studies but most of them include the following steps:

geometry optimization of the molecular species (perhaps saddle-point optimizations as well) and frequency calculations, single-point energy energy evaluations and possibly a spectroscopic property calculation.

It is common in computational studies to be working with 2 different level of theories: a low-level theory and a high-level theory.

The low-level theory is used for geometry optimization and frequencies (frequencies are used to compute corrections like zero-point energy corrections, enthalpy, entropy and free energies) and the high-level theory is used to calculate either accurate energies or accurate spectroscopic properties.

The low-level theory is a method/basis set combination that gives you reasonably accurate results for geometries and frequencies in a reasonably short time. Because the accuracy of geometries and frequencies depend less on basis sets and different methods than energy evaluations yet optimizations and frequencies are expensive, this is a rather good approximation to save time. For example using a fast GGA-DFT method like BP86 and PBE and using the RI approximation (automatic in ORCA) and a small basis set (maybe def2-SVP for organic elements and def2-TZVP for metals) would be a good low-level theory for many molecules allowing you to quickly optimize the geometry of your system.

! BP86 def2-SVP def2-SVP/J TIGHTSCF opt freq

Note: Frequency calculations always have to be performed at the same level of theory as the geometry optimization. You should never do a frequency calculation on a geometry if you are not sure whether the geometry was optimized at the precise same level of theory. Harmonic vibrational frequencies don't make sense if they are not evaluated at a point where the gradient is close to zero (the minimum).

To make sure of this, it is best to always include both Opt and Freq keywords (even though the geometry might already have been optimized) when you are performing frequency calculations. Note that often it makes sense to postpone doing Frequency jobs until later in the project since they are expensive and you don't always need the information (vibrational frequencies and thermochemical corrections).

The high-level theory is where you increase the basis set as much as possible to get rid of basis set artifacts (e.g. def2-TZVPP or def2-QZVPP even) and perhaps choose a more accurate but more expensive method (e.g. a hybrid-DFT method or even a post-HF method).

You use the high-level theory for a single-point energy evaluation (and perhaps spectroscopic/molecular property calculation) on the previously optimized geometry. Because you are only calculating a single energy (and not multiple energy evaluations like in optimizations) you can usually afford to crank up the accuracy here without the calculation becoming too expensive.

! PBE0 def2-TZVPP def2-TZVPP/J TIGHTSCF

Note: You have to use the high-level theory to compute the energies of all your species, not just some of them. Energies change dramatically when you change the basis set or method. You will only get sensible results for your relative energies or reaction energies if you compute everything exactly the same way. You should make sure you are optimizing all species in the same way (i.e. use the low-level theory for all species).

Step by Step Example:

You are computing the reaction: N₂ + 3H₂ -> 2NH₃

You are interested in computing the reaction energy for this reaction with 3 different thermodynamic definitions:

Uncorrected Reaction energy at 0K (also called Potential energy of reaction or Electronic energy of reaction)

Zero-point corrected Reaction energy at 0K

Free energy of Reaction at 298K

1. Low-level theory calculations

You perform geometry optimizations and frequency calculations for all 3 species: N₂, H₂ and NH₃.

Your inputfiles look like (starting coordinates created manually or via Chemcraft):

N₂:

! BP86 def2-SVP def2-SVP/J TIGHTSCF opt freq

*xyz 0 1

N 0 0 0

N 0 0 1

*

_H2:

! BP86 def2-SVP def2-SVP/J TIGHTSCF opt freq

*xyz 0 1

H 0 0 0

H 0 0 1

*

NH₃:

! BP86 def2-SVP def2-SVP/J TIGHTSCF opt freq

*xyz 0 1

N 0.976859589 -0.526053205 3.625355820

H 0.146060739 0.300827253 3.843753573

H 0.651670957 -1.631485881 3.318880766

H 1.302048221 -0.289907124 4.747916905

*

From these calculations you get the following optimized geometries (found towards the bottom of the outputfile or in the jobname.xyz file):

N₂:

N 0.00000000000000 0.00000000000000 -0.05602736815196

N 0.00000000000000 0.00000000000000 1.05602736815196

_H2:

H 0.00000000000000 0.00000000000000 0.11667543233613

H 0.00000000000000 0.00000000000000 0.88332456766387

NH₃:

N 0.98681409702893 -0.43101948464824 3.67509788187030

H 0.12850538257861 0.13891012868357 3.72188182958331

H 0.64393522882764 -1.38533836525978 3.48772904600280

H 1.31738479756482 -0.46917123577554 4.65119830654358

Since you performed a frequency job for each species in these calculations you also have access to the Thermochemistry output that includes corrections that you will use later, such as zero-point energy correction and free energy correction etc.

See Tutorial 6. Energies & Thermochemistry for more details regarding this.

The relevant information here is:

h2.out:Zero point energy ... 0.00976212 Eh 6.13 kcal/mol

n2.out:Zero point energy ... 0.00543247 Eh 3.41 kcal/mol

nh3.out:Zero point energy ... 0.03320734 Eh 20.84 kcal/mol

h2.out:G-E(el) ... -0.00178522 Eh -1.12 kcal/mol

n2.out:G-E(el) ... -0.01302897 Eh -8.18 kcal/mol

nh3.out:G-E(el) ... 0.01513015 Eh 9.49 kcal/mol

2. High-level theory calculations

You now want to perform more accurate energies using the already optimized geometries from the low-level theory. Hence you perform single-point energy evaluations. Here we are going to switch to a larger basis set def2-TZVPP and use the hybrid-DFT method PBE0 as an example. We are also going to use the RIJCOSX approximation which is recommended for hybrid-DFT calculations.

Note: that we could also have employed a post-HF method here like MP2 or coupled cluster theory. If we wanted to compute molecular properties or spectroscopic properties like a UV-VIS spectrum then we would also include more keywords.

Important:

- Remember to use the coordinates of the optimized geometries (from the low-level theory step above) in the inputfiles (as seen below).

- Remember to not have any Opt or Freq keyword in the inputfiles. You just want to do a single energy evaluation. The point is to NOT do expensive geometry optimizations/frequency calculations using the expensive high-level theory.

N₂:

! PBE0 RIJCOSX def2-TZVPP def2-TZVPP/J TIGHTSCF

*xyz 0 1

N 0.00000000000000 0.00000000000000 -0.05602736815196

N 0.00000000000000 0.00000000000000 1.05602736815196

*

_H2:

! PBE0 RIJCOSX def2-TZVPP def2-TZVPP/J TIGHTSCF

*xyz 0 1

H 0.00000000000000 0.00000000000000 0.11667543233613

H 0.00000000000000 0.00000000000000 0.88332456766387

*

NH₃:

! PBE0 RIJCOSX def2-TZVPP def2-TZVPP/J TIGHTSCF

*xyz 0 1

N 0.98681409702893 -0.43101948464824 3.67509788187030

H 0.12850538257861 0.13891012868357 3.72188182958331

H 0.64393522882764 -1.38533836525978 3.48772904600280

H 1.31738479756482 -0.46917123577554 4.65119830654358

*

Note of course that for this simple example of 2-4 atom molecules it is not terribly expensive to simply do all the calculations (opt and freq and energies) at the high-level of theory. But most likely you are interested in computing larger molecules than these and then it will save you a great deal of time to employ a low-level & high-level protocol like this.

From the single-point energy calculations you get the following energies (easiest to grep it like this: grep FINAL *.out where * means "anything")

h2.out:FINAL SINGLE POINT ENERGY -1.168316938077

n2.out:FINAL SINGLE POINT ENERGY -109.444967092759

nh3.out:FINAL SINGLE POINT ENERGY -56.514955364454

These energies are called total energies, potential energies or electronic energies and are the major contribution to a reaction energy.

These energies are now your high-level theory energies evaluated on low-level theory geometries.

This gives us the ∆E_el reaction energy (reaction energy at 0 Kelvin but neglects zero-point motion):

∆Eel = 2*E_el(NH₃) - Eel(N₂) - 3*Eel(H2)

= 2*-56.514955364454 - 1*-109.444967092759 - 3*-1.168316938077 = -0.079992822 hartree = -50.196 kcal/mol

But what if we want a zero-point corrected reaction energy? Then we simply include the previously calculated zero-point energy corrections (from the low-level theory).

We can do it in 2 ways.

a) Compute first the zero-point correction to the reaction energy:

∆ZPEcorr = 2*ZPE(NH₃) - ZPE(N₂) - 3*ZPE(H2)

= 2*0.03320734 - 1*0.00543247 - 3*0.00976212= 0.03169585 hartree = +19.889 kcal/mol

Then add the ∆ZPEcorr to the previously calculated ∆Eel to get our zero-point corrected reaction energy:

∆E⁰ = ∆Eel + ∆ZPEcorr = -50.196 kcal/mol + +19.889 kcal/mol

= -30.307 kcal/mol

b) The other way is to add the zero-point correction first to each electronic energy of each species:

E⁰(N₂) = E_el(N2) + ZPE(N2) = -109.444967092759 + 0.00543247 = -109.4395346 hartree

E⁰(H₂) = E_el(H₂) + ZPE(H₂) = -1.168316938 + 0.00976212 = -1.158554818 hartree

E⁰(NH₃) = E_el(NH₃) + ZPE(NH₃) = -56.51495536 + 0.03320734 = -56.48174802 hartree

Then: ∆E⁰ = 2*E⁰(NH₃) - E⁰(N₂) -3*E⁰(H₂) =2*-56.48174802 - 3*-1.158554818 = -0.048296972 hartree

= -30.307 kcal/mol

You can do it analogously to get the free energy of the reaction at 298 K.

Do your results not make any sense when you compare to experiment?

- It is possible that your protocol is not accurate enough. Most likely, the energy evaluation from the high-level theory has some error. Perhaps you need to use a larger basis set or a more accurate method (perhaps post-HF). This is usually the first thing to check.

- It is also possible but less likely that :

- your geometry is not accurate enough (perhaps a larger basis set is needed there or perhaps a hybrid-DFT method instead of a GGA-DFT method).

- perhaps even the zero-point energy correction from the low-level theory is not good enough.

It is also possible that you have done some arithmetic error somewhere or that the SCF did not converge.

Or perhaps experiment is wrong or the reaction is not as simple as you thought it was.