HartreeFockMethod - Maple Help

All Products Maple MapleSim

Home : Support : Online Help : Toolboxes : Quantum Chemistry : Courses : Lessons : HartreeFockMethod

Hartree-Fock Method

Overview

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.

Hartree-Fock Method

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 &semi;$

$Digits ≔ 15$

(2.1.1)

>	$with (QuantumChemistry) &semi;$

$[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]$

(2.1.2)

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 ≔ map (x \to x / 10 + 1.11, [fsolve (expand (ChebyshevT (6, x)))]) &semi;$

$bond_distances ≔ [1.01340742, 1.03928932, 1.08411810, 1.13588190, 1.18071068, 1.20659258]$

(2.2.1.1)

We define a list of molecular geometries with each geometry corresponding to one of the bond distances

>	$molecules ≔ [seq ([[C, 0, 0, 0], [O, 0, 0, R]], R in bond_distances)] &semi;$

$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]]]$

(2.2.1.2)

The energies for each geometry may be then readily computed with the Energy command in the Quantum Chemistry package.

>	$energies ≔ [seq (Energy (molecule, method = HartreeFock, basis = cc-pVDZ), molecule in molecules)] &semi;$

$energies ≔ [−112.71803406, −112.73395669, −112.74818252, −112.74843700, −112.73846354, −112.72948290]$

(2.2.1.3)

We use polynomial interpolation to generate a polynomial in the bond distance R

>	$pes ≔ interp (bond_distances, energies, R) &semi;$

$pes ≔ - 11.82147324 R^{5} + 74.99529529 R^{4} - 193.48401781 R^{3} + 254.30676890 R^{2} - 169.92125737 R - 66.78245246$

(2.2.1.4)

The potential energy surface (curve) can be plotted

>	$plot (pes (R), R = 1.025 .. 1.2, axes = boxed, labels = [' R',' E'], color = red, thickness = 3) &semi;$

Finally, we differential the potential energy curve with respect to R and set the derivative to zero.

>	$eq ≔ diff (pes, R) = 0 &semi;$

$eq ≔ - 59.10736619 R^{4} + 299.98118117 R^{3} - 580.45205343 R^{2} + 508.61353779 R - 169.92125737 = 0$

(2.2.1.5)

Solving the resulting equation yields the equilibrium bond length

>	$R_eq ≔ fsolve (eq, R = 1.0 .. 1.2) &semi;$

$R_eq ≔ 1.11014350$

(2.2.1.6)

(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]] &semi;$

$molecule ≔ [[C, 0, 0, 0], [O, 0, 0, 1.11014350]]$

(2.2.2.1)

Using the Dipole command, we can compute the dipole from the Hartree-Fock method in the cc-pVDZ basis set

>	$Dipole (molecule, method = HartreeFock, basis = cc-pVDZ) &semi;$

$[\begin{array}{c} X & 0. \\ Y & 0. \\ Z & −0.14581806 \end{array}]$

(2.2.2.2)

(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:

>	$AtomicData (C) &semi;$

$table ([density = 2.20000000 ⟦\frac{g}{{cm}^{3}}⟧, atomicnumber = 6, name = carbon, names = \{carbon\}, atomicweight = 12.01070000 ⟦amu⟧, ionizationenergy = 11.26030000 ⟦eV⟧, electronegativity = 2.55000000, symbol = C, electronaffinity = 1.26290000 ⟦eV⟧])$

(2.2.2.3)

>	$AtomicData (O) &semi;$

$table ([atomicnumber = 8, meltingpoint = 54.36000000 ⟦K⟧, name = oxygen, boilingpoint = 90.20000000 ⟦K⟧, names = \{oxygen\}, atomicweight = 15.99940000 ⟦amu⟧, ionizationenergy = 13.61810000 ⟦eV⟧, electronegativity = 3.44000000, symbol = O, electronaffinity = 1.46111030 ⟦eV⟧])$

(2.2.2.4)

Nitrogen Molecule

We compute the equilibrium bond length, the potential energy curve, and the ionization energy of the nitrogen molecule.

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 ≔ map (x \to x / 10 + 1.11, [fsolve (expand (ChebyshevT (6, x)))]) &semi;$

$bond_distances ≔ [1.01340742, 1.03928932, 1.08411810, 1.13588190, 1.18071068, 1.20659258]$

(2.3.1.1)

We define a list of molecular geometries with each geometry corresponding to one of the bond distances

>	$molecules ≔ [seq ([[N, 0, 0, 0], [N, 0, 0, R]], R in bond_distances)] &semi;$

$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]]]$

(2.3.1.2)

The energies for each geometry may be then readily computed with the Energy command in the Quantum Chemistry package.

>	$energies ≔ [seq (Energy (molecule, method = HartreeFock, basis = cc-pVDZ), molecule in molecules)] &semi;$

$energies ≔ [−108.93853400, −108.94988602, −108.95539406, −108.94471122, −108.92485656, −108.91010098]$

(2.3.1.3)

We use polynomial interpolation to generate a polynomial in the bond distance R

>	$pes ≔ interp (bond_distances, energies, R) &semi;$

$pes ≔ - 11.02261081 R^{5} + 71.10759569 R^{4} - 186.63332264 R^{3} + 249.43991316 R^{2} - 169.02328246 R - 62.79813369$

(2.3.1.4)

The potential energy surface (curve) can be plotted

>	$plot (pes (R), R = 1.025 .. 1.2, axes = boxed, labels = [' R',' E'], color = red, thickness = 3) &semi;$

Finally, we differential the potential energy curve with respect to R and set the derivative to zero.

>	$eq ≔ diff (pes, R) = 0 &semi;$

$eq ≔ - 55.11305403 R^{4} + 284.43038276 R^{3} - 559.89996793 R^{2} + 498.87982632 R - 169.02328246 = 0$

(2.3.1.5)

Solving the resulting equation yields the equilibrium bond length

>	$R_eq ≔ fsolve (eq, R = 1.0 .. 1.2) &semi;$

$R_eq ≔ 1.07729844$

(2.3.1.6)

(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 N₂ 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 ≔ map (x \to x / 3 + 1.33, [fsolve (expand (ChebyshevT (6, x)))]) &semi;$

$bond_distances ≔ [1.00802472, 1.09429774, 1.24372698, 1.41627302, 1.56570226, 1.65197528]$

(2.3.2.1)

We define a list of molecular geometries with each geometry corresponding to one of the bond distances

>	$molecules ≔ [seq ([[N, 0, 0, 0], [N, 0, 0, R]], R in bond_distances)] &semi;$

$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]]]$

(2.3.2.2)

The energies for each geometry may be then readily computed with the Energy command in the Quantum Chemistry package.

>	$energies_hf ≔ [seq (Energy (molecule, method = HartreeFock, basis = cc-pVDZ), molecule in molecules)] &semi;$

$energies_hf ≔ [−108.93528947, −108.95455790, −108.88577598, −108.74797383, −108.62367247, −108.55624763]$

(2.3.2.3)

We use polynomial interpolation to generate a polynomial in the bond distance R

>	$pes_hf ≔ interp (bond_distances, energies_hf, R) &semi;$

$pes_hf ≔ - 3.57881332 R^{5} + 27.16588365 R^{4} - 83.00182262 R^{3} + 127.41624181 R^{2} - 97.30449614 R - 79.62695783$

(2.3.2.4)

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 ≔ [seq (Energy (molecule, method = ActiveSpaceSCF, basis = cc-pVDZ, active = [6, 6]), molecule in molecules)] &semi;$

$energies_ci ≔ [−109.05151002, −109.08966958, −109.05839964, −108.97150902, −108.89798298, −108.86298449]$

(2.3.2.5)

We use polynomial interpolation to generate a polynomial in the bond distance R

>	$pes_ci ≔ interp (bond_distances, energies_ci, R) &semi;$

$pes_ci ≔ - 3.65343945 R^{5} + 27.63278366 R^{4} - 84.16563811 R^{3} + 128.72416429 R^{2} - 98.12886339 R - 79.45353462$

(2.3.2.6)

The potential energy curves from the Hartree-Fock method (red) and the more advanced method (blue) can be plotted together

>	$p_hf ≔ plot (pes_hf, R = 1.0 .. 1.65, axes = boxed, labels = [' R',' E'], color = red, thickness = 3) &colon; p_ci ≔ plot (pes_ci, R = 1.0 .. 1.65, axes = boxed, labels = [' R',' E'], color = blue, thickness = 3) &colon; plots :- display (\{p_hf, p_ci\}) &semi;$

(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 N₂ 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]] &semi;$

$molecule ≔ [[N, 0, 0, 0], [N, 0, 0, 1.07729844]]$

(2.3.3.1)

The energies of the Hartree-Fock molecular orbitals (MOs) can be computed with the MOEnergies command using the Hartree-Fock method

>	$mo_energies ≔ MOEnergies (molecule, method = HartreeFock, basis = cc-pVDZ, symmetry = true) &semi;$

The occupations of the Hartree-Fock MOs can similarly be computed with the MOOccupations command

>	$mo_occupations ≔ MOOccupations (molecule, method = HartreeFock, basis = dzp, symmetry = true) &semi;$

Finally, the point-group symmetries of the Hartree-Fock MOs can be computed with the MOSymmetries command

>	$mo_symmetries ≔ MOSymmetries (molecule, method = HartreeFock, basis = dzp, symmetry = true) &semi;$

We can collect all of this information in a Maple list of lists

>	$mo_data ≔ [seq ([i, mo_occupations [i], mo_symmetries [i], mo_energies [i]], i = 1 .. 10)] &semi;$

$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]]$

(2.3.3.2)

We observe that the 7^th 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 ≔ convert (abs (mo_energies [7]) &ast; Units :- Unit (hartrees), units, eV) &semi;$

$IE ≔ 16.78904012 ⟦eV⟧$

(2.3.3.3)

The electron density of the HOMO can be visualized with the DensityPlot3D command. This command requires the output from the HartreeFock command.

>	$data_hf ≔ HartreeFock (molecule, basis = cc-pVDZ, symmetry = true) &semi;$

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.

>	$DensityPlot3D (molecule, data_hf, orbitalindex = 7, basis = cc-pVDZ, symmetry = true) &semi;$

(k) Use the computed MO energies to construct an MO diagram for the N₂ 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 N₂⁺ is in a doublet, sigma state with a term symbol ²S. 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 N₂ molecule from any standard textbook.

References

1.	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"

2.	V. Fock, Z. Physik 61, 126-148 (1930). "Näherungsmethode zur lösung des quantenmechanischen mehrkörperproblems"

3.	C. C. J. Roothaan, Rev. Mod. Phys. 23, 69–89 (1951). "New developments in molecular orbital theory"

4.	A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books, New York, 1996).

Download Help Document

Maple

Maple Add-Ons

Student Success Platform

Improving Retention Rates

MapleSim

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim

Maple

Powerful math software that is easy to use

Maple Add-Ons

Student Success Platform

Improving Retention Rates

MapleSim

Advanced System Level Modeling

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim