Hartree-Fock Method
Copyright (c) RDMCHEM LLC 2020
Overview
References
The Hartree-Fock method, developed by Daniel Hartree and Vladimir Fock, is an approximation to the many-electron Schrödinger equation in which the electrons are assumed to be statistically independent. In 1928 Hartree approximated the many-electron wave function as a product of orbitals which is equivalent to assuming statistical independence. Such a product, however, does not obey the anti-symmetrization of the wave function, first arrived by Heisenberg and Dirac in 1926. In 1930 Fock generalized Hartree's product wave function to include anti-symmetrization. The antisymmetrized wave function can be easily expressed as follows
ψ(1,2,3,..,N) = φi(1) ∧ φj(2) ∧ φk(3) ∧ ... ∧φn(N)
where each ϕ denotes a spin orbital with the index indicating the specific orbital, each roman number represents the three spatial coordinates and the spin component of an electron, and the symbol ∧ denotes the Grassmann wedge product which is the antisymmetric tensor product named after Hermann Grassmann who developed it in the 1850s in the context of exterior algebra. The antisymmetric wave function can also be expressed as a determinant of orbitals, first recognized by John Slater in the 1930s, and hence, called the Slater determinant.
The antisymmetrized product wave function can be determined by taking the excitation value of the Hamiltonian operator and minimizing the energy with respect to the orbitals. The optimized orbitals at the energy minimum are known as Hartree-Fock (molecular) orbitals. The wedge product of the lowest N spin orbitals generates a Hartree-Fock wave function. Because the wave function is a product of orbitals, the Hamiltonian can be written as an effective one-body Hamiltonian that depends upon the orbitals. The effective Hamiltonian can be solved iteratively with the orbitals updated at each iteration in an algorithm known as the self-consistent-field (SCF) method, which was first proposed and implemented by Clemens Roothan in the 1950s. Modern implementations of the Hartree-Fock method augment the SCF procedure to accelerate convergence. The Hartree-Fock method is also called a mean-field method because with the assumption of statistical independence each electron experiences the mean field of the other electrons. To improve the solution of the many-electron Schrödinger equation beyond the Hartree-Fock method, we must allow the electrons to become statistically dependent (or correlated). Methods that correlate the electrons are known as correlation methods.
After loading the Quantum Chemistry package, we explore several properties of the carbon monoxide and nitrogen molecules with the Hartree-Fock method. In some cases you will find that the properties from Hartree-Fock theory do not agree with experiment. These discrepancies reveal the need for more advanced methods for solving the Schrödinger equation.
Quantum Chemistry
We set the number of Digits to be used in computations to 15 and load the Quantum Chemistry package using Maple's with command.
Digits ≔ 15;
Digits≔15
withQuantumChemistry;
AOLabels,ActiveSpaceCI,ActiveSpaceSCF,AtomicData,BondAngles,BondDistances,Charges,ChargesPlot,CorrelationEnergy,CoupledCluster,DensityFunctional,DensityPlot3D,Dipole,DipolePlot,Energy,ExcitationEnergies,ExcitationSpectra,ExcitationSpectraPlot,ExcitedStateEnergies,ExcitedStateSpins,FullCI,GeometryOptimization,HartreeFock,Interactive,Isotopes,MOCoefficients,MODiagram,MOEnergies,MOIntegrals,MOOccupations,MOOccupationsPlot,MOSymmetries,MP2,MolecularData,MolecularGeometry,NuclearEnergy,NuclearGradient,OscillatorStrengths,Parametric2RDM,PlotMolecule,Populations,RDM1,RDM2,RTM1,ReadXYZ,Restore,Save,SaveXYZ,SearchBasisSets,SearchFunctionals,SkeletalStructure,Thermodynamics,TransitionDipolePlot,TransitionDipoles,TransitionOrbitalPlot,TransitionOrbitals,Variational2RDM,VibrationalModeAnimation,VibrationalModes,Video
Carbon Monoxide
We compute the equilibrium bond length and the dipole moment of carbon monoxide.
Equilibrium Bond Length
To compute the equilibrium bond length, we select a set of bond distances from the roots of the sixth-order Chebyshev polynomial that are suitable for interpolation
bond_distances ≔ mapx → x/10+1.11, fsolveexpandChebyshevT6,x;
bond_distances≔1.01340742,1.03928932,1.08411810,1.13588190,1.18071068,1.20659258
We define a list of molecular geometries with each geometry corresponding to one of the bond distances
molecules ≔ seqC,0,0,0,O,0,0,R, R in bond_distances;
molecules≔C,0,0,0,O,0,0,1.01340742,C,0,0,0,O,0,0,1.03928932,C,0,0,0,O,0,0,1.08411810,C,0,0,0,O,0,0,1.13588190,C,0,0,0,O,0,0,1.18071068,C,0,0,0,O,0,0,1.20659258
The energies for each geometry may be then readily computed with the Energy command in the Quantum Chemistry package.
energies ≔ seqEnergymolecule,method=HartreeFock, basis=cc-pVDZ, molecule in molecules;
energies≔−112.71803406,−112.73395669,−112.74818252,−112.74843700,−112.73846354,−112.72948290
We use polynomial interpolation to generate a polynomial in the bond distance R
pes ≔ interpbond_distances,energies,R;
pes≔−11.82147324R5+74.99529529R4−193.48401781R3+254.30676890R2−169.92125737R−66.78245246
The potential energy surface (curve) can be plotted
plotpesR, R=1.025..1.2, axes=boxed, labels='R','E', color=red, thickness=3;
Finally, we differential the potential energy curve with respect to R and set the derivative to zero.
eq ≔ diffpes, R = 0;
eq≔−59.10736619R4+299.98118117R3−580.45205343R2+508.61353779R−169.92125737=0
Solving the resulting equation yields the equilibrium bond length
R_eq ≔ fsolveeq, R=1.0..1.2;
R_eq≔1.11014350
(a) By changing the basis keyword in the Energy command, repeat the above computations for the following larger basis sets: (i) cc-pVTZ, (ii) cc-pVQZ, and (iii) cc-pV5Z, and report your results as a table. (Note that the basis sets are increasing in size as cc-pVDZ < cc-pVTZ < cc-pVQZ < cc-pV5Z)
(b) Does the bond length of CO from the Hartree-Fock method increase or decrease with larger basis sets?
Dipole Moment
Let us compute the electric dipole moment of CO from Hartree-Fock at the equilibrium bond length in the cc-pVDZ basis set. First, we define the molecule's geometry as a Maple list of lists giving the atom's names and xyz coordinates
molecule ≔ C,0,0,0,O,0,0,1.11014350;
molecule≔C,0,0,0,O,0,0,1.11014350
Using the Dipole command, we can compute the dipole from the Hartree-Fock method in the cc-pVDZ basis set
Dipolemolecule, method=HartreeFock, basis=cc-pVDZ;
X0.Y0.Z−0.14581806
(c) By changing the basis keyword in the Energy command, repeat the above computations for the following larger basis sets: (i) cc-pVTZ, (ii) cc-pVQZ, and (iii) cc-pV5Z, and report your results as a table. (d) Does the dipole moment of CO from the Hartree-Fock method increase or decrease with with larger basis sets? Does it appear to converge in the large-basis-set limit?
(e) Is the dipole moment of CO from the largest basis set in agreement with the electronegativities of carbon and oxygen?
(f) Is the dipole moment of CO from the largest basis set in agreement with the experimental value available at NIST (https://cccbdb.nist.gov/diplistx.asp)? Explain briefly.
Note that the electronegativities as well as other properties of carbon and oxygen can be obtained with the command AtomicData:
AtomicDataC;
tabledensity=2.20000000gcm3,atomicnumber=6,name=carbon,names=carbon,atomicweight=12.01070000amu,ionizationenergy=11.26030000eV,electronegativity=2.55000000,symbol=C,electronaffinity=1.26290000eV
AtomicDataO;
tableatomicnumber=8,meltingpoint=54.36000000K,name=oxygen,boilingpoint=90.20000000K,names=oxygen,atomicweight=15.99940000amu,ionizationenergy=13.61810000eV,electronegativity=3.44000000,symbol=O,electronaffinity=1.46111030eV
Nitrogen Molecule
We compute the equilibrium bond length, the potential energy curve, and the ionization energy of the nitrogen molecule.
molecules ≔ seqN,0,0,0,N,0,0,R, R in bond_distances;
molecules≔N,0,0,0,N,0,0,1.01340742,N,0,0,0,N,0,0,1.03928932,N,0,0,0,N,0,0,1.08411810,N,0,0,0,N,0,0,1.13588190,N,0,0,0,N,0,0,1.18071068,N,0,0,0,N,0,0,1.20659258
energies≔−108.93853400,−108.94988602,−108.95539406,−108.94471122,−108.92485656,−108.91010098
pes≔−11.02261081R5+71.10759569R4−186.63332264R3+249.43991316R2−169.02328246R−62.79813369
eq≔−55.11305403R4+284.43038276R3−559.89996793R2+498.87982632R−169.02328246=0
R_eq≔1.07729844
(g) By changing the basis keyword in the Energy command, repeat the above computations for the following larger basis sets: (i) cc-pVTZ, (ii) cc-pVQZ, and (iii) cc-pV5Z, and report your results as a table. (Note that the basis sets are increasing in size as cc-pVDZ < cc-pVTZ < cc-pVQZ < cc-pV5Z)
(h) Does the bond length of N2 from the Hartree-Fock method increase or decrease with larger basis sets?
Potential Energy Curve
To compute the potential energy curve, we select a set of bond distances from the roots of the sixth-order Chebyshev polynomial that are suitable for interpolation
bond_distances ≔ mapx → x/3+1.33, fsolveexpandChebyshevT6,x;
bond_distances≔1.00802472,1.09429774,1.24372698,1.41627302,1.56570226,1.65197528
molecules≔N,0,0,0,N,0,0,1.00802472,N,0,0,0,N,0,0,1.09429774,N,0,0,0,N,0,0,1.24372698,N,0,0,0,N,0,0,1.41627302,N,0,0,0,N,0,0,1.56570226,N,0,0,0,N,0,0,1.65197528
energies_hf ≔ seqEnergymolecule,method=HartreeFock, basis=cc-pVDZ, molecule in molecules;
energies_hf≔−108.93528947,−108.95455790,−108.88577598,−108.74797383,−108.62367247,−108.55624763
pes_hf ≔ interpbond_distances,energies_hf,R;
pes_hf≔−3.57881332R5+27.16588365R4−83.00182262R3+127.41624181R2−97.30449614R−79.62695783
Similarly, for comparison we can compute the energies from with the ActiveSpaceSCF command which provides a more accurate approximation to the Schrödinger equation, allowing orbitals in an active space to be treated beyond the Hartree-Fock method
energies_ci ≔ seqEnergymolecule,method=ActiveSpaceSCF,basis=cc-pVDZ,active=6,6, molecule in molecules;
energies_ci≔−109.05151002,−109.08966958,−109.05839964,−108.97150902,−108.89798298,−108.86298449
pes_ci ≔ interpbond_distances,energies_ci,R;
pes_ci≔−3.65343945R5+27.63278366R4−84.16563811R3+128.72416429R2−98.12886339R−79.45353462
The potential energy curves from the Hartree-Fock method (red) and the more advanced method (blue) can be plotted together
p_hf ≔ plotpes_hf, R=1.0..1.65, axes=boxed, labels='R','E', color=red, thickness=3:p_ci ≔ plotpes_ci, R=1.0..1.65, axes=boxed, labels='R','E', color=blue, thickness=3:plots:-displayp_hf,p_ci;
(i) Compare the equilibrium bond length from the Hartree-Fock method with the equilibrium bond length from the more advanced (correlated) method.
(j) Where does the Hartree-Fock potential energy curve have its greatest error relative to the curve from correlated method? Explain briefly why?
Ionization Energy
According to Koopman's theorem the ionization energy is equal to the absolute value of the energy of the highest occupied molecular orbital. We can approximate the ionization energy of the nitrogen molecule using Koopman's theorem. First, we define N2 at its equilibrium bond length from the Hartree-Fock method in the cc-pVDZ basis set
molecule ≔ N,0,0,0,N,0,0,1.07729844;
molecule≔N,0,0,0,N,0,0,1.07729844
The energies of the Hartree-Fock molecular orbitals (MOs) can be computed with the MOEnergies command using the Hartree-Fock method
mo_energies ≔ MOEnergiesmolecule, method=HartreeFock, basis=cc-pVDZ, symmetry=true;
The occupations of the Hartree-Fock MOs can similarly be computed with the MOOccupations command
mo_occupations ≔ MOOccupationsmolecule, method=HartreeFock, basis=dzp, symmetry=true;
Finally, the point-group symmetries of the Hartree-Fock MOs can be computed with the MOSymmetries command
mo_symmetries ≔ MOSymmetriesmolecule, method=HartreeFock, basis=dzp, symmetry=true;
We can collect all of this information in a Maple list of lists
mo_data ≔ seqi,mo_occupationsi,mo_symmetriesi,mo_energiesi, i=1..10;
mo_data≔1,2.00000000,A1g,−15.67996903,2,2.00000000,A1u,−15.67611993,3,2.00000000,A1g,−1.48632590,4,2.00000000,A1u,−0.76741734,5,2.00000000,A1g,−0.62793833,6,2.00000000,E1uy,−0.61698585,7,2.00000000,E1ux,−0.61698585,8,0.,E1gx,0.18662326,9,0.,E1gy,0.18662326,10,0.,A1u,0.60065161
We observe that the 7th MO is the highest occupied molecular orbital (HOMO), and hence, the ionization energy from Koopman's theorem is 0.617 hartrees. With Maple we can easily convert the result to another unit such as electron volts (eV). Note that the resulting energy is slightly greater than the ionization energy of the hydrogen atom (13.6 eV).
IE ≔ convertabsmo_energies7*Units:-Unithartrees,units, eV;
IE≔16.78904012eV
The electron density of the HOMO can be visualized with the DensityPlot3D command. This command requires the output from the HartreeFock command.
data_hf ≔ HartreeFockmolecule, basis=cc-pVDZ, symmetry=true;
After running the Hartree-Fock command, we run the DensityPlot3D command with the molecule's geometry, the output table from the Hartree-Fock calculation, the orbital index of the MO to be visualized, and the basis and symmetry keywords of the Hartree-Fock calculation. Click on the plot to rotate the molecule into different orientations.
DensityPlot3Dmolecule,data_hf, orbitalindex=7, basis=cc-pVDZ, symmetry=true;
(k) Use the computed MO energies to construct an MO diagram for the N2 molecule.
(l) What is the term symbol of the ionized state predicted by Koopman's theorem from the Hartree-Fock method? (m) Experimentally, the ionized nitrogen molecule N2+ is in a doublet, sigma state with a term symbol 2S. Does your result in (l) agree with experiment? (n) Explain any differences in (m) by comparing your MO diagram with the MO diagram of the N2 molecule from any standard textbook.
D. R. Hartree, Math. Proc. Camb. Phil. Soc. 24, 111-132 (1928). "The wave mechanics of an atom with a non-Coulomb central field. Part II. Some results and discussion"
V. Fock, Z. Physik 61, 126-148 (1930). "Näherungsmethode zur lösung des quantenmechanischen mehrkörperproblems"
C. C. J. Roothaan, Rev. Mod. Phys. 23, 69–89 (1951). "New developments in molecular orbital theory"
A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books, New York, 1996).
Download Help Document
