Mathematics Laboratory II. 1999

The Hartree-Fock self-consistent field approximation. Calculation of the structure of (Subscript[H, 2]).

This problem should be handed in to your physical chemistry tutors no later than 12.00 on 11/12/99. They will pass the problem on to me for marking. The marks count towards your course work mark for the year. It consists of three tasks each being worth  a maximum of 30 marks.

There will be  maths lab surgeries from 11.00 to 13.00 on Tuesday 9th and Tuesday 30th November. These sessions are reserved for dealing with syntax and computer problems, WHICH YOU HAVE ALREADY ENCOUNTERED. You will NOT be allowed to occupy a computer and work on the problem that session. This means that you should make a serious attempt at most of the problem BEFORE  the surgery .

The notebook containing the problem and the Mathematica code is called SCF Problem V1999.nb and can be found in the Start menu or on the desktop of the PCs and Macs.

Overview

There are three tasks in this problem. In the first task you will explore setting up suitable basis functions for the SCF calculation. In the second task you will minimise the total energy of the molecule for a given nuclear position by allowing the wavefunction of the electrons to vary until the minimum energy is found. By repeatedly employing this minimisation procedure at a range of interatomic separations you will be able to find the inter-atomic separation at which the molecule has minimum energy, and therby determine the bond length in molecular hydrogen. The third and final part requires visualising the molecular orbitals and determining how spherical they are.

There are a number of mathematical procedures that make up a complete SCF protocol. One of the most common things that you will be doing with Mathematica is multiplying matrices together and finding the eigenvectors and eigenvalues of matrices. Look back at your mathematics lecture notes to remind you what these procedures do. The matrices are required because we are constructing our molecular orbital wavefunctions from linear combinations of one electron atomic orbitals. In order to do this we start with a set of atomic wavefunctions called the basis set.  A wavefunction basis set is just a set of  wavefunctions which will be able to represent the molecular orbitals for the molecule if they are combined in the correct proportions. Solving the SCF problem is equivalent to finding out what those proportions are for our particular basis set.

Everything that you will be doing has been covered in your quantum mechanics or mathematics lectures with one exception. In this set of problems you will come across the Fock operator (usually written as \(TraditionalForm\`F\&^\)). The Fock operator essentially calculates the energy of one electron in the presence of all of the other electrons in the molecule.

Quantum mechanics: a brief reminder

For all but the simplest problems, the Schrodinger equation

H Ψ = E Ψ

cannot be solved exactly. Mostly, therefore, we have to be satisfied with approximate solutions. For multi-particle problems (such as the analysis of molecules), the Hartree-Fock method is one of the better ways of finding these.

It makes use of the fact that the true solution, Ψ((Subscript[r, 1]),(Subscript[r, 2]),...  ,(Subscript[r, n])), is the function that minimises the expected energy, given by

((∫(SuperStar[Ψ])H Ψ\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2])...)/(∫\(Ψ\^*\)Ψ\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2])...))

It follows that if we guess a solution, say Φ((Subscript[r, 1]),(Subscript[r, 2]),...  ,(Subscript[r, n])), and calculate

((∫(SuperStar[Φ])H Φ\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2])...)/(∫(SuperStar[Φ])Φ\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2])...)),

then the lower this value the better the guess (this is the variational principle).

In this assessed problem, you study the H\(TraditionalForm\`\_2\) molecule (more strictly, the Born-Oppenheimer approximation to the H\(TraditionalForm\`\_2\) molecule, in which the kinetic energy of the nuclei is held to be negligible and their positions fixed). For simplicity of calculation, we use atomic units, in which many of the constants we need (Planck's h-bar constant, the electronic mass, the electronic charge) have the value 1. In these units, the Hamiltonian (with energy being measured in hartrees, where 1 hartree= 4.3598×(10^(-18)) joules) is

-(1/2)((Subscript[∇, 1])^2)-(1/2)((Subscript[∇, 2])^2)-(1/(|(Subscript[r, 1])-(Subscript[R, 1])|))-(1/(|(Subscript[r, 2])-(Subscript[R, 1])|))-(1/(|(Subscript[r, 1])-(Subscript[R, 2])|))-(1/(|(Subscript[r, 2])-(Subscript[R, 2])|))+(1/(|(Subscript[r, 1])-(Subscript[r, 2])|))+(1/(|(Subscript[R, 1])-(Subscript[R, 2])|)),

where r\(TraditionalForm\`\_1\) and r\(TraditionalForm\`\_2\) are the position vectors of the electrons (measured in bohrs, where 1 bohr = 5.29×(10^(-11)) metres), and (Subscript[R, 1]) and (Subscript[R, 2]) are those of the nuclei. As part of the approximation, we concentrate on the electronic Hamiltonian during the minimisation process, and bring in the contribution of the internuclear distance (a distance which, of course, we're assuming to be constant) later. The electronic Hamiltonian is

(Subscript[H, el])=-(1/2)((Subscript[∇, 1])^2)-(1/2)((Subscript[∇, 2])^2)-(1/(|(Subscript[r, 1])-(Subscript[R, 1])|))-(1/(|(Subscript[r, 2])-(Subscript[R, 1])|))-(1/(|(Subscript[r, 1])-(Subscript[R, 2])|))-(1/(|(Subscript[r, 2])-(Subscript[R, 2])|))+(1/(|(Subscript[r, 1])-(Subscript[r, 2])|)).

This Hamiltonian has three components, (Subscript[H, el]) = (Subscript[H, ke]) + (Subscript[H, en]) + (Subscript[H, ee]).
The kinetic energy:

(Subscript[H, ke])=-(1/2)((Subscript[∇, 1])^2)-(1/2)((Subscript[∇, 2])^2);

the electron-nuclear attraction:

(Subscript[H, en])=-(1/(|(Subscript[r, 1])-(Subscript[R, 1])|))-(1/(|(Subscript[r, 2])-(Subscript[R, 1])|))-(1/(|(Subscript[r, 1])-(Subscript[R, 2])|))-(1/(|(Subscript[r, 2])-(Subscript[R, 2])|));

and the electron-electron interaction:

(Subscript[H, nn])=(1/(|(Subscript[r, 1])-(Subscript[r, 2])|)).

Each of these will make a contribution to the total energy.

The true solution of the Schrodinger equation for this two-electron problem for fixed nuclei will be a wave function Ψ((Subscript[r, 1]),(Subscript[r, 2])), dependent on the positions of each of the two electrons. We make the assumption that this can be expressed as a product of single-electron wave functions, Ψ((Subscript[r, 1]),(Subscript[r, 2]))=ψ((Subscript[r, 1]))ψ((Subscript[r, 2])).  These two functions are the same: there is just one space orbital in the ground state, occupied by two electrons.

Much of the following is done to make the calculations computationally tractable, rather than for any other reason. Much of the complexity in SCF methods is due to the need to use complex mathematical trickery rather than any inherent complexity in the physical picture.

The overall strategy is to express ψ as a linear combination of two carefully chosen known basis functions, \(TraditionalForm\`ψ = \(c\_1\) ϕ\_1 + \(c\_2\) ϕ\_2\), and to find the values of (Subscript[c, 1]) and (Subscript[c, 2]) that give the lowest value for the expected energy. The best physical basis is (Subscript[ϕ, 1])(r)=(e^(-1.24| r-(Subscript[R, 1])|)), (Subscript[ϕ, 2])(r)=(e^(-1.24| r-(Subscript[R, 2])|)): a so-called Slater-type orbital centred at either nucleus. Unfortunately Slater-type orbitals are very hard to work with mathematically, and to evaluate all the integrals we need we're going to have to use Gaussian orbitals instead: that is, orbitals of the form (e^(-ζ (r^2))) rather than (e^(-ζ r)). A reasonable approximation can be achieved by using a sum of three Gaussians for each basis function: the so-called STO-3G (Slater-Type Orbital "3" Gaussian) method. In this introductory example, we're going to be doing something rather cruder, namely representing each basis function as one single Gaussian: the STO-1G method. The basis we'll be using, chosen to give the best overlap with the (physically valid) Slater orbitals, is (Subscript[ϕ, 1])(r)=(e^(-0.417| r-(Subscript[R, 1])\( | \^2\))), (Subscript[ϕ, 2])(r)=(e^(-0.417| r-(Subscript[R, 2])\( | \^2\))).

Actually, we won't be using that basis as such. The maths works out far easier if our basis is chosen to be orthonormal, so that

∫\(ϕ\_1\^*\)(Subscript[ϕ, 1])\[DifferentialD]r=∫\(ϕ\_2\^*\)(Subscript[ϕ, 2])\[DifferentialD]r=1

and

∫\(ϕ\_1\^*\)(Subscript[ϕ, 2])\[DifferentialD]r=∫\(ϕ\_2\^*\)(Subscript[ϕ, 1])\[DifferentialD]r=0.

We can achieve this by making  each basis function a linear superposition of Gaussians. You are used to the idea of making linear superpositions of the type Ψ = (Subscript[ψ, 1]) + (Subscript[ψ, 2]),  but in this case we use the slightly more complex form (Subscript[ϕ, 1])(r)=a (e^(-0.417| r-(Subscript[R, 1])\( | \^2\)))+b (e^(-0.417| r-(Subscript[R, 2])\( | \^2\))) and (Subscript[ϕ, 2])(r)=c (e^(-0.417| r-(Subscript[R, 1])\( | \^2\)))+d (e^(-0.417| r-(Subscript[R, 2])\( | \^2\))), where a, b, c and d are calculated by a method shown below (see Task 1) to produce orthonormal wavefunctions.

In summary, then: we're proposing a solution of the form Ψ=ψ((Subscript[r, 1]))ψ((Subscript[r, 2])), where \(TraditionalForm\`ψ = \(c\_1\) ϕ\_1 + \(c\_2\) ϕ\_2\) and (Subscript[ϕ, 1])(r)=a (e^(-0.417| r-(Subscript[R, 1])\( | \^2\)))+b (e^(-0.417| r-(Subscript[R, 2])\( | \^2\))); (Subscript[ϕ, 2])(r)=c (e^(-0.417| r-(Subscript[R, 1])\( | \^2\)))+d (e^(-0.417| r-(Subscript[R, 2])\( | \^2\))). Therefore, in calculating the  integral

∫\(Ψ\^*\)Ψ\[DifferentialD]r,

(known as the overlap integral) we'll need to evaluate integrals of the form

∫(Subscript[ϕ, i])((Subscript[r, 2]))(Subscript[ϕ, j])((Subscript[r, 1]))(Subscript[ϕ, k])((Subscript[r, 1]))(Subscript[ϕ, l])((Subscript[r, 2]))\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2]),

where i, j, k and l can take on the value 1 or 2. (Note that the complex conjugate isn't needed, since all the wave functions in SCF can be taken as real). Similarly, the kinetic energy, nuclear attraction and electronic interaction components will give us integrals of the form

-(1/2)∫(Subscript[ϕ, i])((Subscript[r, 2]))(Subscript[ϕ, j])((Subscript[r, 1]))((Subscript[∇, 1])^2)(Subscript[ϕ, k])((Subscript[r, 1]))(Subscript[ϕ, l])((Subscript[r, 2]))\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2])    (etc),

-∫(((Subscript[ϕ, i])((Subscript[r, 2]))(Subscript[ϕ, j])((Subscript[r, 1]))(Subscript[ϕ, k])((Subscript[r, 1]))(Subscript[ϕ, l])((Subscript[r, 2])))/(|(Subscript[r, 1])-(Subscript[R, a])|))\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2])     (etc)

and

∫(((Subscript[ϕ, i])((Subscript[r, 2]))(Subscript[ϕ, j])((Subscript[r, 1]))(Subscript[ϕ, k])((Subscript[r, 1]))(Subscript[ϕ, l])((Subscript[r, 2])))/(|(Subscript[r, 1])-(Subscript[r, 2])|))\[DifferentialD](Subscript[r, 1])\[DifferentialD](Subscript[r, 2])

respectively. You should recognise these integrals from your notes on quantum mechanics. All but the last of them can be split into the product of two single-electron integrals, whose values are given by known formulae (which were not given to you in your notes). Even the non-separable electron-electron interaction integral can be found by a known formula.

The Mathematica code for evaluating the various contributions to the total energy is given below, for the non-normalized Gaussian basis functions. You do not need to investigate this section of code unless you want to; it uses standard integral formulae from quantum chemistry textbooks.

Initialisation code: block 1

\( (*\ 
      Miscellaneous\ constants\ and\ functions\ used\ by\ the\ SCF\ procedure\ *) 
\(α = 0.417;\)\)
sameSignQ[{x_, y_}] := TrueQ[x y  > 0];
f0[t_] = (1/(Sqrt[t])) Erf[(Sqrt[t])];
f0[0] = Limit[f0[t], t -> 0];
\(\(f0[0. ] = Limit[f0[t], t -> 0];\)

 (*\ 
      The\ \(integrals : \ 
          non - normalized\ Gaussian\ basis\ functions\)\ *) \)
overlap1[i_, j_] := ((((π/(2  α))))^(3/2)) 
          Exp[-(α/2) (rN[i] - rN[j]) . (rN[i] - 
                    rN[j])];
kineticEnergy1[i_, 
          j_] := (α/2) (3 - α (rN[i] - 
                      rN[j]) . (rN[i] - 
                      rN[j])) ((((π/(2  α))))^(3/2)) 
          Exp[-(α/2) (rN[i] - rN[j]) . (rN[i] - 
                    rN[j])];
\(\(electronNuclearInteraction1[i_, j_, 
          a_] := 
		\((rP = \(\(\ \)\(rN[i]\  + \ rN[j]\)\)\/2; 
		\(-\(π\/α\)\)\ Exp[\(-\(α\/2\)\) \((rN[i] - 
                      rN[j])\) . \((rN[i] - rN[j])\)] zN[a] 
            f0[2  α \((rP - rN[a])\) . \((rP - 
                      rN[a])\)])\);\)\)
\(\(electronElectronInteraction1[i_, j_, k_, 
          l_] := 
			\((rP = \(rN[i]\  + \ rN[j]\)\/2; 
          rQ = \(rN[k] + \ rN[l]\)\/2; 
			\(π\^\(5/2\)\/\(4  α\^\(5/2\)\)\) 
            Exp[\(-\(α\/2\)\) \((rN[i] - rN[j])\) . \((rN[i] - 
                        rN[j])\) - \(α\/2\) \((rN[k] - 
                        rN[l])\) . \((rN[k] - rN[l])\)] 
            f0[α \((rP - rQ)\) . \((rP - rQ)\)])\);\)\)

Orthonormalization

∫\(Ψ\^*\)Ψ\[DifferentialD]r,

∫(Subscript[ϕ, 1])(r)(Subscript[ϕ, 2])(r)\[DifferentialD]r,

\(\(\(
\)\(orthonormalizationMatrix := {{1\/\@overlap1[1, 1], 
          0}, {\(\(-overlap1[1, 2]\)\/overlap1[1, 1]\) 
            1\/\@\(overlap1[2, 2] - overlap1[1, 2]\^2\/overlap1[1, 1]\), 
          1\/\@\(overlap1[2, 2] - overlap1[1, 2]\^2\/overlap1[1, 1]\)}}\)\)\)

It's not hard to show that for the modified basis the overlap integrals are given by

overlap2[i_, j_] := 
      Sum[overlap1[μ, ν] orthonormalizationMatrix[[i, μ]] 
          orthonormalizationMatrix[[j, ν]], {μ, 1, 2}, {ν, 1, 
          2}]

Task 1: Comparing basis sets (20 marks)

R = 1.0;
rN[1] = {0, 0, R/2};
rN[2] = {0, 0, -R/2};

Initialisation code: block 2

We need to define energy integrals for the modified, orthonormal basis. This is done in much the same way as for the overlap integral; you needn't investigate this code in detail.

kineticEnergy2[i_, j_] := 
        Sum[kineticEnergy1[μ, ν] 
            orthonormalizationMatrix[[i, μ]] 
            orthonormalizationMatrix[[j, ν]], {μ, 1, 2}, {ν, 
            1, 2}];
electronNuclearInteraction2[i_, j_, a_] := 
        Sum[electronNuclearInteraction1[μ, ν, a] 
            orthonormalizationMatrix[[i, μ]] 
            orthonormalizationMatrix[[j, ν]], {μ, 1, 2}, {ν, 
            1, 2}];
electronElectronInteraction2[i_, j_, k_, l_] := 
        Sum[electronElectronInteraction1[μ, ν, ξ, η] 
            orthonormalizationMatrix[[i, μ]] 
            orthonormalizationMatrix[[j, ν]] 
            orthonormalizationMatrix[[k, ξ]] 
            orthonormalizationMatrix[[l, η]], {μ, 1, 2}, {ν, 
            1, 2}, {ξ, 1, 2}, {η, 1, 2}];

We're now ready to go.

Worked example: R = 1.0

Up to now  you have examined the contribution of part of the hamiltonian to the overall molecular energy. The minimum energy of the molecule will correspond to the correct electron distribution, but how to find it? The answer is to differentiate the expression for the energy with respect to the coefficients of the basis wavefunctions (in this case orthonormalised atomic orbitals of hydrogen atoms). This wavefunction is then used as the starting point for recalculating the energies, finding the minimum.....and so on until a global minimum is found. This iterative procedure lies at the heart of the Self-Consistent Field (SCF) method. The SCF method consists of the following stages. First, we choose where to put our nuclei. Clearly, the best thing would be to put them one bond lenth apart, but since one of the tasks in this problem is to calculate this bond length, we won't do that. Instead, for the purposes of this worked example, we'll choose 1 bohr as a "reasonable" value.

Next, we make an initial guess as to the values of (Subscript[c, 1]) and (Subscript[c, 2]) (the values (Subscript[c, 1]) = (Subscript[c, 2]) = 1/(Sqrt[2]) usually work quite well). Using these values, we set up a two-by-two matrix representing the so-called Fock operator, \(TraditionalForm\`\(\(F . \)\&^\)\) The Fock operator is closely related to the full Hamiltonian, but whereas this operates on the full wave function Ψ((Subscript[r, 1]),(Subscript[r, 2])), the Fock operator is used to calculate the energy of the one-electron orbital wave function ψ(r), in the presence of the other electron and the nuclei. The matrix that we set up (known as the Fock matrix, F) represents the effect of the Fock operator on our basis functions (Subscript[ϕ, 1]) and (Subscript[ϕ, 2]):

(Subscript[F, ij])=\(∫\(ϕ\_i\) \(F\&^\) \(ϕ\_j\) 
              d\)r.

One very significant fact about the Fock operator (and therefore the Fock matrix) is that because it takes into account the other electron's wave function, it depends on the values of (Subscript[c, 1]) and (Subscript[c, 2]), and they're exactly what we're trying to find. This is why we need an initial guess for these values: a guess that we gradually refine as follows.

We use the Fock matrix to calculate a formula, in terms of (Subscript[c, 1]) and (Subscript[c, 2]), for the energy of one electron in the presence of the other. We then minimise this energy, subject to the constraint ((Subscript[c, 1])^2) + ((Subscript[c, 2])^2) = 1, which gives us fresh values for (Subscript[c, 1]) and (Subscript[c, 2]). These fresh values can then be used as a new "guess", and so the process continues. We stop when the values of (Subscript[c, 1]) and (Subscript[c, 2]) no longer change; at this point, when the field is "self-consistent" and further attempts to minimise the expected energy lead to no improvement, we have arrived at our approximate solution and can calculate, for example, the total energy of the system.

Start by setting up two nuclei, distance 1 bohr apart, and making an initial guess as to the values of (Subscript[c, 1]) and (Subscript[c, 2]):

R = 1.0;
rN[1] = {0, 0, R/2};
rN[2] = {0, 0, -R/2};
zN[1] = 1;
zN[2] = 1;
coeffs = {(1/(Sqrt[2])), (1/(Sqrt[2]))};

We'll now set up our Fock operator. The following code sets up  a combined expression for the kinetic energy and electron-nuclear interaction terms (which we're calling core), and expressions for the two main contributors to the electronic interaction energy (known as the Coulomb and exchange integrals; see your Quantum Mechanics notes). Finally, the Fock operator itself is defined; note its dependence on the value of the coefficients (Subscript[c, 1]) and (Subscript[c, 2]).

Clear[core, twoElectron];
core[i_, j_] := core[i, j] = 
          Expand[kineticEnergy2[i, j] + 
              electronNuclearInteraction2[i, j, 1] + 
              electronNuclearInteraction2[i, j, 2]];
twoElectron[i_, j_, k_, l_] := twoElectron[i, j, k, l] = 
          electronElectronInteraction2[i, j, k, l];
coulomb[i_, j_, coeffs_] := 
        Sum[coeffs[[μ1]] coeffs[[μ2]] 
            twoElectron[i, j, μ1, μ2], {μ1, 1, 2}, {μ2, 1, 
            2}];
exchange[i_, j_, coeffs_] := 
        Sum[coeffs[[μ1]] coeffs[[μ2]] 
            twoElectron[i, μ1, j, μ2], {μ1, 1, 2}, {μ2, 1, 
            2}];
fock[i_, j_, coeffs_] := 
        Expand[core[i, j] + 2  coulomb[i, j, coeffs] - 
            exchange[i, j, coeffs]];

Mathematica note: a standard programming trick has been used in the definitions of core and twoElectron. When we define a function whose values remain constant throughout a run and are used repeatedly, execution time is speeded up several-fold if we type, for example, 
core[i_, j_] := core[i, j] = ...
instead of simply
core[i_, j_] := ... .
This is because Mathematica stores all  the values of core[1, 1],  core[1, 2] etc after the first calculation, and from then on it accesses the stored values whenever they're called instead of recalculating them each time. If we want at any stage to  reset these values we must Clear them first.

Repeating block (calculus method)

Now for the repeating block. First, we calculate the Fock matrix, F, from the guessed coefficients. Then, the quantity we are trying to minimise--the energy of one electron in the presence of the other--is equal to

\(TraditionalForm
      \`\((\(c\_1\^2\) F\_11 + \(c\_1\) \(c\_2\) F\_12 + 
            \(c\_2\) \(c\_1\) F\_21 + \(c\_2\^2\) F\_22)\)/
        \((c\_1\^2 + c\_2\^2)\)\).

(Note the structural similarity between this, the total one-electron energy, and the expression that you saw in Task 1 for the contribution to it due to electron-nuclear interaction 1.) It can be expressed more succinctly as

(((({{(Subscript[c, 1]), (Subscript[c, 2])}}))F(({{(Subscript[c, 1])}, {(Subscript[c, 2])}})))/((({{(Subscript[c, 1]), (Subscript[c, 2])}}))(({{(Subscript[c, 1])}, {(Subscript[c, 2])}})))).

(Remember, though, that Mathematica doesn't distinguish between row and column vectors.) Having calculated this expression, we then work out its derivatives with respect to (Subscript[c, 1]) and (Subscript[c, 2]), and thus minimise it subject to the constraint ((Subscript[c, 1])^2) + ((Subscript[c, 2])^2) = 1.

fockMatrix = 
        Table[N[fock[i, j, coeffs]], {i, 1, 2}, {j, 1, 2}];
Clear[c1, c2, energy];
\(\(energy[c1_, 
          c2_] = {c1, c2} . fockMatrix . {c1, c2}\/{c1, c2} . {c1, c2};\)\)
deriv1 = Together[D[energy[c1, c2], c1]];
deriv2 = Together[D[energy[c1, c2], c2]];
solutions = 
      NSolve[{deriv1 == 0, deriv2 == 0, (c1^2) + (c2^2) == 1}, {c1, 
          c2}]

A quick check of the associated energies...

energies = energy[c1, c2] /. solutions

... shows that it's the "like sign" term we want: the others represent an energy maximum. Either one will do:

coeffs = {c1, c2} /. solutions[[4]]

At this point, you should go back to the start of the repeating block, and execute it all again (all of it, including the last two input cells). Keep doing this until the values of the coefficients converge. (NB: you'll need to keep a record, either in Mathematica or on paper, of what happens to the coefficients.)

Repeating block (linear algebra method)

The whole business of differentiating and solving can be streamlined into a simple bit of linear algebra: it turns out that the coefficients, (Subscript[c, 1]) and (Subscript[c, 2]), that we get from the minimisation process are precisely the components of the "same sign" eigenvector of the Fock matrix.

Here is an alternative repeating block, which should show exactly the same behaviour as the one above (if you want to test it, you'll need to reset coeffs to its initial guessed value).

fockMatrix = 
        Table[N[fock[i, j, coeffs]], {i, 1, 2}, {j, 1, 2}];
coeffMatrix = Eigenvectors[fockMatrix];
coeffs = Abs[Select[coeffMatrix, sameSignQ][[1]]]

Total energy

The following gives the total expected energy (including the nuclear repulsion component) associated with the approximation you've just calculated. Naively, we might expect this to be twice the one-electron energy, plus the contribution of the nuclei. However, if we double the total one-electron energy, we end up counting the electron-electron contribution twice. To avoid this, we add the total one-electron energy to the core component, from which the electron-electron component is excluded, and then add the nucleus-nucleus repulsion term. Following the Born-Oppenheimer approximation, we haven't taken account of nucleus-nucleus interactions until now, but they obviously affect the total energy.

totalEnergy = 
      Sum[coeffs[[i]] 
            coeffs[[j]] (core[i, j] + fock[i, j, coeffs]), {i, 
            2}, {j, 2}] + 
        (1/Sqrt)[(rN[1] - rN[2]) . (rN[1] - rN[2])]

This is referred to as the ground state Hartree-Fock energy.

Task 2: bond length and dissociation energy. (40 Marks)

Task 3: plots of electron density. (30 Marks)

gaussians = {(E^(-0.417 ({x, y, z} - rN[1]) . ({x, y, z} - rN[1]))), 
        (E^(-0.417 ({x, y, z} - rN[2]) . ({x, y, z} - rN[2])))}

basisFunctions = orthonormalizationMatrix . gaussians

ψ = Expand[coeffs . basisFunctions]

Needs["Graphics`ContourPlot3D`"];

ContourPlot3D[ψ, {x, -1.5, 1.5}, {y, -1.5, 
        1.5}, {z, -1.5, 1.5}, Contours -> {0.2}, Axes -> True, 
      PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-1.5, 
            1.5}}]

CODA: Real time 3D graphics

<< Experimental`

<< RealTime3D`

<< Default3D`