CorrelationEnergies - Maple Help

Electron Correlation and Correlation Energy

 Overview 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:     Based on the variational principle,  is lower in energy than , so ${\mathrm{ϵ}}_{\mathrm{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:

Bond Dissociation of Hydrogen $\left({\mathbf{H}}_{\mathbf{2}}\right)$

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'.

 >
 > $\mathrm{AObasis}≔"6-31g":$
 > $\mathrm{rmin}≔0.4:$
 > $\mathrm{rmax}≔5.0:$
 > $\mathrm{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:

 >
 >
 >
 >
 >
 > 

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.

 > $\mathrm{atom}≔\left[\left["H",0,0,0\right]\right];$
 ${\mathrm{atom}}{≔}\left[\left[{"H"}{,}{0}{,}{0}{,}{0}\right]\right]$ (2.1)
 >
 ${\mathrm{dissociation}}{≔}{-0.99646582}$ (2.2)

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:

 >
 >
 >
 > $\mathrm{plots}\left[\mathrm{display}\right]\left(\mathrm{plot_HF},\mathrm{plot_MP2}\right)$

Parametric 2-RDM:

 >
 >
 >
 > $\mathrm{plots}\left[\mathrm{display}\right]\left(\mathrm{plot_HF},\mathrm{plot_P2RDM}\right)$

Configuration Interaction

 > 
 >
 >
 > $\mathrm{plots}\left[\mathrm{display}\right]\left(\mathrm{plot_HF},\mathrm{plot_CI}\right)$
 > 

To compare, combine all 4 potential energy curves in a single plot:

 > $\mathrm{plots}\left[\mathrm{display}\right]\left(\mathrm{plot_HF},\mathrm{plot_MP2},\mathrm{plot_P2RDM},\mathrm{plot_CI}\right);$

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.

Molecular Orbitals of Nitrogen (${\mathbf{N}}_{\mathbf{2}}$)

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.

 >
 > $\mathrm{molec}≔\left[\left["N",0,0,0\right],\left["N",0,0,1.09\right]\right];$
 ${\mathrm{molec}}{≔}\left[\left[{"N"}{,}{0}{,}{0}{,}{0}\right]{,}\left[{"N"}{,}{0}{,}{0}{,}{1.09000000}\right]\right]$ (3.1)
 > $\mathrm{AObasis}≔"sto-3g":$
 > $\mathrm{data}≔\mathrm{HartreeFock}\left(\mathrm{molec}\right);$
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744096972586142}\right)$ (3.2)

For  N2 with 14 electrons all paired, the HOMO is the 7th MO, as indicated by the occupation numbers of the molecular orbitals:

 > $\mathrm{data}\left[\mathrm{mo_occ}\right]\left[1..9\right];$
 ${{\mathrm{_rtable}}}_{{18446744097429920934}}$ (3.3)

List the molecular orbital energies and plot the HOMO, HOMO-1, and HOMO-2 molecular orbitals using the Maple input below:

 > $\mathrm{data}\left[\mathrm{mo_energy}\right]\left[1..9\right];$
 ${{\mathrm{_rtable}}}_{{18446744097429923230}}$ (3.4)
 > $\mathrm{DensityPlot3D}\left(\mathrm{molec},\mathrm{data},\mathrm{basis}=\mathrm{AObasis},\mathrm{orbitalindex}=5\right);$
 > $\mathrm{DensityPlot3D}\left(\mathrm{molec},\mathrm{data},\mathrm{basis}=\mathrm{AObasis},\mathrm{orbitalindex}=6\right);$
 > $\mathrm{DensityPlot3D}\left(\mathrm{molec},\mathrm{data},\mathrm{basis}=\mathrm{AObasis},\mathrm{orbitalindex}=7\right);$

Does Hartree-Fock produce the correct ordering?

 Answer 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.

 > $\mathrm{data2}≔\mathrm{DensityFunctional}\left(\mathrm{molec},\mathrm{basis}=\mathrm{AObasis}\right);$
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744097081198942}\right)$ (3.5)
 > $\mathrm{data2}\left[\mathrm{mo_occ}\right]\left[1..9\right];$
 ${{\mathrm{_rtable}}}_{{18446744096959174526}}$ (3.6)
 > $\mathrm{data2}\left[\mathrm{mo_energy}\right]\left[1..9\right];$
 ${{\mathrm{_rtable}}}_{{18446744096959197046}}$ (3.7)
 > $\mathrm{DensityPlot3D}\left(\mathrm{molec},\mathrm{data2},\mathrm{basis}=\mathrm{AObasis},\mathrm{orbitalindex}=5\right);$
 >
 > $\mathrm{DensityPlot3D}\left(\mathrm{molec},\mathrm{data2},\mathrm{basis}=\mathrm{AObasis},\mathrm{orbitalindex}=7\right);$

Does DFT produce the correct ordering?

 Answer Yes!   The HOMO is non-degenerate with σ-symmetry, as expected!

Trimethylenemethane (TMM)

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) (${\mathrm{C}}_{4}{\mathrm{H}}_{6}$) 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  symmetry:

Use the Maple input below to calculate the molecular orbitals of TMM using Hartree-Fock and the 6-31g basis:

 >
 > FileTools:-Exists("TMM.xyz");
 ${\mathrm{false}}$ (4.1)
 > $\mathrm{molec}≔\mathrm{ReadXYZ}\left("TMM.xyz"\right);$
 ${\mathrm{molec}}{≔}\left[\left[{"C"}{,}{-0.93706970}{,}{1.14088611}{,}{-0.41334155}\right]{,}\left[{"C"}{,}{-1.19654514}{,}{-0.27351694}{,}{-0.24674748}\right]{,}\left[{"C"}{,}{-0.09444354}{,}{-1.19591357}{,}{-0.07493269}\right]{,}\left[{"H"}{,}{0.92734618}{,}{-0.83696815}{,}{-0.08216846}\right]{,}\left[{"H"}{,}{-0.28064807}{,}{-2.25607106}{,}{0.04697584}\right]{,}\left[{"C"}{,}{-2.55791700}{,}{-0.76491556}{,}{-0.25258610}\right]{,}\left[{"H"}{,}{-2.75944612}{,}{-1.82309142}{,}{-0.13812208}\right]{,}\left[{"H"}{,}{-3.38868420}{,}{-0.07727560}{,}{-0.35264345}\right]{,}\left[{"H"}{,}{0.07403739}{,}{1.52392775}{,}{-0.35060242}\right]{,}\left[{"H"}{,}{-1.75113596}{,}{1.82832517}{,}{-0.60781581}\right]\right]$ (4.2)
 > $\mathrm{PlotMolecule}\left(\mathrm{molec}\right)$
 > $\mathrm{AObasis}≔"sto-3g":$
 > $\mathrm{HFdata}≔\mathrm{HartreeFock}\left(\mathrm{molec},\mathrm{basis}=\mathrm{AObasis}\right);$
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744098030518718}\right)$ (4.3)
 > $\mathrm{HFdata}\left[\mathrm{mo_occ}\right]\left[10..19\right]$
 ${{\mathrm{_rtable}}}_{{18446744096949539414}}$ (4.4)
 > $\mathrm{HFdata}\left[\mathrm{mo_energy}\right]\left[10..19\right]$
 ${{\mathrm{_rtable}}}_{{18446744096949534958}}$ (4.5)
 > $\mathrm{DensityPlot3D}\left(\mathrm{molec},\mathrm{HFdata},\mathrm{basis}=\mathrm{AObasis},\mathrm{orbitalindex}=16\right);$

An inspection of the HF MOs indicates that orbitalindex = 14 corresponds to the $1{\mathrm{b}}_{1}$ orbital (see Figure 2 for symmetry label) with all C2pz orbitals in phase. The HOMO is orbitalindex = 15, which corresponds to the $2{\mathrm{b}}_{1}$ orbital and is doubly occupied. The LUMO then is orbitalindex = 16, which corresponds to the  $1{\mathrm{a}}_{2}$ orbital.

 > $\mathrm{molec_spin}≔1:$
 >
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744097077751006}\right)$ (4.6)
 > $\mathrm{dataHF2}\left[\mathrm{mo_occ}\right]\left[10..19\right];$
 ${{\mathrm{_rtable}}}_{{18446744096959324998}}$ (4.7)
 > $\mathrm{dataHF2}\left[\mathrm{mo_energy}\right]\left[10..19\right];$
 ${{\mathrm{_rtable}}}_{{18446744096959327878}}$ (4.8)
 >
 > $\mathrm{RDMdata}≔\mathrm{Variational2RDM}\left(\mathrm{molec},\mathrm{basis}=\mathrm{AObasis},\mathrm{spin}=0,\mathrm{active}=\left[2,2\right],\mathrm{casscf}=\mathrm{true}\right);$
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744097040015806}\right)$ (4.9)
 > $\mathrm{DensityPlot3D}\left(\mathrm{molec},\mathrm{HFdata},\mathrm{basis}=\mathrm{AObasis},\mathrm{orbitalindex}=16\right)$

HMM, it seems the correlation energy is essentially 0, so UHF is good enough. Multireference is not really the issue here?

 References