Electron Correlation and Correlation Energy
Copyright (c) RDMCHEM LLC 2019
Overview
Bond Dissociation of Hydrogen H2
Molecular Orbitals of Nitrogen (N2)
Trimethylenemethane (TMM)
References
For an N-electron system, the Hartree-Fock approximation reduces the N-body Schrodinger equation down to N one-body Schrodinger-like equations, each describing the motion of a single electron in the mean field due to all other electrons. Furthermore, it approximates the N-body wavefunction as a single Slater determinant. These approximations treat the electrons as being uncorrelated. In reality, electrons are correlated through instantaneous repulsions and the possibility of quantum entanglement. The correlation energy is defined as the difference between the "exact" energy for a given basis set and the Hartree-Fock energy:
ϵcorr=ϵexact−ϵHF
Based on the variational principle, ϵexact is lower in energy than ϵHF , so ϵcorr is always negative. In this activity, you will consider three pedagogical examples where Hartree-Fock theory fails to give even qualitative agreement with experiment: 1) bond dissociation of H2, 2) molecular orbital energies of N2, and 3) molecular orbital energies of trimethylenemethane (TMM).
For each example, you will first calculate energies using Hartree-Fock and then explore some electronic structure methods capable of treating electron correlation in varying degrees of approximation: Moller-Plesset perturbation theory (MP2), configuration interaction, the parametric two-electron reduced density matrix (P2RDM) method, and the variational two-electron reduced density matrix (V2RDM) method.
Moller-Plesset perturbation theory (MP2):
Coupled Cluster:
Configuration Interaction:
Variational 2RDM:
Density Functional Theory:
Consider what happens as the H2 bond is stretched. The energy should decrease until the distance between the two H atoms is equal to the equilibrium bond distance. Beyond this distance, the energy should increase until the bond effectively breaks, giving rise to 2 H atoms at an 'infinite' distance.
To visualize this process, calculate the potential energy surface as a function of bond distance, r, for H2 using Hartree-Fock energy. Begin by choosing an atomic orbital basis, minimum and maximum bond distances, and the spacing between each point on the potential surface. [Note: It is not necessary to go to 'infinity'.
restart:withQuantumChemistry:withplots:Digits≔15:
AObasis≔6-31g:
rmin≔0.4:
rmax≔5.0:
dr≔0.1:
The PES can be constructed by creating a sequence of bond distances, rdata, and corresponding molecules and using the map function to calculate the energy of each molecule in the sequence:
rdata ≔ seqr, r=rmin..rmax,dr:
molecules ≔ seqH,0,0,0,H,0,0,i, i=rmin..rmax,dr:
edata_HF ≔ mapEnergy,molecules, method=HartreeFock,basis=AObasis:
pes_HF ≔ splinerdata,edata_HF,R,1:
plot_HF ≔ plot pes_HF, R=rmin..rmax, axes=BOXED, font=ROMAN,16, labelfont=ROMAN,16, labels='R','E', color=blue;
Now, let's calculate the Hartree-Fock energy of an individual H atom. The energy at dissociation should be 2× energy of a single H atom.
atom≔H,0,0,0;
atom≔H,0,0,0
dissociation≔2⋅ Energyatom,HartreeFock,basis=AObasis;
dissociation≔−0.99646582
Based on the energy of two separated H atoms, does the PES curve produce qualitative results?
Answer
NO! At 5 Angstroms, the energy is almost 0.25 Hartrees too large, or 656 kJ/mol too large!!!
Now, calculate the potential energy curve for each of the MP2, CC, and CI methods:
Moller-Plesset:
edata_MP2≔ mapEnergy,molecules, method=MP2,basis=AObasis:
pes_MP2 ≔ splinerdata,edata_MP2,R,1:
plot_MP2 ≔ plot pes_MP2, R=rmin..rmax, axes=BOXED, font=ROMAN,16, labelfont=ROMAN,16, labels='R','E', color=red:
plotsdisplayplot_HF,plot_MP2
Parametric 2-RDM:
edata_P2RDM≔ mapEnergy,molecules, method=Parametric2RDM,basis=AObasis:
pes_P2RDM ≔ splinerdata,edata_P2RDM,R,1:
plot_P2RDM ≔ plot pes_P2RDM, R=rmin..rmax, axes=BOXED, font=ROMAN,16, labelfont=ROMAN,16, labels='R','E', color=green:
plotsdisplayplot_HF,plot_P2RDM
Configuration Interaction
edata_CI≔ mapEnergy,molecules, method=FullCI,basis=AObasis:
pes_CI ≔ splinerdata,edata_CI,R,1:
plot_CI ≔ plot pes_CI, R=rmin..rmax, axes=BOXED, font=ROMAN,16, labelfont=ROMAN,16, labels='R','E', color=black:
plotsdisplayplot_HF,plot_CI
To compare, combine all 4 potential energy curves in a single plot:
plotsdisplayplot_HF,plot_MP2,plot_P2RDM,plot_CI;
Notice that all three correlation methods are below the HF curve, indicating a negative correlation energy, as expected. We find that while MP2 does a good job calculating the correlation energy at small bond distances, it also fails to describe dissociation. The parametric 2-RDM method, on the other hand, is essentially equal to the configuration interaction result, the exact numerical result for the AObasis of interest. Both the P2RDM and CI methods do qualitatively capture the key physics of dissociation at E = −1, as expected!!
Why did Hartree-Fock do so poorly?? The answer lies in the fact that HF 'restricts' electrons to be doubly occupied orbitals. So at dissociation, Hartree-Fock gives the energy of H- + H+, not H + H! An 'unrestricted' version of HF does a better job at describing dissociation but with poor accuracy.
The molecular orbital diagram in Figure 1 shows the correct order and symmetry of the molecular orbitals for dinitrogen (N2). The HOMO has σ-symmetry and is non-degenerate. (Note that the σ and σ* orbitals formed from the overlap of N1s orbitals are not shown.)
Figure 1: Molecular orbital diagram for dinitrogen. (TC Reuter, CC BY-SA 4.0)
Use the Maple input below to calculate the molecular orbital energies of N2 using Hartree-Fock and the 6-31G basis set.
withQuantumChemistry:Digits≔15:
molec≔N,0,0,0,N,0,0,1.09;
molec≔N,0,0,0,N,0,0,1.09000000
AObasis≔sto-3g:
data≔HartreeFockmolec;
For N2 with 14 electrons all paired, the HOMO is the 7th MO, as indicated by the occupation numbers of the molecular orbitals:
datamo_occ1..9;
List the molecular orbital energies and plot the HOMO, HOMO-1, and HOMO-2 molecular orbitals using the Maple input below:
datamo_energy1..9;
DensityPlot3Dmolec,data,basis=AObasis,orbitalindex=5;
DensityPlot3Dmolec,data,basis=AObasis,orbitalindex=6;
DensityPlot3Dmolec,data,basis=AObasis,orbitalindex=7;
Does Hartree-Fock produce the correct ordering?
No, the degenerate πu,x and πu,y orbitals are the HOMO, not the expected 2σg.
Now, recalculate the molecular orbitals using another correlated method, density functional theory (DFT) and list the molecular orbital occupation numbers and energies.
data2≔DensityFunctionalmolec,basis=AObasis;
data2mo_occ1..9;
data2mo_energy1..9;
DensityPlot3Dmolec,data2,basis=AObasis,orbitalindex=5;
DensityPlot3Dmolec,data2, basis=AObasis,orbitalindex=6;
DensityPlot3Dmolec,data2,basis=AObasis,orbitalindex=7;
Does DFT produce the correct ordering?
Yes! The HOMO is non-degenerate with σ-symmetry, as expected!
In this last example, we explore a particular type of electron correlation: multireference correlation, which arises when a single determinant wavefunction is not a good approximation to the true wavefunction. Rather, more than one determinant must be included in the reference wavefunction to account for contributions from (nearly) degenerate electron configurations.
Consider the trimethylenemethane (TMM) (C4H6) diradical:
Figure 2: TMM diradical (Jorgi Stolfi, CC BY-SA 3.0)
TMM has 4 π-electrons in a π-system composed of 4 C2pz orbitals. Based on Huckel Theory, one would expect the following molecular orbital energy diagram for TMM with D3 h symmetry:
Figure 3: Trimethylenemethane diradical molecular orbital diagram (Copyright RDMCHEM LLC 2019)
Use the Maple input below to calculate the molecular orbitals of TMM using Hartree-Fock and the 6-31g basis:
withQuantumChemistry:withplots:Digits≔15:
FileTools:-Exists("TMM.xyz");
false
molec≔ReadXYZTMM.xyz;
molec≔C,−0.93706970,1.14088611,−0.41334155,C,−1.19654514,−0.27351694,−0.24674748,C,−0.09444354,−1.19591357,−0.07493269,H,0.92734618,−0.83696815,−0.08216846,H,−0.28064807,−2.25607106,0.04697584,C,−2.55791700,−0.76491556,−0.25258610,H,−2.75944612,−1.82309142,−0.13812208,H,−3.38868420,−0.07727560,−0.35264345,H,0.07403739,1.52392775,−0.35060242,H,−1.75113596,1.82832517,−0.60781581
PlotMoleculemolec
HFdata≔HartreeFockmolec,basis=AObasis;
HFdatamo_occ10..19
HFdatamo_energy10..19
DensityPlot3Dmolec,HFdata,basis=AObasis,orbitalindex=16;
An inspection of the HF MOs indicates that orbitalindex = 14 corresponds to the 1b1 orbital (see Figure 2 for symmetry label) with all C2pz orbitals in phase. The HOMO is orbitalindex = 15, which corresponds to the 2b1 orbital and is doubly occupied. The LUMO then is orbitalindex = 16, which corresponds to the 1a2 orbital.
molec_spin≔1:
dataHF2≔HartreeFockmolec,basis=AObasis,spin=2⋅molec_spin;
dataHF2mo_occ10..19;
dataHF2mo_energy10..19;
#RDMdata≔Variational2RDMmolec,basis=AObasis,spin=2⋅ molec_spin,active=2,2,casscf=true;
RDMdata≔Variational2RDMmolec,basis=AObasis,spin=0,active=2,2,casscf=true;
DensityPlot3Dmolec,HFdata,basis=AObasis,orbitalindex=16
HMM, it seems the correlation energy is essentially 0, so UHF is good enough. Multireference is not really the issue here?
Download Help Document