# Introduction to Quantum Chemistry

## Contents

# Introduction to Quantum Chemistry#

The main objective of this course is to write a working restricted Hartree–Fock program.

Contents

Note

Start by extracting the `course-material.zip`

which you downloaded earlier.
You should find a `scf/`

directory containing the starting code for your
program (see `scf/app/main.f90`

).

You can just continue to use fpm to work on your program now.

Alternatively, to make building your program easier with just the Fortran compiler
and give you more time to focus on the actual programming we put some build
automation there as well.
To use it just invoke `make`

and find everything else taken care for.

For more details on building software and how we automated this process you can find a guide at the fortran-lang learning page.

## Restricted Hartree–Fock#

The first task in this course is to code a working restricted Hartree–Fock program.

### Getting Input#

First, we have to gather which information is important to perform a
Hartree–Fock calculation.
We need to know about the molecular geometry. There are several ways
to represent the geometry, the most simple is to use Cartesian coordinates
as a vector of *x*, *y* and *z* coordinates.
What maybe is not that obvious, we have to decide on a unit system for
the coordinates, usually we would follow IUPAC recommendations and use
SI units like meters, but if you think about a molecule on the length
scale in meters it becomes quite inconvenient.
To get into the right ballpark we could choose Ångström which has an
obvious relation to the SI unit, but here we will use Bohr which is
the length unit of the atomic units system.
Atomic units have the advantage that a lot of constants drop out of the
equations we want to code up and we can convert them easily back,
in fact, we provide a conversion routine back to SI units for you.

These are the considerations we have to put in the geometry, now we have to identify the atoms/nuclei in our molecule. Chemical identity like the element is given by the nuclear charge and the number of electrons. It is quite popular to specify both with the element symbol which map to the atomic number (which corresponds to the nuclear charge and the number of electrons of the neutral atom) since it is the intuitive choice. We will not follow this approach since it would require you to work with character type variables, instead we will separate the nuclear charge information and the number of electrons. For each element, we will read the nuclear charge in multiples of the electron charge (which is the atomic unit of charge) and specify the number of electrons separately.

Having put some thoughts in the geometric representation of the system, we now have to tackle the electronic representation, that is, we need a basis set to expand our wavefunction. There are many possible choices, like atom centered basis functions (Slater-type, Gaussian-type, …) plain waves, wavelets, … This is one of the most important choices for every quantum chemistry code, usually, a single kind of functions is supported which is limiting the chemical systems that can be calculated with this. The most common distinction is made between codes that support periodic boundary conditions or not, while periodic boundary conditions are naturally included in plain wave and wavelet based programs, extra effort has to be put into codes using Gaussian-type basis functions to support this kind of calculation. Also, most wavefunction centered programs use atom centered orbitals since the resulting integrals are easier to solve. Here the exception from the rule is quite common in our field of research and usually offers a unique competitive edge.

For writing your Hartree-Fock program you do not have to bother with this choice, because we already made it for you by providing the implementation for integrals over Gaussian-type basis functions. We will limit you here to spherical basis functions only, so you can concentrate on coding the self-consistent field procedure and do not have to worry about mapping shells to orbitals to basis functions to primitives.

We will start with the input for the dihydrogen molecule in a minimal basis set. With the input format we provide, the geometrical structure of the system and the basis set are tied together.

```
12 2 2
20.0 0.0 -0.7 1.0 1
31.20
40.0 0.0 0.7 1.0 1
51.20
```

This input contains the information necessary to code up your program,
we start with the first line, which contains the number of atoms,
the number electrons and the total number of basis functions as integer
values.
None of this information is necessary to read the input file, but is
included for convenience, *e.g.* such that you can allocate memory before
starting to read the rest of the file.

Starting from the second line we expect a tuple of the three Cartesian coordinates in Bohr, the nuclear charges in multiples of electron charge and the number of basis functions this particular atom. In the lines after position and identity, we find the exponents of our basis functions, the number of lines following corresponds to the number of basis functions for this particular atom.

Important

You can compile your program using

```
make
...
gfortran ... -o ./scf
```

The resulting binary will be placed in the current working directory. Run it with

```
./scf
Here could start a Hartree-Fock calculation
```

The starting code provides the possibility to pass the input file as command line argument to your program with

```
./scf molecules/h2.in
Here could start a Hartree-Fock calculation
```

The input file will always be opened to the `input`

unit in `app/main.f90`

.

Important

The starting code is already setup as fpm project. The package manifest should contain the following content

```
name = "scf"
[[executable]]
name = "scf"
main = "prog.f90"
```

You can run your program with

```
fpm run
...
+ build/gfortran_debug/app/scf
Here could start a Hartree-Fock calculation
```

The starting code provides the possibility to pass the input file as command line argument to your program with

```
fpm run -- molecules/h2.in
...
+ build/gfortran_debug/app/scf "molecules/h2.in"
Here could start a Hartree-Fock calculation
```

The input file will always be opened to the `input`

unit in `app/main.f90`

.

Exercise 1

Before you start coding the input reader for this format, try to write with the specifications inputs for the following systems:

A helium atom in a double zeta basis set with exponents 2.5 and 1.0.

A helium-hydrogen cation with a bond distance of 2 Bohr in a minimal basis set on both atoms and an exponent of 1.25 for each atom.

Solutions 1

For a single atom the choice of the position is unimportant, so we can write something like

```
11 2 2
20.0 0.0 0.0 2.0 2
32.5
41.0
```

The helium-hydrogen cation looks similar to the dihydrogen, except for the changed nuclear charge.

```
12 2 2
20.0 0.0 -1.0 2.0 1
31.25
40.0 0.0 1.0 1.0 1
51.25
```

Exercise 2

Code an input reader that can read all the provided input files.

Make sure the input file is read correctly by printing all data read to the terminal.

### Classical Contributions#

First, we start by computing the classical nuclear repulsion energy, *i.e.*
the Coulomb energy between the nuclei.

Exercise 3

Which data is needed in the computation?

Code up the nuclear repulsion energy in a separate

`subroutine`

. Write the resulting energy in a meaningful way to the terminal.Evaluate the Coulomb-law for the dihydrogen molecule at 1.4 Bohr distance and compare it with your program. Can you run your

`subroutine`

multiple times and get the same result? If not recheck your code.Check what happens if you calculate the nuclear repulsion energy for a single atom. Do you get the expected result?

Classical contributions to the total energy do not dependent on the density or wavefunction and can already be calculated before starting with the self-consistent field procedure. It is usually a good idea to evaluate these contributions first, because they are less expensive to compute and, for the purpose of this course, easier to implement.

### Basis Set Setup#

This Hartree-Fock program will use contracted Gaussian-type orbitals to
expand the molecular orbitals in atomic orbitals.
We will use an STO-6G basis set, *i.e.* we use the best representation of a
Slater-type orbital by six primitive Gaussian-type orbitals.

This is the first time you will use an external library function, therefore
we will clarify to you how to read and use an `interface`

.
In your program, you will call a provided `subroutine`

to perform
the expansion from the Slater orbital to six primitive Gaussian-type orbitals.

The final call in your program might look somewhat similar to this:

```
call expand_slater(ng, zeta, exponents, coefficients)
```

To understand why the `subroutine expand_slater`

takes four arguments,
we look up its `interface`

:

```
interface
subroutine expand_slater(ng, zeta, exponents, coefficients)
import wp
integer, intent(in) :: ng !< number of primitive Gaussian functions
real(wp), intent(in) :: zeta !< exponent of slater function
real(wp), intent(out) :: exponents(:) !< of primitive Gaussian functions
real(wp), intent(out) :: coefficients(:) !< of primitive Gaussian functions
end subroutine expand_slater
end interface
```

An `interface`

provides the necessary information on how to
invoke a `subroutine`

or `function`

without concerning
you with the implementation details.

Note

for programmers coming from C or C++, it is similar to an `extern`

declaration in a header file for a function.

Usually, you do not have to write an `interface`

since they
are conveniently created and handled for you by your compiler.

Exercise 4

Create a

`parameter`

storing the information about the number of primitive Gaussian functions you want to expand your Slater function in.What is the optimal layout for saving the exponents and coefficients of the primitive Gaussian functions? Which quantities from the input do you need to determine the amount of memory to store them?

Allocate enough space to store all the primitive exponents and coefficients from the expansion for calculating the integrals later.

Loop over all basis functions, perform the expansion for each and save the resulting primitive Gaussians to the respective arrays.

### One-Electron Integrals#

Note that the basis set we have chosen is very simple, we only allow
spherical basis function (*s*-functions), also the contraction depth of
each function is the same. Usually one would choose more sophisticated
basis sets for quantitative calculations, but the basic principle remains
the same.

Integral calculations can quickly be very obscure depending on the way a basis set is stored and mainly handle implementation-specific details of the program to perform the integral evaluation in some clever way. We use a simple basis set here to teach you the basic principle of integral evaluation.

We start with the simple one-electron integrals, for Hartree-Fock we need
two-center overlap integrals, two-center kinetic energy integrals and
three-center nuclear attraction integrals.
To make things easier we provide the implementation for all three integrals
over contracted Gaussian orbitals, let’s check out the `interface`

:

```
interface
!> one electron integrals over spherical Gaussian functions
subroutine oneint(xyz, chrg, r_a, r_b, alp, bet, ca, cb, sab, tab, vab)
import wp
real(wp), intent(in) :: xyz(:, :) !< position of all atoms in atomic units, dim: [3, nat]
real(wp), intent(in) :: chrg(:) !< nuclear charges, dim: nat
real(wp), intent(in) :: r_a(:) !< aufpunkt of orbital a, dim: 3
real(wp), intent(in) :: r_b(:) !< aufpunkt of orbital b, dim: 3
real(wp), intent(in) :: alp(:) !< Gaussian exponents of the primitives at a
real(wp), intent(in) :: bet(:) !< Gaussian exponents of the primitives at b
real(wp), intent(in) :: ca(:) !< contraction coeffients of primitives at a
real(wp), intent(in) :: cb(:) !< contraction coeffients of primitives at b
real(wp), intent(out) :: sab !< overlap integral <a|b>
real(wp), intent(out) :: tab !< kinetic energy integral <a|T|b>
real(wp), intent(out) :: vab !< nuclear attraction integrals <a|Σ z/r|b>
end subroutine oneint
end interface
```

The most important information is we need *two* centers for the calculation,
meaning we have to implement it as a loop over all orbital pairs (= pairs of basis
functions).

Exercise 5

Which matrices can you compute from the one-electron integrals?

Allocate space for the necessary matrices.

Loop over all pairs and calculate all the matrix elements.

Exercise 6

Symmetric matrices, like the overlap, can be stored in two ways,
as full *N×N* matrix with `dimension(n,n)`

or as packed matrix
in a one-dimensional array with `dimension(n*(1+n)/2)`

, like:

Check if the matrices you have calculated are symmetric.

Pack all symmetric matrices.

The eigenvalue solver provided can only work with symmetric packed matrices, therefore you should have an unpacked to packed conversion routine ready for the next exercise.

### Symmetric Orthonormalizer#

The eigenvalue problem we have to solve is a general one given by

where **F** is the Fock matrix, **C** is the matrix of the eigenvectors,
**S** is the overlap matrix and **ε** is a diagonal matrix of the eigenvalues.
While this is in principle possible with a more elaborated solver routine,
we want to solve the following problem instead:

For this reason, we have to find a transformation from **F** to **F’**
and back from **C’** to **C**.
We choose the Loewdin orthonormalization or symmetric orthonormalization
by calculating **X** = **S**^{-1/2}, to invert (or perform any operation
with) a matrix we have to diagonalize it first and perform the operation on the
eigenvalues **s**.

Exercise 7

Use

**F’**=**X**^{T}**FX**and**C**=**XC’**to show that the general eigenvalue problem and the transformed one are equivalent.Write a

`subroutine`

to calculate the symmetric orthonormalizer**X**from the overlap matrix.Make sure that

**X**^{T}**SX**=**1**and that you are not overwriting your overlap matrix in the diagonalization.Do you see why

**X**is called symmetric orthonormalizer?

To solve the actual eigenvalue problem we will use the linear algebra package
(LAPACK), namely the subroutine `dspev`

, which stands for double precision,
symmetric packed eigenvalue problem, for more details you can look up the
documentation.
Since LAPACK routines can be somewhat unintuitive to work with on the first
encounter we provide a wrapper called `solve_spev`

.
As usual we checkout the interface in `src/linear_algebra.f90`

:

```
interface
subroutine solve_spev(matrix, eigval, eigvec, stat)
import wp
!> symmetric packed matrix to diagonalize, dim: n*(n+1)/2
! matrix is overwritten by the LAPACK solver
real(wp), intent(inout) :: matrix(:)
!> contains eigenvalues on exit, dim: n
real(wp), intent(out) :: eigval(:)
!> contains eigenvectors on exit, dim: [n, n]
real(wp), intent(out) :: eigvec(:, :)
!> error status, if not provided the routine will stop the program
integer, intent(out), optional :: stat
end subroutine solve_spev
end interface
```

Note

The linear algebra package (LAPACK) and the basic linear algebra subprograms (BLAS) are commonly used libraries in many quantum chemical programs, as they provide a computationally efficient and uniform interface to many linear algebra problems.

### Initial Guess#

There are now two possible ways to initialize your density matrix **P**.

Provide a guess for the orbitals

Provide a guess for the Hamiltonian

If you happen to have converged orbitals around, the first method would be
the most suitable.
Alternatively, you can provide a model Hamiltonian, usually a tight-binding
or extended Hückel Hamiltonian is used here.
The simplest possible model Hamiltonian we have readily available is the
core Hamiltonian **H**_{0} = **T** + **V**.

The initial density matrix **P** is obtained from the orbital coefficients by

where **n**_{occ} is a diagonal matrix with the occupation numbers of the
orbitals.
For the restricted Hartree-Fock case *n*_{*i*,occ} is two for occupied orbitals
and zero for virtual orbitals.

Exercise 8

By using

**F**=**H**_{0}as an initial guess for the Fock matrix we effectively set the density matrix**P**to zero. What does this mean from a physical point of view?Add the initial guess to your program.

Diagonalize the initial Fock matrix (use the symmetric orthonormalizer) to obtain a set of guess orbital coefficients

**C**.Calculate the initial density matrix

**P**resulting from those orbital coefficients.

### Hartree-Fock Energy#

The Hartree-Fock energy is given by

where Tr denotes the trace.

Exercise 9

You should have now an initial Fock matrix **F** and an initial density
matrix **P**.

Write a

`subroutine`

to calculate the Hartree-Fock energy.Calculate the initial Hartree-Fock energy.

### Two-Electron Integrals#

We have ignored the two-electron integrals for a while now, up to now they were not important, but we will now need to calculate them to perform a self-consistent field procedure. The four-center two-electron integrals are the most expensive quantity in every Hartree-Fock calculation and there exist many clever algorithms to avoid calculating them all together, we will again go the straight-forward way and calculate them the naive way.

Again we will check the `interface`

of the `twoint`

routine
in `src/integrals.f90`

.

Exercise 10

Allocate enough space to store the two-electron integrals.

Create a (nested) loop to calculate all two-electron integrals.

Exercise 11

The four-center integrals have an eight-fold symmetry relation we want to use when calculating such an expensive integral to speed up the calculation. Rewrite your loops to only calculate the unique integrals:

Write down all the indices for the four-center integrals and figure out unique relation between the indices (it is similar to packing matrices).

Implement this eight-fold-symmetry in your calculation.

You can also pack the two-electron integrals like you packed your matrices (you need to pack them three times, the bra and the ket separately, and then the packed bra and ket again). Note that this will make it later more difficult to access the values again, therefore it is optional.

### Self Consistent Field Procedure#

Now you have everything together to build the self-consistent field loop. Remember, the Fock matrix in terms of the two-electron integrals is given by

Make sure your definition of the density matrix matches the definition of the Fockian.

Exercise 12

First, you need to construct a new Fock matrix from the density matrix.

Construct a loop that performs your self-consistent field calculation.

Copy or move (whatever you find more appropriate) the necessary steps inside the SCF loop.

Define a convergence criterion using a conditional statement (

`if`

) based on the energy change and/or the change in the density matrix and`exit`

the SCF loop.Compare your final HF energies with

Input

E(RHF) / E

_{h}H

_{2}–1.127785613

He

–2.860251227

Be

–14.568567143

Calculate the proton affinity of H

_{2}. H_{3}^{+}has a D_{3h}structure with*R*_{HH}= 1.7 Bohr

### Dissociation Curves#

To calculate a dissociation curve we could either implement this functionality in the program itself or use an external script to automatically write input files and read the program output to collect the necessary data. Since this approach is common for computional chemistry, we will use a scripting approach here as well.

Modify the provided `bash`

script to calculate a dissociation curve:

```
#!/usr/bin/env bash
# stop on errors
set -eu
# Ensure non-localized printout
expot LC_NUMERIC=en_US.UTF-8
# put the name of your program here ("./scf" or "fpm run --")
program=echo
# unique pattern to find the final energy
pattern='final SCF energy'
# output file for plotting
datafile=plot.dat
# scan distances
start_distance=1.4
last_distance=5.0
step=0.1
read -r -d 'END' input <<-EOF
2 2 2
0.0 0.0 0.0 1.0 1
1.20
0.0 0.0 DIST 1.0 1
1.20
END
EOF
tmpinp=temporary.inp
tmpout=temporary.out
# cleanup
[ -f $datafile ] && rm -v $datafile
steps=$(seq $start_distance $step $last_distance | wc -l)
printf "Scanning from %.3f Bohr to %.3f Bohr in %d steps\n" \
$start_distance $last_distance $steps
for distance in $(seq $start_distance $step $last_distance | sed s/,/./)
do
# generate the input file
echo "$input" | sed s/DIST/$distance/ > $tmpinp
# perform the actual calculation on the input file
2>&1 $program $tmpinp > $tmpout
# get the energy from the program output
energy=$(grep "$pattern" $tmpout | awk '{printf "%f",$(NF)}' | tail -1)
# if there is no energy to be found, we complain
if [ -z "$energy" ]
then
1>&2 printf "ERROR!\n"
1>&2 printf "'%s' cannot be found in '%s' output\n" "$pattern" "$program"
1>&2 printf "please inspect '%s' and '%s'\n" "$tmpinp" "$tmpout"
exit 1
fi
# otherwise we write to the logfile
printf "Current energy is %.8f Hartree for distance %.3f Bohr\n" \
$energy $distance
printf "%8.3f %12.8f\n" $distance $energy >> $datafile
done
# cleanup
[ -f $tmpinp ] && rm $tmpinp
[ -f $tmpout ] && rm $tmpout
```

Exercise 13

The above script only contains dummies and can be executed without performing a calculation. Perform such a dry-run to understand how the script is working.

Modify it to match your program and plot the resulting dissociation curve.

Note

In case the script exits an error message like

```
script.bash: Zeile 34: printf: 1.4: Ungültige Zahl.
script.bash: Zeile 34: printf: 5.0: Ungültige Zahl.
Scanning from 1,000 Bohr to 5,000 Bohr in 37 steps
```

your language settings might be set to something other then English (in this case German). To overwrite the localisation settings run the script with

```
LC_NUMERIC=en_US.UTF-8 bash script.bash
```

## Properties#

With a working RHF program we have access to a wavefunction and we can use it to calculate certain properties.

### Partial Charges#

The total number of electrons in a system is given by

If your SCF wavefunction does not return its input number of electrons, the wavefunction is most certainly not properly normalized. In this case recheck your SCF and integral implementation.

Based on the above formula, Mulliken concluded that the number of electrons
associated with a particular nucleus is equal to the number of electrons
associated with its basis functions. Thus the partial Mulliken charge of atom *A*
is defined as:

Exercise 13

Code a subroutine to perform a Mulliken atomic population analysis

Determine the Mulliken atomic partial charges in LiH using the input file provided.

### Charge Density#

In terms of real contracted (Gaussian) basis functions, the electron density is given by

To calculate the product of the two basis functions the Gaussian product theorem can be used

where *K*_{AB} is given by

Note that the norming constants of the Gaussians have been absorbed into the contraction coefficients already.

Exercise 14

Code a subroutine to calculate the electron density

*ρ*at a given point in cartesian space.Calculate and plot the electron density along the axis connecting the bonds in H

_{2}and through the Be and He atoms. Plot them together.

## Geometry Optimization#

The next task is to expand your program to perform a simple geometry optimization.

### Numerical Derivatives#

The simplest way to perform geometry optimizations is by using the information from the energy derivative, since we do not want to code analytical derivatives of the Hartree–Fock energy expression, we will resort to numerical derivatives instead. This requires to evalulate several SCF energies in one program run, since you coded your SCF in a subroutine this should not be an issue. Nevertheless, try to run the subroutine several times to check if the code is correctly allocating and initializing its variables, they might show up now and you can fix them before adding much more code to your program.

Copy and modify your input files such that for each parameter (cartesian atomic
coordinates and Slater exponents) *θ*_{i} it contains an integer indicating
whether a numerical gradient with respect to *θ*_{i} is to be calculated.
Create a new subroutine that will allow the variation of one parameter
*θ*_{i} at a time as necessary for each element of the numerical gradient

Exercise 15

Code a numerical gradient that allows the calculation of the derivative of the HF energy with respect to atom positions and Slater exponents.

Save your gradient components to arrays analogous to the ones for the atom positions and Slater exponents. We shall call the conceptual combination of them the gradient vector

**g**.By definition, in which direction does a gradient point? How is this different from the force?

### Steepest Decent#

Similar to the SCF, parameter optimizations are performed iteratively and depend on a convergence criterion concerning energy and/or size of gradient.

Code a variant of the “steepest descent” optimization routine as given by:

*k* denotes the number of the optimization cycle, *Θ* is the parameter set for an
iteration. Choose *η* to get smooth, fast convergence.

Exercise 16

Optimize the geometry of HeH

^{+}in a minimal basis set.Optimize the Slater exponents of Be and H

_{2}(*R*_{HH}= 1.4 Bohr) in a full double-ζ basis set. Compare to energies for the est. HF basis set limit:System

Energy/E

_{h}H

_{2}–1.134

Be

–14.573

## Unrestricted Hartree-Fock#

Copy and modify your RHF subroutine to create an UHF program. The input files will
need to be modified such that they contain the number of *α* and the number of
*β* electrons. By definition, the number of *α* electrons is bigger then the number
of *β* electrons. You can check the correctness of your results against the Li atom
(Slater exponents 3.5, 2.0, 0.7 and 0.3: E = –7.419629 E_{h}) and the H
atom.

For UHF you use one set of **F**, **C** and **P** matrices for each spin to solve
the two eigenvalue problems

concurrently. They are coupled only through the formation of the Fock matrix
*via* the Coulomb interaction as demonstrated for **F**^{α} below:

Note that the calculation of **P** differs in the occupation number and for
closed-shell test cases provided your UHF program must give RHF results.
For this reason, you have to break the spatial symmetry of your system to obtain
the UHF solution, if it is available. The easiest way to do this is through the
initial guess. If you use a symmetric guess (like **P** = 0), you will only find
the RHF solution.

The UHF Energy is given by:

Exercise 17

Calculate and plot the dissociation/potential curves for

^{3}H_{2},^{1}H_{2}and^{1}Li_{2}.The exchange reaction H

_{2}+ H → H + H_{2}has a linear symmetric transition state. Find it and its relative energy.Calculate the first ionization potential of Be. The experimental value is 9.3 eV.

Are He

_{2}^{+}or^{3}He_{2}bonded?

### Spin Contamination#

Recall that spin contamination can occur in UHF calculations, *i.e.* deviations
of the expectation value of the square of the total spin angular momentum operator
S^{2} from the ideal value:

Calculate the spin contamination according to

where *i*, *j* indicate *α* and *β*-MOs, respectively.

Exercise 18

For the system:

h 0.0 0.0 -1.0 h 0.0 0.0 0.0 h 0.0 0.2 1.0

and a Slater exponent of 1.24, you should find a spin contamination of 0.004682 and a UHF energy of -1.265643.

Plot the spin contamination along both RHF and UHF dissociation curves of H

_{2}.Why is the spin contamination of RHF always 0?

## Møller–Plesset Perturbation Theory#

Using your RHF program, code a subroutine to calculate the RMP2 energy after
your RHF procedure has converged. The RMP2 correction to the HF energy
(assuming real MOs *a, b, r, s*) can be written as

As for all correlation methods, you will need to have two-electron integrals
over MOs (denoted by Roman letters) instead of over AOs (Greek letters).
The process to obtain them is called AO-MO-Transformation.
There are different algorithms possible, among them are one to scale
with *M*^{8} and one with *M*^{5} that takes approximately twice
as much memory.
For now, code up the straigthforward algorithm to transform the two-electron
integrals. This is the one that scales with *M*^{8}.

Input |
E(RHF) / E |
E |
---|---|---|

H |
–1.127785613 |
–0.012541746 |

He |
–2.860251227 |
–0.012686549 |

Be |
–14.568567143 |
–0.014939565 |

In contrast to the *M*^{8} variant, where you transform from
\(\left({\mu \nu}|{\lambda \kappa}\right)\) to
\(\left({ar}|{bs}\right)\) at once, you need to transform the
two-electron integrals step-by-step:

Exercise 19

Copy your RMP2 routine and modify it such that your integral transformation scales with

*M*^{5}.For the provided examples H

_{2}and LiH compare the efficiency of both algorithms by counting the number of cycles passed through in the most inner loops. What happens if you increase the number of Slater functions by a factor of two?Calculate the dissociation curve of H

_{2}with RMP2. Plot both the RHF and RMP2 curves on a relative energy scale (in kcal/mol) in a minimal STO basis with Slater exponent of*ζ = 1.2*. As in the provided input for H_{2}, the Hartree-Fock energy of a hydrogen atom is –0.4798356 E_{h}. What behavior do you observe compared to RHF close to the equilibrium and at far distances*R*= 10 Bohr?Calculate the dissociation curve of He

_{2}between 4 and 10 Bohr with RHF and RMP2 and plot them on a relative energy scale (in kcal/mol). You will find that there is a minimum in each curve. From your knowledge about HF and MP2, did you expect this behavior? What effect could be the cause for the minimum in the RHF curve?