\documentstyle{article}
\setlength{\oddsidemargin}{0in}
\setlength{\topmargin}{-0.5in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{9in}
\begin{document}

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}

\begin{center}
{\bf Hartree-Fock: One  Case Study of  Numerical Linear Algebra
Interacting with Computational Chemistry}\\
{\bf Alan Edelman } \\
UC Berkeley \\
May, 1993
\end{center}

Scientific Computing, now an academic discipline of its own at
major universities, may be defined as the specialty that
solves those computational problems on the thin boundary
that separates the infeasible from the trivial.
Success is the ability to solve a problem once before thought
impossible.
The rocky road to
progress  also includes problems that cross the border in 
the opposite direction.  New computer architectures 
provide the bumps in the road. 

One characteristic of this maturing science is the emphasis on
interdisciplinary teams.  
It is a sign of the times that our lives have become 
increasingly more interdependent in so many ways.
The scientist (chemist, physicist, or engineer, say)
needs the computer specialist to efficiently use
a machine architecture. The well-equipped
numerical analyst can serve in the middle.  
One way he may do this is by suggesting 
appropriate algorithms.  More significantly, however, 
the well-equipped numerical analyst is 
bilingual enough to comfortably discuss continuous
problems and also modern architectures. 

What are the road blocks to interdisciplinary communication?
The most obvious hindrance is jargon.  However, a more
subtle hindrance is becoming increasingly clear.
Scientist A  thinks in physical terms; computer specialist B
is in the trenches of architecture or
software.  Pay attention next time to the "A Says--B  Says"
of a conversation between a physical scientist and a computer
specialist.    

In the next section, every effort has been made to explain
the important aspects of Hartree-Fock Chemistry Computations
from a numerical analysts perspective.
Our treatment differs from those of many excellent chemistry
references in that we discuss the computational aspects
from the start and include a toy problem to illustrate
the ideas. 

\section{Hartree-Fock Computations: A Simple Example} 

From a linear algebra point of view, Hartree-Fock computations
repeat an inner loop consisting of
\begin{enumerate}
\item The solution of a linear system of equations corresponding
to a Poisson equation and
\item The solution of an eigenvalue problem corresponding to
an approximation to Schr\"{o}dinger's equation.
\end{enumerate} 


In this ``chicken and egg'' computation, the 
solution to the Poisson equation is an input to the
eigenvalue problem, and the solution to the eigenvalue
problem is an input to the Poisson equation.
The hope is that the process will converge, and 
when we are sufficiently satisfied that we have done
so, we say that we have reached self-consistency.


The simplest example would be the lowst energy state
of the Helium atom.   One reason that it is simple
is that by symmetry, the problem becomes one dimensional.
The Poisson equation is one dimensional with unknown
$u(r)$ and known function $n(r)$:
\be
\frac{d^2}{dr^2} u(r)=-4\pi r n(r).
\ee
The eigenvalue problem has unknown smallest eigenvalue
$E$ and corresponding eigenfunction $\theta(r)$:
\be
\left[  -\frac{1}{2}\frac{d^2}{dr^2}-\frac{2}{r}+\frac{u(r)}{r}
 \right] \theta(r)
 =E \theta(r),
\ee
where
$$
4\pi \int_{r=0}^\infty \theta(r)^2dr =1, \
\theta(r)=r n(r)^{1/2},
u(0)=0,
{\mbox{ and }} 
 u(\infty)=1.
$$   
The above equations specify the normalization of the
eigenfunction, the connection between the eigenvalue
problem and the Poisson equation,
and the boundary conditions for the Poisson equation
respectively.
Since the  function $u(r)$ depends on the eigenfunction
$\theta(r)$, we have a non-linear eigenvalue problem.

The most straightforward discretization of this problem
would break the real line into equal intervals of
length \verb+d+ from $r=0$ to \verb+rmax+.

A small Matlab program that does this is

\begin{verbatim}
% input d and rmax
r=[d:d:rmax]'; nr=r.*exp(-r); k=length(r);    %Discretize r,intitial guess n(r)
m=-eye(k)+diag(ones(1:k-1),1);m=(m+m')/d^2;   %Second Order difference matrix

while(1)
   nr=nr/(4*pi*d*sum(r.^2 .* nr));            % Normalize nr
   b=-4*pi*r.*nr;b(k)=b(k)-1/d^2;             % right hand side
   u=m\b;                                     % solve Poisson equation
   h=diag((u-2)./r)-m/2;                      % Define Hamiltonian
   e=min(eig(h)),th=(h-e*eye(k))\((1:k)'==k); % Find smallest eig pair
   nr=(th./r).^2;                             % next n(r)
end
\end{verbatim}

The boundary condition $u(\infty)=1$ is built into
the right hand side of the Poisson equation.
We chose to perform inverse iteration for the eigenvector
corresponding to the smallest eigenvalue.  


\section{Scientific Background}


The Hartree-Fock approximation for the lowest
energy of a  system with $k$ electrons
replaces
the linear  Schr\"{o}dinger equation which involves
a function $\psi$ of
$3k$ spatial coordinates (three for each of the two
%electrons) with an approximation that requires
only three coordinates.  The cost of this approximation
is that the linear Schr\"{o}dinger equation
becomes non-linear  and also the solution to the
problem is now merely an approximation.


Let $x_1,x_2\in R^3$ denote the positions of the
two electrons in a Helium atom.  Ideally, we would
like to obtain $\psi(x_1,x_2)$ the wave function
for the electrons, but we settle for an approximation 
$\psi_1(x_1)$ and $\psi_2(x_2)$.  The function
$\psi_1(x_1)$ is obtained by assuming $\psi_1(x)=\psi_2(x)$
(there is no difference between the two electrons) and
that the potential that the second electron ``feels''
is the average potential obtained from the first electron.




To be precise, the average  potential energy $U(x)$ due to the
first electron is the solution to
\be
\nabla^2 U(x) = -4\pi |\psi_1(x)|^2.
\ee
This is exactly the potential that would occur if 
space were filled with a charge density of $|\psi_1(x)|^2$.
If we instead, think of $|\psi_1(x)|^2$ as a probability
density ($\int |\psi_1(x)|^2 dx=1$), for the location of
one electron, we then are computing the average potential.
This is a Poisson equation, whose numerical discretization
is a linear system of equations.

An electron (labeled 2)
that is attracted to a nucleus by a fixed
Coulomb potential of $2/r$ (two protons) and is
under the influence of  a potential $U(x)$ has a wave function
and energy given by
\be
\left[ -\frac{1}{2}\nabla^2- 2/\|x\|+U(x) \right] \psi_2(x) = E \psi_2(x).
\ee
When discretized this is a matrix eigenvalue problem.

Given a $\psi_1(x)$ we can compute a $U(x)$ (Equation (1)) from which
we can compute a $\psi_2(x)$ (Equation (2)). If we use $\psi_2(x)$ as
our new $\psi_1(x)$ we can repeat the procedure iterating
hoping to converge to the symmetric situation $\psi_1(x)=\psi_2(x)$.

An extra simplicity that makes this problem suitable for
discussion in elementary numerical linear algebra classes occurs
when we are interested in the lowest energy state.
In this case, the wave functions and potentials are all
spherically symmetric, and thus only depends on the
radial parameter $r$.  The equations in terms of $r$
were given in the previous section.


\section{Alternative Formulation}
The problem could be stated as an optimization problem:
minimize an energy expression $E(x)$ subject to $\|x\|=1$.
If it happened to be the case that $E(x)=x^TAx$ for some
symmetric linear operator $A$, then we would say that
we have an eigenvalue problem.  Our problem 
is a non-linear eigenvalue problem.


\end{document}


