## A Simple Introduction to Computational Quantum Mechanics

This piece is intended as a gentle introduction to basic concepts in computational quantum mechanics. It's meant for a reader with a basic understanding of linear algebra and derivatives, or at least matrices. If you wish to follow along with the code, then you also need Python 3 and the popular scientific library numpy. We will restrict ourselves to the consideration of one-dimensional systems. I have chosen here to explain rather extensively the concepts used in the (very short) code, so that the reader may apply similar principles also to other problems.

First, let us recall to mind some basic features of (non-relativistic) quantum mechanics.
The basic equation to be solved is the celebrated *Schrödinger equation*,

Here, the letter H stands for the Hamiltonian ("the energy operator"), and the letter E for the energy of the solution.
The Schrödinger equation is an *eigenvalue equation*: the solution \( \psi \) is that state which is an eigenvector of
the Hamiltonian \( \hat{H} \) with eigenvalue \( E \).

In quantum mechanics, operators like \( \hat{H} \) correspond to things that can be observed. Thus, what this equation says
is that, when you operate on some state with the Hamiltonian ("measure the energy"), it doesn't change (except by a constant factor).
Thought the Schrödinger equation can be motivated by wave mechanics, it is usually taken as a *postulate*, accepted on
experimental grounds without further theoretical justification.

I called \( \hat{H} \) an energy operator. Well, what is in an "energy operator"? The answer is, for our purposes, that it contains similar elements as the classical variant: it has a term representing kinetic energy (though with a minus sign), and a term representing potential. In classical mechanics, we would write this \( H=\frac{1}{2} m v^2 + V(x) \). As an operator, we write

\begin{align} \hat{H} = -\frac{\hbar ^2 }{2m}\frac{\mathrm{d}^2}{\mathrm{d}x^2} + v(x) \end{align}The only difference to the classical case is really the kinetic energy term. For the purposes of following along, it's not really necessary to understand the following explanation, but I'll give it for completeness. First, we note that \( \frac{1}{2}mv^2 = \frac{p^2}{2m} \) where \( p = mv \). Now using the quantum-mechanical momentum operator \( \hat{p} = -i\hbar \frac{\partial}{\partial x} \) we immediately arrive to the quantum mechanical form of the kinetic energy, as now \( \hat{p}^2 = -\hbar ^2\frac{d^2}{dx^2} \).

We are faced with a rather natural question: what does it mean for something to be an eigenvector of a derivative operator?
If you've been taught about matrices, it certainly doesn't seem like a derivative should have anything to do with a matrix. Indeed, it is not
possible to represent a derivative operator such as the kinetic energy as a matrix (you would end up with an infinite matrix), and the mathematical
theory that justifies talk of "eigenvalues" and "eigenvectors" in this context is rather complicated. These sorts of objects are called infinite-dimensional
operators. However, computers can't deal with infinite numbers of points anyways. In order to make progress, we make our first approximation: we take
our space to consist of a *finite number of intervals* of length \( \Delta x \). That is to say, if the length of our one-dimensional system was 10 units,
then we would divide those 10 units to (for example) lengths of \(\Delta x= 0.01\) units.

How does this help us? Well, recall how derivatives are defined:

\begin{align} \frac{d f(x)}{dx} = \lim _{h\rightarrow 0}\frac{f(x+h)-f(x)}{h} \end{align}In our space consisting of finite intervals, instead of the h going to 0, it will instead go to our minimum length, \( \Delta x \). It is as though we have replaced the continuous universe with our own universe, in which we can't speak of length scales smaller than \( \Delta x\).

Now, how do we cast this derivative in to the form of a matrix? First, in our universe, also the wave function \( \psi \) will be discrete, being defined on the finite intervals, so that we may write it as a column vector

\begin{align} \psi = \begin{bmatrix} \psi (x_1) \\ \psi (x_2) \\ \vdots \\ \psi (x_f) \end{bmatrix} \end{align}If we were to cast the first derivative in matrix form, it would need to be such a matrix (call it A) that when we operate to this column vector, we get the result

\begin{align} A\psi = \frac{1}{\Delta x}\begin{bmatrix} \psi (x_2)-\psi (x_1) \\ \psi (x_3)-\psi (x_2) \\ \psi (x_4)-\psi (x_3) \\ \vdots \\ \end{bmatrix} \end{align}since in our universe \( \psi(x_i+h) = \psi (x_i+\Delta x) = \psi (x_{i+1} )\) and \( h = \Delta x\). To find such a matrix, we need to remember the multiplication rule for matrices, which is (for matrix A and column vector \(\psi \): \( (A\psi)_i = \sum _{j=1}^n A_{ij}\psi _j \). So for the first element, we need to have \( (Av)_1 = \sum _{j=1}^n A_{1j}\psi_j =\frac{\psi (x_2) - \psi (x_1) }{\Delta x} \). This means that \( A_{11} = -1/\Delta x \) and \( A_{12} = 1/\Delta x \) and the rest of the elements on that row must be zero, as then

\begin{align} \sum _{j=1}^n A_{1j}\psi_j &= A_{11}\psi _1 + A_{12}\psi _2\\ &= A_{11}\psi(x_1 ) + A_{12}\psi (x_2 ) \\ &= \frac{1}{\Delta x} \bigg( \psi (x_2 ) - \psi (x_1 ) \bigg) \end{align}Using similar reasoning, we deduce also the rest of the elements, which creates a matrix like

\begin{align} A = \frac{1}{\Delta x}\begin{bmatrix} -1 & 1 & 0 &\cdots \\ 0 & -1 & 1 & 0 & \cdots \\ 0 & 0 & -1 & 1 & 0 & \cdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots \end{bmatrix} \end{align}However, we need the second derivative, not the first. This is easily fixed, since we can just apply the derivative operator to itself to get the second
derivative. Notice, however, that we defined the derivative as the "forward difference", where we took \( f(x+h)-f(x) \) instead of \(f(x)-f(x-h) \). In a continuum
case for a continuous function this wouldn't matter, but in our discrete space it might. To balance things out, we apply the derivative operator
to the (negative of the)* transpose* of itself, thus multiplying forwards- and backwards differences (convince yourself that the negation of the transpose is the backwards difference!)
Now we arrive at the matrix

This is a matrix with the diagonal and the elements just above and below the diagonal filled, with the rest zero. Great! Now all we need to do is multiply this with \( -\frac{\hbar ^2 }{2m} \) and we have our kinetic energy operator.

How big is this matrix exactly? If we choose our system to have 1000 points, then it must apply to every point of \( \psi \), that is, to all values \( \psi (x_1) ... \psi (x_{1000} ) \). It is then easy to see the matrix must be 1000 by 1000! That's a rather big matrix. Although it is not a problem for modern computers, we can see that our naive method doesn't scale very well at all.

What of the potential \( v(x) \)? Suppose that we have the harmonic potential, \( x^2 \). How do we represent it? Fortunately, this is much easier than the derivative. The expression \( x^2 \psi (x) \) just means we multiply \( \psi (x) \) with \( x^2 \) at every point, so the matrix is just diagonal! Let us forget for now, however, about he potential and set it o zero. We'll first investigate just the eigenvectors of the kinetic energy operator.

This would ordinarily mean we're dealing with a "free particle" that doesn't interact and is not affected by
external forces. However, our method has a drawback: it obeys *boundary conditions* which set \(\psi(x) \)
to zero at the edges of our 1-dimensional world (hint: look at the first and last elements, see the supplementary page if you want details). Thus, if we solve the
eigenvalues of our kinetic energy operator on the computer, we should actually get the result of
particle in a box . This is not a serious drawback: usually
we want to study a limited system anyway, and we could make our system quite large, put the interesting part in the middle, and
then our infinite walls at the beginning and end wouldn't matter much.

Now it's time to get to the code. We import numpy and create our universe (the x variable):

```
import numpy as np
```

#x goes from 0 to 4 with 1000 points linearly spaced

x = np.linspace(0,4,1000)

Creating the kinetic energy operator. We now move to units where \( \hbar = 1 \) and \( m = 1\), the so-called "atomic units":

` ````
h = x[1]-x[0]
```

kinetic_energy = 1/h**2 * 1/2 *(2*np.eye(len(x)) + np.diagflat(-np.ones(len(x)-1),1) + np.diagflat(-np.ones(len(x)-1),-1))

And now we can find and plot the energy eigenvalues/vectors:

```
energies,vec = np.linalg.eigh(kinetic_energy)
```

print(energies[0:4])

The printing should give you the first 4 values as 0.307, 1.228, 2.76, and 4.91, which agrees quite well with an exact solution of the particle in a box problem! If you have matplotlib installed, you can do the following to see the wave functions:

```
import matplotlib.pyplot as plt
```

plt.plot(x, vec[:,0]) #To see the lowest energy wave function; the next one is vec[:,1], etc

It should look like this:

Now we might like to return to the problem of an external potential affecting the system. As I said, those are diagonal: for the harmonic potential, we would get

```
harmonic = np.diagflat(x**2)
```

H = kinetic_energy + harmonic

which we now would sum to the kinetic energy operator. I leave it as an exercise to the reader to see the results of this potential ;). You might also want to explore other potentials. What happens when the potential is constant for, say, first half of the box? (Add the potential only to the first 500 elements of the diagonal). What if there's a constant potential in the middle but not at edges? (Add the constant to the middle of the diagonal, say, elements 400-600). You can explore an endless number of different ideas with this same method.

That's all for now. Click here for the index of this site, where I may write more articles in the future. (There's also a link to the RSS feed and a twitter account).