1s Program

Quantum Chemistry

What is a 1s program?

Well, this is a terminology that I developed myself to describe a quantum chemistry program which only uses 1s slater orbitals.

Why 1s Slater Orbitals?

In this case we have restricted ourselves to 1s slater orbitals for simplicity. The 1s functions have the unique property that they have the easiest integrals which allows us to focus on the logic of the program as opposed to the implementation of some crazy integrals. I would say that an apt comparison would be implementing floating point addition in a numerical methods class. While it could be done you would probably spend one semester debating the IEE 754 specification instead of actually calculating anything useful.

There are so many moving parts in a quantum chemistry which require a deep understanding and slowing the presentation down by presenting advanced integral calculations does the learner a disservice. If you've ever watched a movie with too many flashbacks you know exactly what I'm talking about. Something interesting is about to happen -- oh wait 2 episode flashback that was not important to the story.

The logic of the Hartree-Fock program is hard enough without the 100+ pages on molecular integrals in Helgaker's book. I believe this is why Szabo and Ostlund, as well as Thjissen implement 1s programs. The point is that this program is very limited in scope on purpose. If one wanted to extend the program to do more complicated calculations you would just need to change how the matrix elements are calculated.

Linear Combination of Atomic Orbitals

This is exactly what it sounds like. You take an atomic orbital centered on each atom, and represent the system using this basis. The representative orbitals are usually exponentially decaying functions called slater orbitals.

An example

If we take Slater orbitals centered on the point $T$ and the point $U$ . The reason for not using $A$ and $B$ is that these are the generic labels used in the matrix element formulas.

ψ_{T}^{STO} (r) = exp (- ζ ∣ r - R_{T} ∣)

ψ_{U}^{STO} (r) = exp (- ζ ∣ r - R_{U} ∣)

STO-nG

Now that we have established we want to represent our wavefunction as a linear combination of slater orbitals we might ask why we did all this work with Gaussians. Well, the integrals are very easy for Gaussians. We will represent each slater orbital with a set of gaussians. We now have a linear combination of atomic orbitals which are a linear combination of Gaussians. Just to be sure we understand we would have 2 total basis functions with slater-type orbitals, but with GTOs we have for example 3 for each slater, for 3 times 2 basis functions. This leaves us with a larger matrix, but the matrix calculation is much easier with gaussians.

An example

If we take Slater orbitals centered on the point $T$ and the point $U$ . The reason for not using $A$ and $B$ is that these are the generic labels used in the matrix element formulas.

ψ_{T}^{STO} (r) = exp (- ζ ∣ r - R_{T} ∣)

ψ_{U}^{STO} (r) = exp (- ζ ∣ r - R_{U} ∣)

Now we can approximate each Slater orbital as a sum of gaussians. For an arbitrary orbital centered at the point $Z$ we have.

ψ_{Z}^{STO} (r) = exp (- ζ ∣ r - R_{Z} ∣) \approx ψ_{Z}^{GTO} (r) = i = 1 \sum n c_{i} exp (- ρ_{i} ∣ r - R_{i} ∣^{2})

This is where the STO-nG nomenclature shows up because we approximate an STO with $n$ Gaussians. We apply the gaussian approximation to $ψ_{A}^{STO}$ and $ψ_{B}^{STO}$ , yielding

ψ_{T}^{STO} (r) \approx ψ_{T}^{GTO} (r) = c_{1} exp (- ρ_{1} ∣ r - R_{1} ∣^{2}) + c_{2} exp (- ρ_{2} ∣ r - R_{2} ∣^{2}) + c_{3} exp (- ρ_{3} ∣ r - R_{3} ∣^{2})

ψ_{U}^{STO} (r) \approx ψ_{U}^{GTO} (r) = c_{4} exp (- ρ_{4} ∣ r - R_{4} ∣^{2}) + c_{5} exp (- ρ_{5} ∣ r - R_{5} ∣^{2}) + c_{6} exp (- ρ_{6} ∣ r - R_{6} ∣^{2})

For a molecule like hydrogen we need both wavefunctions to form a linear combination of atomic orbitals. This means we need to sum two sets of three GTOs. This is why the index runs from $1$ to $6$ . We have two STO-3G approximations. Indices $1, 2, 3$ correspond to $A$ and indices $4, 5, 6$ correspond to $B$ .

How to Construct an STO-nG basis set

Turns out that it can be done using a nonlinear least-squares procedure.

Θ (ψ_{A B}) = [[ρ_{1}, R_{1}], [ρ_{2}, R_{2}], [ρ_{3}, R_{3}], [ρ_{4}, R_{4}], [ρ_{5}, R_{5}], [ρ_{6}, R_{6}]]

Forming the matrix

H_{ij} = ⟨ ψ_{i} ∣ H ∣ ψ_{j} ⟩

S_{ij} = ⟨ ψ_{i} ∣ ψ_{j} ⟩

H C = E SC

For the detailed formulas and Python implementations of all 1s Gaussian-type orbital integrals, see the 1s Integrals page.

Putting it all together

We now need to form our matrix elements using the above formulas. Our Hamiltonian will look like this:

H = V_{ee} + V_{n e} + T

All we are missing now are our basis functions. Our model system will be h2+ and then h2

Practical considerations

Each matrix element is a function of six parameters, three for each basis function.

2 contraction coefficients
2 gaussian exponents
2 gaussian centers

One key observation we can make is that the contraction coefficients do not factor into the integrals. They can be multiplied at the end using a hadammard product.

We introduce the function $Θ$ to denote the extraction of the parameters from the basis function, $ψ$ . $α, R_{A}$ will denote the parameters for $ψ_{i}$ :

Θ (ψ_{i}) = [α, R_{A}]

and $β, R_{B}$ for $ψ_{j}$ :

Θ (ψ_{j}) = [β, R_{B}]

This is just a helpful shorthand to make it so I don't have to completely change the notation. If you look above most of the derivations make use of the variables $α, R_{A}, β, R_{B}$ which refer to two different generic gaussians. I could have labeled them by some other indexed variables like

ρ_{1}, R_{1}, ρ_{2}, R_{2}

but I am saving this convention to index the gaussian which we will see later. The upshot is that the integrals only depend on certain parameters of the gaussians, you don't have to explicitly form the basis functions. From a programming perspective, this function $Θ$ would be like a getter method. So when we are practically calculating an arbitrary matrix element

M_{ij} = F ([α, A], [β, B])

The math can look a little closer to the code. The function $F$ will depend on which part of the hamiltonian we form.

A note on the matrix element formula

All too often in quantum mechanical calculations we run across the operation of calculating matrix elements. I just wanted to take a second to think about what this operation means in terms of other operations. This notation for an integral is foreign to the non quantum mechanic and there is no native support for this operation in programming languages. I guess the point I'm trying to make is that you don't get to use this neat bra-ket notation when you program it so it's not as natural to think about as programming $f (x) = sin (x)$ .

We have a list of functions and then must perform some kind of integral outer product in order to form a matrix. In mathematical terms we have a two-argument functional. A functional in the context of quantum chemistry or density functional theory means a function which takes a function and outputs a number. A very basic example is shown below. This functional takes a function and integrates it from zero to one.

f [f (x)] = \int_{0}^{1} f (x) d x

A concrete example:

f [x^{2}] = \int_{0}^{1} x^{2} d x = \frac{1}{3}

Another example:

f [x^{3}] = \int_{0}^{1} x^{3} d x = \frac{1}{4}

For the case of matrix elements, we take two functions as arguments, insert those functions into a functional and from there we get a value. Something that looks like this:

F [ψ_{i}, ψ_{j}] = F_{ij} = ⟨ ψ_{i} ∣ \hat{O} ∣ ψ_{j} ⟩ = \int_{- \infty}^{\infty} ψ_{i}^{*} \hat{O} ψ_{j} d x

This functional $F$ is the combination of the usual inner product and the operator. Honestly, this just boils down to personal preference, but I find it much less jarring to think about the matrix elements as functionals because the notation looks closer to what is programmed. Additionally, if we program this method in any programming logic the integral will be inside some function.

def get_matrix_element_for_operator_o(psi_i, psi_j):
    op_psi_j = operator(psi_j)
    o_ij = do_integral(psi_i, op_psi_j)
    return o_ij

You could do this all without explicitly splitting the operator and integral portion inside the function. The only problem with this is that there are different strategies or methods for realizing the operator and the integral. To me there are two main possibilities:

There is an analytical form in which you can just program the matrix element as a function of the basis functions (1s GTO)
You evaluate each basis function on a grid then numerically evaluate operators and integrals (think finite differences, trapezoidal rule, simpson's method)

M_{ij} = F ([α, A], [β, B])

⟨ ψ_{i} ∣ \hat{O} ∣ ψ_{j} ⟩

The Fully Numerical Approach

One fascinating property of the trapezoidal rule is that it is exponentially convergent for gaussian-type functions.

Check your understanding

Hopefully after reading this article you should be able to answer a few questions:

Questions

What is a generalized eigenvalue problem?
Why do we use 1s basis functions?
What is a bond order matrix?
How do we fit an STO-nG basis?
Describe the procedure for calculating matrix elements.
Describe two different ways to think about calculating matrix elements?
What parameters do the matrix elements in a gaussian type orbital program depend on?
Why can we solve the helium atom using the Hartree-Fock method? What is special about it?
What is the role of the zeta in the slater type orbital?
If I have an STO-6G single zeta basis for an h2+ atom how many basis functions are there?
What is meant by linear combination of atomic orbitals?
What is the variational principle?
What integrals are difficult to calculate using slater orbitals?
What methods are used to calculate matrix elements for basis functions with higher angular momentum states?