\documentclass[10pt]{article}
\usepackage{a4wide}
\usepackage{listings}
\begin{document}

\section*{Format for delivery of report and programs}
%
The format of the project is that of a printed file or hand-written report.
The programs should also be included with the report. Write only your candidate number on the first
page of the report and state clearly that this is your report for project 5 of FYS3150, fall 2011. 
There will be a box marked 'FYS3150' at the reception of the Department of Physics (room FV128). 

\section*{Project 5, Variational Monte Carlo studies of light atoms, deadline  December 12, 3pm}
For this project you can build upon program programs/chapter14/program1.cpp (or the f90 version).
You will need to parallelize exercises b-d and you should therefore use parts of project 3. 
You can also use the corresponding python code under the link for project 5. This has been written by Brynjar 
Arnfinsson.

The aim of this project is to investigate the variational Monte Carlo method applied
to light atoms such as helium and beryllium.
Various trial wave functions are to be tested and compared.
The aim is to find wave functions which reproduce the best possible theoretical ones (which are close to the experiment) 
energies as best as possible. We will in the text call these energies 'experimental ones', since they are very close 
to data. 
  
\subsection*{Ground state energy of helium}

Helium consists of two electrons and a nucleus with charge
$Z=2$. We are going to use the Born-Oppenheimer approximation in modelling the system, 
assuming thereby that we can neglect nucleonic degrees of freedom. 
The nucleus, whose extension is on the order of $\sim 10^{-15}$ m
(roughly six order of magnitude smaller than interatomic distances), is taken to be a point charge with mass
much larger than that of the electrons.
Electrostatic forces constitute then the essential contribution to the potential energy,
given in this case by the attraction experienced by every electron from the nucleus
and the repulsion between the two electrons.  

We label $r_1$ the distance from electron 1 to the nucleus and similarly 
$r_2$ the distance between electron 2 and the nucleus.
The contribution to the potential energy from the interactions between the 
electrons and the nucleus is
\begin{equation}
   -\frac{2}{r_1}-\frac{2}{r_2},
\end{equation} 
and if we add the electron-electron repulsion with
$r_{12}=|{\bf r}_1-{\bf r}_2|$, the total potential energy 
$V(r_1, r_2)$ is
\begin{equation}
 V(r_1, r_2)=-\frac{2}{r_1}-\frac{2}{r_2}+
               \frac{1}{r_{12}},
\end{equation}
yielding the total Hamiltonian
\begin{equation}
   \hat{H}=-\frac{\nabla_1^2}{2}-\frac{\nabla_2^2}{2}
          -\frac{2}{r_1}-\frac{2}{r_2}+
               \frac{1}{r_{12}},
\end{equation}
and Schr\"odinger's equation reads
\begin{equation}
   \hat{H}\psi=E\psi.
\end{equation}
All equations are in so-called atomic units. The distances
$r_i$ and $r_{12}$ are dimensionless. To have energies in electronvolt
you need to multiply all results with 
$2\times E_0$,
where $E_0=13.6$ eV.
The 'experimental' energy for the ground state of helium in atomic units is 
a.u. is $E_{\mathrm{He}}=-2.9037$ a.u..
The basic wave functions we will employ in this exercise are 
\begin{equation}
   \psi_{T1}({\bf r_1},{\bf r_2}, {\bf r_{12}}) = 
   \exp{\left(-\alpha(r_1+r_2)\right)}, 
\end{equation}
and
\begin{equation}
   \psi_{T2}({\bf r_1},{\bf r_2}, {\bf r_{12}}) = 
   \exp{\left(-\alpha(r_1+r_2)\right)}
   \exp{\left(\frac{r_{12}}{2(1+\beta r_{12})}\right)}, 
\end{equation}
with  $\beta$ as a new  variational parameter.
Your task is to perform a Variational Monte Carlo calculation
using the Metropolis algorithm to compute the integral
\begin{equation}
   \langle H \rangle =
   \frac{\int d{\bf R}\psi^{\ast}_T({\bf R})H({\bf R})\psi_T({\bf R})}
        {\int d{\bf R}\psi^{\ast}_T({\bf R})\psi_T({\bf R})},
\end{equation}
 with the above trial wave functions.   
\begin{enumerate}
\item[a)] Find closed form expressions for the local energy (see below) for the above 
two trial wave functions and explain shortly how these 
trial functions satisfy 
the cusp condition when $r_1\rightarrow 0$ or
$r_2\rightarrow 0$ or  $r_{12}\rightarrow 0$.
The first wave function
\[
   \Psi_T(\mathbf{r}_1,\mathbf{r}_2) = e^{-\alpha(r_1+r_2)}
\]
gives a closed-form expression
\[ 
E_{L1} = \left(\alpha-Z\right)\left(\frac{1}{r_1}+\frac{1}{r_2}\right)+\frac{1}{r_{12}}-\alpha^2,
\]
while the second wave function results in
\[ 
E_{L2} = E_{L1}+\frac{1}{2(1+\beta r_{12})^2}\left\{\frac{\alpha(r_1+r_2)}{r_{12}}(1-\frac{\mathbf{r}_1\mathbf{r}_2}{r_1r_2})-\frac{1}{2(1+\beta r_{12})^2}-\frac{2}{r_{12}}+\frac{2\beta}{1+\beta r_{12}}\right\}
\]
\item[b)] Compute    
\begin{equation}
   \langle H \rangle =
   \frac{\int d{\bf R}\Psi^{\ast}_T({\bf R})H({\bf R})\Psi_T({\bf R})}
        {\int d{\bf R}\Psi^{\ast}_T({\bf R})\Psi_T({\bf R})},
\end{equation}
for the  helium atom using a variational Monte Carlo method employing
the Metropolis algorithm to sample over different states. 
You will have to calculate
 \begin{equation}
  \langle H \rangle =\int P({\bf R})E_L({\bf R}) d{\bf R},
\end{equation}
where $E_L$ is the local energy.  Here all calculations are performed with the 
trial wave function $\psi_{T1}({\bf r_1},{\bf r_2}, {\bf r_{12}})$ only.
Study the stability of your calculation as function of the number of Monte Carlo
samples and compare these results with the exact variational result
 \begin{equation}
  \langle H \rangle =\alpha^2-2\alpha\left(Z-\frac{5}{16}\right).
\end{equation}
Your Monte Carlo moves are determined by
\begin{equation}
   {\bf R}' = {\bf R} +\delta \times r,
\end{equation}
where $r$ is a random number from the uniform distribution and $\delta$
a chosen step length.
In solving this exercise you need to devise an algorithm which finds
an optimal value of $\delta$ for each variational parameter $\alpha$,
resulting in roughly $50\%$ accepted moves.

Give a physical  interpretation of the best value of $\alpha$.
Make a plot of the variance as a function of the number of Monte Carlo
cycles.  You should parallelize your code as you did in project 3.
\item[c)] 
Use thereafter the  optimal value for  $\alpha$ as a starting point for 
computing the ground state energy of the helium atom using the trial wave functions
$\psi_{T2}({\bf r_1},{\bf r_2}, {\bf r_{12}})$. 
In this case you need to vary both $\alpha$ and $\beta$.
The strategy here is to use 
$\alpha$ from the previous exercise, [1b)] and then 
vary $\beta$ in order to find the lowest energy as function of
$\beta$. Thereafter you change 
$\alpha$ in order to see whether you find an even lower energy and so forth.

Which one of the wave functions  $\psi_{T1}({\bf r_1},{\bf r_2}, {\bf r_{12}})$ 
and  $\psi_{T2}({\bf r_1},{\bf r_2}, {\bf r_{12}})$ 
would you prefer? Give arguments for your choices.

\end{enumerate}

\subsection*{Ground state of beryllium atom}
\begin{enumerate}
\item[d)] 
The new item you need to pay attention to is the calculation of the so-called 
Slater Determinant. This is an additional complication
to your VMC calculations.
If we stick to hydrogen-like wave functions,
the trial wave function for Beryllium can be written as 
\begin{equation}
   \psi_{T}({\bf r_1},{\bf r_2}, {\bf r_3}, {\bf r_4}) = 
   Det\left(\phi_{1}({\bf r_1}),\phi_{2}({\bf r_2}),
   \phi_{3}({\bf r_3}),\phi_{4}({\bf r_4})\right)
   \prod_{i<j}^{4}\exp{\left(\frac{ar_{ij}}{(1+\beta r_{ij})}\right)}, 
\end{equation}
where $Det$ is a Slater determinant and the single-particle wave functions
are the hydrogen wave functions for the $1s$ and $2s$ orbitals. Their form
within the variational ansatz are given by
\begin{equation}
\phi_{1s}({\bf r_i}) = e^{-\alpha r_i},
\end{equation}
and 
\begin{equation}
\phi_{2s}({\bf r_i}) = \left(1-\alpha r_i/2\right)e^{-\alpha r_i/2}.
\end{equation}
You can approximate the Slater determinant for the ground state of the Beryllium atom
by writing it out as
\begin{equation}
   \psi_{T}({\bf r_1},{\bf r_2}, {\bf r_3}, {\bf r_4}) \propto 
\left(\phi_{1s}({\bf r_1})\phi_{2s}({\bf r_2})-\phi_{1s}({\bf r_2})\phi_{2s}({\bf r_1})\right)
\left(\phi_{1s}({\bf r_3})\phi_{2s}({\bf r_4})-\phi_{1s}({\bf r_4})\phi_{2s}({\bf r_3})\right).
\end{equation}
Here you can see a simple code example which implements the above expression
\begin{lstlisting}
 for (i = 0; i < number_particles; i++) {
    argument[i] = 0.0;
    r_single_particle = 0;
    for (j = 0; j < dimension; j++) {
      r_single_particle  += r[i][j]*r[i][j];
    }
    argument[i] = sqrt(r_single_particle);
  }
// Slater determinant, no factors as they vanish in Metropolis ratio
wf  = (psi1s(argument[0])*psi2s(argument[1])
       -psi1s(argument[1])*psi2s(argument[0]))*
      (psi1s(argument[2])*psi2s(argument[3])
       -psi1s(argument[3])*psi2s(argument[2]));
\end{lstlisting}
For beryllium we can easily implement the explicit evaluation of the Slater determinant.
The derivatives of the single-particle wave functions can be computed analytically and you should consider
using the closed form expression for the local energy (not mandatory, you can use numerical derivatives as well although a closed form expressions speeds up your code).

For the correlation part 
\[
\Psi_C=\prod_{i< j}g(r_{ij})= \exp{\left\{\sum_{i<j}\frac{ar_{ij}}{1+\beta r_{ij}}\right\}},
\]
we need to take into account whether electrons have equal or opposite spins since we have to obey the
electron-electron cusp condition as well.  For Beryllium you can fix electrons 1 and 2 to have spin up while
electrons 3 and 4 have spin down.
When the electrons have  equal spins 
\[
a= 1/4,
\]
while for opposite spins (as for the ground state of  helium)
\[
a= 1/2.
\] 
Set up the Hamiltonian of the system and
compute the ground state energy and compare with the 'experimental' value $E_{\mathrm{Be}}=-14.667$a.u. Run calculations with and without the correlation part (Jastrow factor). Comment your results. 
\item[e)]  This part is optional but gives you an additional 20\% on the final score!  
Compute the one-body density function for helium with an electron in the $1s$ orbital using the 
optimal (at the energy minimum) trial wave functions $\psi_{T1}$ and $\psi_{T2}$. 
The one-body density is given by 
\[
\rho({\bf r}_1)_{1s} = \int d {\bf r}_2 \left |\psi_{Ti}({\bf r}_1,{\bf r}_2)\right|^2.
\]
The variable ${\bf r}_1$ is in cartesian coordinates. Since the $1s$ state is isotropic (no angle-dependence, as applies for all $s$ waves), there is no angle dependence.  Make a plot of the one-body density as function
of $r_1$  only for the two trial wave functions.
Comment your results and discuss also the results in case your wave function is based on only the product
of two hydrogen-like wave functions. In this case, there is no repulsion between the two electrons and
we have what is called a non-interacting case. The difference between this case and the above trial wave functions tells a story about the role of the interaction between two electrons and correlations built upon that.
The non-interacting wave function is
\[
\psi_{\mathrm{non-interact}}({\bf r}_1,{\bf r}_2)= (2\left(\frac{Z}{\sqrt{4\pi}a_0}\right)^{3/2}\exp{(-Zr_1/a_0)}) (2\left(\frac{Z}{\sqrt{4\pi}a_0}\right)^{3/2}\exp{(-Zr_2/a_0)}).
\]
Here $Z=2$ and $a_0$ is the Bohr radius, which is set to one in our calculations above due to our
choice of dimensionless variables.
Plot the resulting one-body density from this wave function as well and comment your results.
\item[f)]  This part is also optional but gives you an additional 10\% on the final score!  
Compute the one-body density for the optimal wave function of the beryllium atom, with and without
the correlation part in the wave function. You will now see two peaks, one corresponding to the $1s$ part of the Slater determinant and one corresponding to the $2s$ part. 
Comment your results. Are there differences between the role 
of correlations in the helium atom and the beryllium atom?
\end{enumerate}
\section*{Literature}
\begin{enumerate}
\item S.A.~Alexander and R.L.~Coldwell,
Int.~Journal of Quantum Chemistry, {\bf 63} (1997) 1001.
\item B.~L.~Hammond, W.~A.~Lester and P.~J.~Reynolds, Monte Carlo methods
in Ab Inition Quantum Chemistry, World Scientific, Singapore, 1994, chapters
2-5 and appendix B.
\item B.H.~Bransden and C.J.~Joachain, Physics of Atoms and molecules,
Longman, 1986. Chapters 6, 7 and 9.
\item C.J.~Umrigar, K.G.~Wilson and J.W.~Wilkins, Phys.~Rev.~Lett.~{\bf 60}
(1988) 1719. 



\end{enumerate}


\end{document}








