This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
In this project we will determine the ground state correlation energy between two electrons in an helium atom by calculating a six-dimensional integral that appears in many quantum mechanical applications. The methods we will use are Gauss-Legendre and Gauss-Laguerre quadrature and Monte-Carlo integration.
@@ -73,7 +73,7 @@ \section{Introduction}
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
Gaussian quadrature (hereafter GQ) is a method that solves integrals with excellent results, giving high precision for few integration points compared to simpler integration methods such as Newton-Cotes quadrature. The basic idea behind the method is to approximate the given function by a polynomial of degree $2N-1$ so that we can calculate our integral in the following way:
where $\omega$ are the weights given by the weight function $W(x)$ and $x$ are the chosen mesh points. The theory behind GQ is to obtain an arbitrary weight $\omega$, which will not be equally spaced, through the use of orthogonal polynomials, such as Legende and Laguerre polynomials.
Similarly, the Laguerre polynomials are solutions to the differential equation arising in for example the solution of the \textit{radial} Schr\"{o}dinger's equation as described above. The Laguerre polynomials are defined as
Our integral in Eq. (\ref{eq:correlationenergy}) can be rewritten with spherical coordinates to better deal with the infinite limits. Then the integral in Eq. (\ref{eq:correlationenergy}) becomes:
The radial part of the integral has limits [0,$\infty$) and we can use Laguerre polynomials with a weight function $W(x) = r^{\alpha} e^{-r}$. But first we need to do a variable change where we substitute $r_1 = u_1/2\alpha$ and $r_2 = u_2/2\alpha$, and $\mathrm{d}r_1 = \mathrm{d}u_1/2\alpha$ and $\mathrm{d}r_2 = \mathrm{d}u_2/2\alpha$. In our case $\alpha = 2$. Thus the integral becomes:
The radial part of the integral has limits [0,$\infty$) and we can use Laguerre polynomials with a weight function $W(x) = r^{\alpha} e^{-r}$. But first we need to do a variable change where we substitute $r_1 = u_1/2\alpha$ and $r_2 = u_2/2\alpha$, and $dr_1 = du_1/2\alpha$ and $dr_2 = du_2/2\alpha$. In our case $\alpha = 2$. Thus the integral becomes:
\begin{equation}
I = \frac{1}{1024} \int_0^{\pi} \sin\theta_1 \sin\theta_2 \mathrm{d}\theta_1 \mathrm{d}\theta_2
The Monte Carlo method is a statistical simulation method which can be used for systems that are described by their probability distribution functions (PDFs). The Monte Carlo method proceeds by random samplings from the PDF. The final result is taken as an average over the number of simulations, and multiplied with the Jacobi determinant of the change of variables.
\begin{equation}
-l + 2lx_i,
\end{equation}
where $l$ is the limit. The Jacobi determinant for this is $(2l)^6$, as we have 6 variables.
The Monte Carlo method is a statistical simulation method which can be used for systems that are described by their probability distribution functions (PDFs). The Monte Carlo method proceeds by random samplings from the PDF. The final result is taken as an average over the number of simulations.
We also want to run the Monte Carlo method by using importance sampling, as this should give better precision.
For this case we return to spherical coordinates,
\begin{equation}
I = \int\limits_0^{\infty}\mathrm{d}r_1\mathrm{d}r_2\int\limits_0^{\pi}\mathrm{d}\theta_1\mathrm{d}\theta_2\int\limits_0^{2\pi}\mathrm{d}\phi_1\mathrm{d}\phi_2
\end{equation}, where $r_{12}$ is given in eq.~(\ref{eq:r12}).
For $r_i$ we want to use exponentially distributed random numbers.
%In the regular Monte Carlo method we use uniformly distributed random numbers, but when considering importance sampling we will instead use exponentially distributed random numbers.
%We will also switch to spherical coordinates.
This is done by changing variables as follows
\begin{align*}
r_i &\rightarrow -\frac{1}{4}\ln(1-x_i) \\
\theta_i &\rightarrow\pi x_i \\
\phi_i &\rightarrow 2\pi x_i,
\end{align*}
where $x_i$ is a random number between 0 and 1.
The Jocobi determinant is in this case given by $(1/4)^2\times\pi^2\times(2\pi)^2 = \pi^4/4$.
Lastly, we calculate the variance, which is given by
We also want to run the Monte Carlo method by using importance sampling. In the regular Monte Carlo method we use uniformly distributed random numbers, but when considering importance sampling we will instead use exponentially distributed random numbers. We will also switch to spherical coordinates.