\documentclass[10pt]{article}
\usepackage{a4wide}
\usepackage{listings}
\usepackage{graphicx}
\usepackage{times}
\usepackage{amsfonts}
\usepackage{amssymb}


\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 {\bf only your candidate number} on the first
page of the report and state clearly that this is your report for project 5 of FYS3150/FYS4150, fall 2013. 
There will be a box marked 'FYS3150/FYS4150' at the reception of the Department of Physics (room FV128). 

\section*{Project 5, Numerical integration and variational Monte Carlo methods, deadline   Monday December 9, 12pm (noon)}


The first part of this project deals with numerical integration of a six-dimensional
integral. This integral represents the expectation value of the Coulomb interaction between two electrons in a harmonic oscillator well. The integral appears in many quantum mechanical applications.
We will employ both Gauss-Legendre and Gauss-Hermite 
quadrature and Monte-Carlo integration.
Furthermore, you will need to parallelize your codes.

The second part extends these calculations to a variational Monte Carlo (VMC) calculation of the energy for two
interacting electrons in a three-dimensional quantum dot well. You can then use project 2 to benchmark your 
VMC calculations against numerically exact diagonalization results. 

The task of the first part is to integrate first in a brute force manner a six-dimensional integral which is used
to determine the ground state correlation energy between two electrons 
confined to  a three-dimensional harmonic oscillator potential.  


We assume that the wave function of each electron can be modelled like the single-particle
wave function of an electron in a three-dimensional harmonic oscillator. 
The single-particle wave function  for an electron $i$ in the 
lowest harmonic oscillator state $n_x=n_y=n_z=0$
is given in terms of a dimensionless position variable    (the wave function is not properly normalized)
\[
   {\bf r}_i =  x_i {\bf e}_x + y_i {\bf e}_y +z_i {\bf e}_z ,
\]
as
\[
   \psi_{0}({\bf r}_i)  =   \exp{-(\alpha^2 r_i^2/2)},
\]
where $\alpha=\sqrt{m\omega/\hbar}$ is a parameter related to the mass of the particle $m$ and the oscillator frequency
$\omega$ and 
\[
r_i^2 = x_i^2+y_i^2+z_i^2.
\]
We will simply fix $\alpha=1$ in the first part. We will also stay with cartesian coordinates for all exercises.
It is fully possible to switch to spherical coordinates. One would then need to replace the Hermite polynomials
with Laguerre polynomials.

The ansatz for the wave function for two electrons is then given by the product of two 
harmonic oscillator  wave functions in their lowest oscillator energy
\[
   \Psi({\bf r}_1,{\bf r}_2)  =   \exp{-(\alpha^2 (r_1^2+r_2^2))/2}.
\]
This function is not properly normalized. 
Note that it is not possible to find a closed-form  solution to Schr\"odinger's equation for 
two interacting electrons in a three-dimensional harmonic oscillator well. 

The integral we need to solve is the quantum mechanical expectation value of the correlation
energy between two electrons which repel each other via the classical Coulomb interaction, namely
\begin{equation}\label{eq:correlationenergy}
   \langle \frac{1}{|{\bf r}_1-{\bf r}_2|} \rangle =
   \int_{-\infty}^{\infty} d{\bf r}_1d{\bf r}_2  \exp{-(\alpha^2 (r_1^2+r_2^2))}\frac{1}{|{\bf r}_1-{\bf r}_2|}.
\end{equation}
Note again that our wave function is not normalized. There is a normalization factor missing, but for this project
we don't need to worry about that.

\begin{enumerate}
\item[a)] Use Gauss-Legendre quadrature and compute the integral by integrating 
for each variable $x_1$, $y_1$, $z_1$, $x_2$, $y_2$ and $z_2$ from $-\infty$ to $\infty$.
How many mesh points do you need before the results converges at the level of the third 
leading digit?  Hint:  the single-particle wave function $\exp{-(\alpha^2 r_i^2)/2}$  is more or less zero at
$r_i \approx ?$ (find the appropriate limit).  
You can therefore replace the integration limits $-\infty$ and $\infty$ with 
$-?$ and $?$, respectively.  You need to check that this approximation is satisfactory, that is, make a plot
of the function and check if the abovementioned limits are appropriate.
\item[b)]   The Legendre polynomials are defined for $x\in [-1,1]$. The previous exercise gave a very unsatisfactory ad hoc procedure. We wish to improve our results. It can therefore be useful to change to another type of polynomials
and employ the Hermite polynomials. 
The Hermite polynomials are defined for $x\in (-\infty,\infty)$ and carry a weigth function $\exp{-(ar_i^2)}$, with $a$ being a constant.
Replace the Gauss-Legendre 
approach in a) with Hermite polynomials.  The function GaussHermite.cpp can be found at the webpage of the course, see under the project link.
Do your results improve? Compare with the results from a).

{\bf Important notice for c++ programmers:}, the function which computes the Gauss-Hermite integration points and weights returns arrays which start at $1$ and end $n$ instead of the default values $0$ and $n-1$. You need to declare an array of length $n+1$. 

\item[c)] Compute the same integral but now with brute force Monte Carlo integration
and compare your results with those from the previous points. Discuss the differences.
With bruce force we mean that you should use the uniform distribution.
\item[d)] Improve your brute force Monte Carlo calculation by using importance sampling.
Hint: use a gaussian distribution. Chapter 11.5 of the lecture notes has an example for a six-dimensional 
integral with a Gaussian distribution. The function GaussianDeviate.cpp included with this project can be used.
Does the variance decrease? Does the CPU time used compared with the brute force 
Monte Carlo decrease in order to achieve the same accuracy? Comment your results
and make a list over the time each method uses. Compare the results also.
\item[e)] Finally, for the last exercise you should 
parallelize your code using MPI or openMP and run either on your laptop or the machines at the computer laboratory. Comment these results as well. In particular, we want to see whether you achieve an optimal speed-up or not.
\end{enumerate}

Our next step is to perform a VMC calculation of the ground state energy of two interacting electrons
confined to a three-dimensional harmonic oscillator well.
We label $r_1$ the position of electron 1 and similarly
$r_2$ the position of electron 2. We will set the mass of the electrons $m=1$, the oscillator frequency $\omega=1$,
and $\hbar=c=1$. We will thus use so-called atomic units.  
The contribution to the potential energy from the harmonic oscillator potential for the two electrons is
\begin{equation}
\frac{1}{2}\left(r_1^2+r_2^2\right).
\end{equation}
Adding 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{1}{2}\left(r_1^2+r_2^2\right)+\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{1}{2}\left(r_1^2+r_2^2\right)+\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. 

The basic wave functions we will employ in this part are 
\begin{equation}
   \psi_{T1}({\bf r_1},{\bf r_2}, {\bf r_{12}}) = 
   \exp{\left(-\alpha^2(r_1^2+r_2^2)/2\right)}, 
\end{equation}
and
\begin{equation}
   \psi_{T2}({\bf r_1},{\bf r_2}, {\bf r_{12}}) = 
   \exp{\left(-\alpha^2(r_1^2+r_2^2)/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[f)] 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 two electrons in a three-dimensional Harmonic oscillator well 
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 numerically exact calculations from project 2.
In order to compare the numbers you need to figure out the relation between $\alpha$ and
the frequencies you used in project 2. Another useful test is to compare the results 
without the Coulomb repulsion between two electrons. Then you have simple closed form expressions for the 
total energy. The energy is then simply the sum of two non-interacting electrons trapped in a three-dimensional
harmonic oscillator well. 

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.  For which values of $\alpha$ should you perform the variations?
Plot the variational energy as a function of $\alpha$ and localize the variational minimum.

Make a plot of the variance as a function of the number of Monte Carlo
cycles.  You should parallelize your code.
\item[g)] 
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, [1f)] 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.

\item[h)] This part is optional but gives you an additional 30\% on the final score!  
Find closed form expressions for the local energy  for the above 
two trial wave functions. Implement the closed-form expressions for the local energy and compare the speed-up in CPU time
with respect to the numerical derivation involved in the calculation of kinetic energy. How a large a speed-up do you obtain?

\end{enumerate}

Variational Monte Carlo methods will be discussed during week 47. 
\begin{thebibliography}{100} 
\bibitem{hammond} 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.
\end{thebibliography} 
\end{document}



\end{document}



