Hard and easy cases for DFT

7.2. Hard and easy cases for DFT#

7.2.1. Dispersion correction#

Dispersion interactions are caused by long-range electron correlation and decay asymptotically with \(E_{\text{disp}} \propto - C_6/R^6\) (London interaction). As the distance between two electron clouds is reduced this dispersion interaction becomes part of the normal correlation energy for short e-e distances. The London part of the dispersion energy is not included in standard local density functionals. The inclusion of short/medium range dispersion depends on the functional. Semi-local functionals that include the gradient of the density or the kinetic energy density tend to do better in this regard. A proper description of dispersion requires non-local correlation functionals.

Dispersion is often synomously used with Van der Waals interaction, but technically this is not correct. Dispersion is a time-independent phenomen that one can observe e.g. for non-overlapping, spherically-symmetric charge densities (e.g two argon atoms A and B at points \(\vec{r}\) and \(\vec{r}'\)). The dispersion interaction shows up as large density gradient in regions of space where the density itself is small.

This can be better described by newer exchange-correlation potentials having the following form

\[ v_c^{\text{nl}}(\vec{r}) = \int f(\vec{r}, \vec{r}') dr' \]

where a perturbation at \(\vec{r}'\) due to atom B, induces an exchange-correlation potential at \(\vec{r}\) on atom A, even at a great distance. Examples for such functionals are VV10 and vdW-DF-10.

A cheaper and simple solution is to use an empirical potential of the form \(- C_6/R^6\). A popular dispersion correction is DFT-D3 by Grimme and coworkers. Like other corrections of the DFT-D family, this is an additive correction.

\[E_{\text{B3LYP-D3}} = E_{\text{B3LYP}} + E_{\text{D3, 2-body}} + E_{\text{ATM, 3-body}}\]

The dispersion correction consists of a two-body term plus a three-body term of the Axilrod-Teller-Muto (ATM) triple-dipole variety. If the BJ approach is used, a slightly different damping function is employed to make the potential finite (but non-zero) for \(R_{AB} \rightarrow 0\). This is the preferred dispersion correction for most uses.

In Psi4 we can also add non-local density based correlation energy from the VV10 functional to e.g the b3lyp functional by using the following keyword b3lyp-nl. The VV10 functional is less empirical than the DFT-D approaches, as it contains only two fitted parameters.

7.2.2. Accuracy of DFT: Swiss army knife for computational chemistry?#

DFT although fast and available in many quantum chemistry programs is not a universal solution for any system. In principle, DFT is exact. However, in practice, as we rely on density functional approximations (DFA), the results can vary widely, depending on the DFA used. DFA are often fit using reference data from databases containing atomization energies and various other data. In this example, we will look at three reactions using geometries obtained from quantum chemical reference datasets. For reasons of efficiency, we will only perform single point calculations using a variety of functionals to learn about some points that we need to be careful about when using DFT.

The reactions of interest are:

Homolytic bond dissociation energy of CaO into atomic calcium and oxygen (taken from Minnesota 2015B database, 96.15 kcal mol\(^{-1}\))
Interaction energy of two ethane molecules (taken from GMTKN55 database, -1.34 kcal mol\(^{-1}\))
Energy difference between the staggered and eclipsed methanol conformations (taken from GMTKN55 database, -1.01 kcal mol\(^{-1}\))

Note, that one of the reactions is contained in the training data of one of the functionals that we will be using.
These exercises are based on Pierpaolo Morgante and Roberto Peverati’s 2020 Review: (https://doi.org/10.1002/qua.26332).

Exercise 5

Before running the calculations, give your best estimate as to which of the reaction energies (i.e. the energy difference between product and reactants) do you think is most “difficult” to calculate? Rank them in order of increasing difficulty and explain the reasoning behind your choices.

Difficult in this case refers to getting an accurate reaction energy with low absolute error to true energy (reference energy from the dataset). As chemical accuracy we usually use 1 kcal mol\(^-1\)

Exercise 6

Run the calculations below and check whether your assessment of the hard, medium and easy reactions is correct. Include the tables containing the calculated energies for the three problems in your report. Which functional do you think contained the reaction in its training data?

Bonus Exercise 7

How could fortitous error cancellation help in getting good results when calculating reaction energies?

Exercise 8

What is the impact of using the dispersion correction for the ethane complex? For what phenomena are such corrections relevant? Give a real world example system.

import psi4
import pandas as pd


import sys
sys.path.append("..")
from helpers import *

psi4.set_num_threads(2)
psi4.set_memory('2 GB')

  Threads set to 2 by Python driver.

  Memory set to   1.863 GiB by Python driver.

2000000000

7.2.3. BDE CaO into atomic calcium and oxygen#

cao=psi4.geometry("""
0 1
Ca 0.0000000 0.0000000 1.8210000
O 0.0000000 0.0000000 0.0000000
""")

ca = psi4.geometry("""
0 1
Ca   0.0   0.0   0.0
""")

o= psi4.geometry("""
0 3
O   0.0   0.0   0.0
""")

methods =[ 'pbe', 'b3lyp', 'mn15','mp2', 'scan'] #These are the methods that we are going to use
energies = {} #Here we will store the energies
psi4.core.set_output_file('dft_cao.log', False)

for method in methods:
    # we clone all the molecules as we want to start the geometry optimization from the same starting structure
    cao_opt = cao.clone()
    ca_opt = ca.clone()
    o_opt = o.clone()

    # call psi4  single point energy calculation once per method
    # use def2-QZVP basis set    
    psi4.set_options({'dft_radial_points':99, 'dft_spherical_points': 590, 'reference':'uks'})
    E_cao = None
    E_o = None
    E_ca = None
    # calculate the raction energy in kcal/mol 
    # you can use psi4.constants.hartree2kcalmol for the unit conversion    
    E_reaction = None 

    energies[method]= [E_cao, E_o, E_ca, E_reaction]

# insert the reference here
energies['reference']= ['','','', 96.15]

pd.DataFrame.from_dict(energies, orient='index', columns=['E_CaO', 'E_O', 'E_Ca', 'E_reaction (kcal/mol)'])

	E_CaO	E_O	E_Ca	E_reaction (kcal/mol)
pbe	None	None	None	NaN
b3lyp	None	None	None	NaN
mn15	None	None	None	NaN
mp2	None	None	None	NaN
scan	None	None	None	NaN
reference				96.15

psi4.core.clean()
psi4.core.clean_options()

7.2.4. Interaction of two Ethane molecules#

dimer = psi4.geometry("""
symmetry c1
0 1
C  0.684160244353  0.337751120631  -1.9114826827
C  -0.684160244353  -0.337751120631  -1.9114826827
H  1.48964553204  -0.401434143375  -1.9114826827
H  0.812773290289  0.969660346322  -2.794291998
H  0.812773290289  0.969660346322  -1.0286733674
H  -1.48964553204  0.401434143375  -1.9114826827
H  -0.812773290289  -0.969660346322  -2.794291998
H  -0.812773290289  -0.969660346322  -1.0286733674
C  0.684160244353  -0.337751120631  1.9114826827
C  -0.684160244353  0.337751120631  1.9114826827
H  1.48964553204  0.401434143375  1.9114826827
H  0.812773290289  -0.969660346322  2.794291998
H  0.812773290289  -0.969660346322  1.0286733674
H  -1.48964553204  -0.401434143375  1.9114826827
H  -0.812773290289  0.969660346322  2.794291998
H  -0.812773290289  0.969660346322  1.0286733674
""")
monomer = psi4.geometry("""
symmetry c1
0 1
C  0.0  0.0  0.762988243308
C  0.0  0.0  -0.762988243308
H  -0.882809099823  0.509690071426  1.1580401983
H  0.0  -1.01938014285  1.1580401983
H  0.882809099823  0.509690071426  1.1580401983
H  0.882809099823  -0.509690071426  -1.1580401983
H  -0.882809099823  -0.509690071426  -1.1580401983
H  0.0  1.01938014285  -1.1580401983
""")

drawXYZSideBySide(monomer, dimer)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

methods = ['b3lyp', 'b3lyp-d3(bj)', 'b3lyp-nl' ,'pbe', 'pbe-d3', 'mp2', 'svwn'] #These are the methods that we are going to compare
energies = {} #Here we will store the energies
psi4.core.set_output_file('dft_ethane.log', False)

for method in methods:
    # we clone all the molecules as we want to start the geometry optimization from the same starting structure
    monomer_opt = monomer.clone()
    dimer_opt = dimer.clone()

    # call psi4  single point energy calculation once per method
    # use def2-TZVP basis set
    E_monomer = None
    E_dimer = None
    # calculate the raction energy in kcal/mol 
    # you can use psi4.constants.hartree2kcalmol for the unit conversion    
    E_reaction = None
    
    energies[method] = [E_dimer, E_monomer, E_reaction]

# insert the reference here
energies['reference']= ['','', -1.34]

pd.DataFrame.from_dict(energies, orient='index', columns=['E_dimer', 'E_monomer', 'E_interaction (kcal/mol)'])

	E_dimer	E_monomer	E_interaction (kcal/mol)
b3lyp	None	None	NaN
b3lyp-d3(bj)	None	None	NaN
b3lyp-nl	None	None	NaN
pbe	None	None	NaN
pbe-d3	None	None	NaN
mp2	None	None	NaN
svwn	None	None	NaN
reference			-1.34

psi4.core.clean()
psi4.core.clean_options()

7.2.5. Eclipsed/Staggered Methanol#

methanol_eclipsed = psi4.geometry("""
0 1
C     0.1354501    0.3335445    0.0000000 
O     0.1744673   -1.0774079    0.0000000 
H     0.6439155    0.7202929    0.8794450 
H     0.6439155    0.7202929   -0.8794450 
H    -0.8812042    0.7207713    0.0000000 
H    -0.7165441   -1.4174937    0.0000000 
""")

methanol_staggered = psi4.geometry("""0 1
O     0.1640874    1.0771522    0.0000000 
C     0.1319097   -0.3310036    0.0000000 
H     1.1612146   -0.6758875    0.0000000 
H    -0.3605642   -0.7355506    0.8858582 
H    -0.3605642   -0.7355506   -0.8858582 
H    -0.7360833    1.4008402    0.0000000 
   """)

drawXYZSideBySide(methanol_eclipsed, methanol_staggered)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

methods = [ 'pbe', 'b3lyp', 'mn15','mp2', 'scan']#These are the methods that we are going to compare
energies = {} #Here we will store the energies
psi4.core.set_output_file('dft_methanol.log', False)

for method in methods: 
    # we clone all the molecules as we want to start the geometry optimization from the same starting structure
    eclipsed_opt = methanol_eclipsed.clone()
    staggered_opt = methanol_staggered.clone()

    # call psi4  single point energy calculation once per method
    # use def2-TZVP basis set
    psi4.set_options({'dft_radial_points':99, 'dft_spherical_points': 590})
    E_eclipsed = None
    E_staggered = None
    # calculate the raction energy in kcal/mol 
    # you can use psi4.constants.hartree2kcalmol for the unit conversion    
    E_reaction = None
    
    energies[method]= [E_eclipsed, E_staggered, E_reaction]

# insert the reference here
energies['reference']= ['','',-1.01]

pd.DataFrame.from_dict(energies, orient='index', columns=['E_eclipsed', 'E_staggered', 'E_reaction (kcal/mol)'])

	E_eclipsed	E_staggered	E_reaction (kcal/mol)
pbe	None	None	NaN
b3lyp	None	None	NaN
mn15	None	None	NaN
mp2	None	None	NaN
scan	None	None	NaN
reference			-1.01