This website contains the QMC tutorial of the 2023 LTTC winter school Tutorials in Theoretical Chemistry.
We propose different exercises to understand quantum Monte Carlo (QMC) methods. In the first section, we start with the computation of the energy of a hydrogen atom using numerical integration. The goal of this section is to familiarize yourself with the concept of local energy. Then, we introduce the variational Monte Carlo (VMC) method which computes a statistical estimate of the expectation value of the energy associated with a given wave function, and apply this approach to the hydrogen atom. Finally, we present the diffusion Monte Carlo (DMC) method which we use here to estimate the exact energy of the hydrogen atom and of the H_2 molecule, starting from an approximate wave function.
Code examples will be given in Python3, C and Fortran. You can use whatever language you prefer to write the programs.
We consider the stationary solution of the Schrödinger equation, so
the wave functions considered here are real: for an
All the quantities are expressed in atomic units (energies, coordinates, etc).
For a given system with Hamiltonian
where
The electronic energy of a system,
\begin{eqnarray*}
E & = & \frac{\langle Ψ| \hat{H} | Ψ\rangle}{\langle Ψ |Ψ \rangle}
= \frac{∫ Ψ(\mathbf{r})\, \hat{H} Ψ(\mathbf{r})\, d\mathbf{r}}{∫ |Ψ(\mathbf{r}) |^2 d\mathbf{r}}
& = & \frac{∫ |Ψ(\mathbf{r})|^2\, \frac{\hat{H} Ψ(\mathbf{r})}{Ψ(\mathbf{r})}\,d\mathbf{r}}{∫ |Ψ(\mathbf{r}) |^2 d\mathbf{r}}
= \frac{∫ |Ψ(\mathbf{r})|^2\, E_L(\mathbf{r})\,d\mathbf{r}}{∫ |Ψ(\mathbf{r}) |^2 d\mathbf{r}}
\end{eqnarray*}
For few dimensions, one can easily compute
To this aim, recall that the probabilistic expected value of an
arbitrary function
where the probability density function
Similarly, we can view the energy of a system,
where the probability density is given by the square of the wave function:
If we can sample $N\rm MC$ configurations
Here is an animation of the local energy in the N2 molecule along a QMC trajectory:
In this section, we consider the hydrogen atom with the following wave function:
We will first verify that, for a particular value of
To do that, we will compute the local energy and check whether it is constant.
You will now program all quantities needed to compute the local energy of the Hydrogen atom for the given wave function.
Write all the functions of this section in a single file :
hydrogen.py if you use Python, hydrogen.F90 if you use
Fortran, or hydrogen.c and hydrogen.h if you use C.
Test-driven development (TDD) is a software development process in which automated tests are written before the code they are testing. The tests are then run to ensure that the code is working as expected, and the developer continues to write code until all tests pass. The benefits of TDD include:
- Early detection of bugs: By writing tests before the code, developers can identify and fix bugs early on in the development process, before the code becomes too complex to easily debug.
- Improved code quality: Writing tests forces developers to think about the behavior and functionality of the code they are writing, leading to more robust and maintainable code.
- Better documentation: Tests act as a form of documentation for the code, making it clear what the code is supposed to do and how it should behave.
#!/usr/bin/env python3
import numpy as np
def potential(r):
# TODO
<<test_p>>If the extension of your file is .F90 and not .f90, the C
preprocessor will be called before compiling. This enables the
possibility to have conditional compilation with #ifdef
statements, activated with the -D compiler option.
To compile your code and activate the test, use:
gfortran -DTEST_H hydrogen.F90 -o test_hydrogendouble precision function potential(r)
implicit none
double precision, intent(in) :: r(3)
! TODO
end function potential
<<test_f>>Compile your code with the -DTEST_H option. It will activate
the creation of the main function that will test your functions.
Don’t forget to use -lm to link with the math library.
gcc -DTEST_H hydrogen.c -lm -o test_hydrogen#include <stdio.h> // printf
#include <math.h> // sqrt
#include <stdlib.h> // exit
#include <assert.h> // assert
double potential(double r[3]) {
// TODO
}
<<test_c>>#!/usr/bin/env python3
import numpy as np
def potential(r):
distance = np.sqrt(np.dot(r,r))
if distance == 0:
print("potential at r=0 diverges")
return -float("inf")
return -1. / distance
<<test_p>>double precision function potential(r)
implicit none
double precision, intent(in) :: r(3)
double precision :: distance
distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
if (distance > 0.d0) then
potential = -1.d0 / distance
else
print *, 'Warning: potential at r=0.d0 diverges'
potential = -huge(1.d0)
end if
end function potential
<<test_f>>#include <stdio.h> // printf
#include <math.h> // sqrt
#include <stdlib.h> // exit
#include <assert.h> // assert
double potential(double r[3]) {
double distance;
distance = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
if (distance > 0) {
return -1.0 / distance;
} else {
printf("Warning: potential at r=0 diverges\n");
return -HUGE_VAL;
}
}
<<test_c>>def psi(a, r):
# TODOdouble precision function psi(a, r)
implicit none
double precision, intent(in) :: a, r(3)
! TODO
end function psidouble psi(double a, double r[3]) {
// TODO
}def psi(a, r):
return np.exp(-a*np.sqrt(np.dot(r,r)))double precision function psi(a, r)
implicit none
double precision, intent(in) :: a, r(3)
psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
end function psidouble psi(double a, double r[3]) {
double distance;
distance = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
return exp(-a * distance);
}The local kinetic energy is defined as
We differentiate
\[ Ψ(\mathbf{r}) = exp(-a\,|\mathbf{r}|) \] \[\frac{∂ Ψ}{∂ x} = \frac{∂ Ψ}{∂ |\mathbf{r}|} \frac{∂ |\mathbf{r}|}{∂ x} = - \frac{a\,x}{|\mathbf{r}|} Ψ(\mathbf{r}) \]
and we differentiate a second time:
The Laplacian operator
Therefore, the local kinetic energy is $$ T_L (\mathbf{r}) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) $$
def kinetic(a,r):
# TODOdouble precision function kinetic(a,r)
implicit none
double precision, intent(in) :: a, r(3)
! TODO
end function kineticdouble kinetic(double a, double r[3]) {
//TODO
}def kinetic(a,r):
distance = np.sqrt(np.dot(r,r))
if distance > 0:
dinv = 1./distance
else:
print ('Warning: kinetic energy diverges at r=0')
dinv = float("inf")
return a * (dinv - 0.5 * a)double precision function kinetic(a,r)
implicit none
double precision, intent(in) :: a, r(3)
double precision :: distance
distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
if (distance > 0.d0) then
kinetic = a * (1.d0 / distance - 0.5d0 * a)
else
print *, 'Warning: kinetic energy diverges at r=0'
kinetic = a * (huge(1.d0) - 0.5d0 * a)
end if
end function kineticdouble kinetic(double a, double r[3]) {
double distance;
distance = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
if (distance > 0) {
return a * (1.0 / distance - 0.5 * a);
} else {
printf("Warning: kinetic energy diverges at r=0\n");
return a * (HUGE_VAL - 0.5 * a);
}
}def e_loc(a,r):
#TODOdouble precision function e_loc(a,r)
implicit none
double precision, intent(in) :: a, r(3)
double precision, external :: kinetic
double precision, external :: potential
! TODO
end function e_locdouble e_loc(double a, double r[3]) {
// TODO
}def e_loc(a,r):
return kinetic(a,r) + potential(r)double precision function e_loc(a,r)
implicit none
double precision, intent(in) :: a, r(3)
double precision, external :: kinetic
double precision, external :: potential
e_loc = kinetic(a,r) + potential(r)
end function e_locdouble e_loc(double a, double r[3]) {
return kinetic(a, r) + potential(r);
}\begin{eqnarray*}
E &=& \frac{\hat{H} Ψ}{Ψ} = - \frac{1}{2} \frac{Δ Ψ}{Ψ} -
\frac{1}{|\mathbf{r}|}
&=& -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) -
\frac{1}{|\mathbf{r}|} \
&=&
-\frac{1}{2} a^2 + \frac{a-1}{\mathbf{|r|}}
\end{eqnarray*}
The program you will write in this section will be written in
another file (plot_hydrogen.py, plot_hydrogen.F90 or plot_hydrogen.c for
example).
It will use the functions previously defined.
If you use C, don’t forget to write the header file corresponding
to the functions defined in the previous section.
In Python, you should put at the beginning of the file
#!/usr/bin/env python3
from hydrogen import e_locto be able to use the e_loc function of the hydrogen.py file.
In Fortran, you will need to compile all the source files together:
gfortran hydrogen.F90 plot_hydrogen.F90 -o plot_hydrogenSimilarly, in C
gcc hydrogen.c plot_hydrogen.c -lm -o plot_hydrogen#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
from hydrogen import e_loc
x=np.linspace(-5,5)
plt.figure(figsize=(10,5))
# TODO
plt.tight_layout()
plt.legend()
plt.savefig("plot_py.png")program plot
implicit none
double precision, external :: e_loc
double precision :: x(50), dx
integer :: i, j
dx = 10.d0/(size(x)-1)
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
! TODO
end program plotTo compile and run:
gfortran hydrogen.F90 plot_hydrogen.F90 -o plot_hydrogen
./plot_hydrogen > data#include <stdio.h>
#include <math.h>
#include "hydrogen.h"
#define NPOINTS 50
#define NEXPO 6
int main() {
double x[NPOINTS], energy, dx, r[3];
double a[NEXPO] = { 0.1, 0.2, 0.5, 1.0, 1.5, 2.0 };
int i, j;
dx = 10.0/(NPOINTS-1);
for (i = 0; i < NPOINTS; i++) {
x[i] = -5.0 + i*dx;
}
// TODO
return 0;
}To compile and run:
gcc hydrogen.c plot_hydrogen.c -lm -o plot_hydrogen
./plot_hydrogen > dataPlotting for Fortran of C
Plot the data using Gnuplot:
set grid
set xrange [-5:5]
set yrange [-2:1]
plot './data' index 0 using 1:2 with lines title 'a=0.1', \
'./data' index 1 using 1:2 with lines title 'a=0.2', \
'./data' index 2 using 1:2 with lines title 'a=0.5', \
'./data' index 3 using 1:2 with lines title 'a=1.0', \
'./data' index 4 using 1:2 with lines title 'a=1.5', \
'./data' index 5 using 1:2 with lines title 'a=2.0'#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
from hydrogen import e_loc
x=np.linspace(-5,5)
plt.figure(figsize=(10,5))
for a in [0.1, 0.2, 0.5, 1., 1.5, 2.]:
y=np.array([ e_loc(a, np.array([t,0.,0.]) ) for t in x])
plt.plot(x,y,label=f"a={a}")
plt.tight_layout()
plt.legend()
plt.savefig("plot_py.png")
program plot
implicit none
double precision, external :: e_loc
double precision :: x(50), energy, dx, r(3), a(6)
integer :: i, j
a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
dx = 10.d0/(size(x)-1)
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
r(:) = 0.d0
do j=1,size(a)
print *, '# a=', a(j)
do i=1,size(x)
r(1) = x(i)
energy = e_loc( a(j), r )
print *, x(i), energy
end do
print *, ''
print *, ''
end do
end program plotdouble potential (double r[3]);
double psi (double a, double r[3]);
double kinetic (double a, double r[3]);
double e_loc (double a, double r[3]);#include <stdio.h>
#include <math.h>
#include "hydrogen.h"
#define NPOINTS 50
#define NEXPO 6
int main() {
double x[NPOINTS], energy, dx, r[3];
double a[NEXPO] = { 0.1, 0.2, 0.5, 1.0, 1.5, 2.0 };
int i, j;
dx = 10.0/(NPOINTS-1);
for (i = 0; i < NPOINTS; i++) {
x[i] = -5.0 + i*dx;
}
for (i = 0; i < 3; i++) {
r[i] = 0.0;
}
for (j = 0; j < NEXPO; j++) {
printf("# a=%f\n", a[j]);
for (i = 0; i < NPOINTS; i++) {
r[0] = x[i];
energy = e_loc(a[j], r);
printf("%f %f\n", x[i], energy);
}
printf("\n\n");
}
return 0;
}
If the space is discretized in small volume elements
#!/usr/bin/env python3
import numpy as np
from hydrogen import e_loc, psi
interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3
r = np.array([0.,0.,0.])
for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
# TODO
print(f"a = {a} \t E = {E}") program energy_hydrogen
implicit none
double precision, external :: e_loc, psi
double precision :: x(50), w, delta, energy, dx, r(3), a(6), normalization
integer :: i, k, l, j
a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
dx = 10.d0/(size(x)-1)
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
do j=1,size(a)
! TODO
print *, 'a = ', a(j), ' E = ', energy
end do
end program energy_hydrogenTo compile the Fortran code and run it:
gfortran hydrogen.F90 energy_hydrogen.F90 -o energy_hydrogen
./energy_hydrogen #include <stdio.h>
#include <math.h>
#include "hydrogen.h"
#define NPOINTS 50
#define NEXPO 6
int main() {
double x[NPOINTS], energy, dx, r[3], delta, normalization, w;
double a[NEXPO] = { 0.1, 0.2, 0.5, 1.0, 1.5, 2.0 };
dx = 10.0/(NPOINTS-1);
for (int i = 0; i < NPOINTS; i++) {
x[i] = -5.0 + i*dx;
}
for (int j = 0; j < NEXPO; j++) {
// TODO
printf("a = %f E = %f\n", a[j], energy);
}
}To compile the C code and run it:
gcc hydrogen.c energy_hydrogen.c -lm -o energy_hydrogen
./energy_hydrogen Hints if you are stuck
The program starts by defining some variables and arrays, including
an array a that contains 6 different values of the parameter a
which will be used in the e_loc and psi functions to calculate
the local energy and wave function respectively.
The program then calculates the value of dx, which is the step size in
x that contains 50 equally spaced points
between -5 and 5. The program sets all elements of the r array to 0,
and then enters a nested loop structure. The outer loop iterates over
the values of a in the a array, and the next three loops iterate
over the values of x in the x array for the three dimensions. For
each value of a and x, the program sets the first element of the
r array to the current value of x, calls the psi function to
calculate the wave function, calls the e_loc function to calculate
the local energy, and then accumulates the energy and the
normalization factor.
At the end of the outer loop, the program calculates the final energy
by dividing the accumulated energy by the accumulated normalization
factor, and prints the value of a and the corresponding energy.
#!/usr/bin/env python3
import numpy as np
from hydrogen import e_loc, psi
interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3
r = np.array([0.,0.,0.])
for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
E = 0.
normalization = 0.
for x in interval:
r[0] = x
for y in interval:
r[1] = y
for z in interval:
r[2] = z
w = psi(a,r)
w = w * w * delta
E += w * e_loc(a,r)
normalization += w
E = E / normalization
print(f"a = {a} \t E = {E}") a = 0.1 E = -0.24518438948809218 a = 0.2 E = -0.26966057967803525 a = 0.5 E = -0.3856357612517407 a = 0.9 E = -0.49435709786716214 a = 1.0 E = -0.5 a = 1.5 E = -0.39242967082602226 a = 2.0 E = -0.08086980667844901
program energy_hydrogen
implicit none
double precision, external :: e_loc, psi
double precision :: x(50), w, delta, energy, dx, r(3), a(6), normalization
integer :: i, k, l, j
a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
dx = 10.d0/(size(x)-1)
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
delta = dx**3
r(:) = 0.d0
do j=1,size(a)
energy = 0.d0
normalization = 0.d0
do i=1,size(x)
r(1) = x(i)
do k=1,size(x)
r(2) = x(k)
do l=1,size(x)
r(3) = x(l)
w = psi(a(j),r)
w = w * w * delta
energy = energy + w * e_loc(a(j), r)
normalization = normalization + w
end do
end do
end do
energy = energy / normalization
print *, 'a = ', a(j), ' E = ', energy
end do
end program energy_hydrogena = 0.10000000000000001 E = -0.24518438948809140 a = 0.20000000000000001 E = -0.26966057967803236 a = 0.50000000000000000 E = -0.38563576125173815 a = 1.0000000000000000 E = -0.50000000000000000 a = 1.5000000000000000 E = -0.39242967082602065 a = 2.0000000000000000 E = -8.0869806678448772E-002
#include <stdio.h>
#include <math.h>
#include "hydrogen.h"
#define NPOINTS 50
#define NEXPO 6
int main() {
double x[NPOINTS], energy, dx, r[3], delta, normalization, w;
double a[NEXPO] = { 0.1, 0.2, 0.5, 1.0, 1.5, 2.0 };
dx = 10.0/(NPOINTS-1);
for (int i = 0; i < NPOINTS; i++) {
x[i] = -5.0 + i*dx;
}
delta = dx*dx*dx;
for (int i = 0; i < 3; i++) {
r[i] = 0.0;
}
for (int j = 0; j < NEXPO; j++) {
energy = 0.0;
normalization = 0.0;
for (int i = 0; i < NPOINTS; i++) {
r[0] = x[i];
for (int k = 0; k < NPOINTS; k++) {
r[1] = x[k];
for (int l = 0; l < NPOINTS; l++) {
r[2] = x[l];
w = psi(a[j], r);
w = w*w*delta;
energy += w*e_loc(a[j], r);
normalization += w;
}
}
}
energy = energy/normalization;
printf("a = %f E = %f\n", a[j], energy);
}
}a = 0.100000 E = -0.245184 a = 0.200000 E = -0.269661 a = 0.500000 E = -0.385636 a = 1.000000 E = -0.500000 a = 1.500000 E = -0.392430 a = 2.000000 E = -0.080870
The variance of the local energy is a functional of
$$ σ^2(E_L) = \frac{∫ |Ψ(\mathbf{r})|^2\, \left[ E_L(\mathbf{r}) - E \right]^2 \, d\mathbf{r}}{∫ |Ψ(\mathbf{r}) |^2 d\mathbf{r}} $$ which can be simplified as
If the local energy is constant (i.e.
\begin{eqnarray*}
\langle (E - \bar{E})^2 \rangle & = &
\langle E^2 - 2 E \bar{E} + \bar{E}^2 \rangle
&=& \langle E^2 \rangle - 2 \langle E \bar{E} \rangle + \langle \bar{E}^2 \rangle \
&=& \langle E^2 \rangle - 2 \langle E \rangle \bar{E} + \bar{E}^2 \
&=& \langle E^2 \rangle - 2 \bar{E}^2 + \bar{E}^2 \
&=& \langle E^2 \rangle - \bar{E}^2 \
&=& \langle E^2 \rangle - \langle E \rangle^2 \
\end{eqnarray*}
#!/usr/bin/env python3
import numpy as np from hydrogen import e_loc, psi
interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3
r = np.array([0.,0.,0.])
for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
# TODO
print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")program variance_hydrogen
implicit none
double precision :: x(50), w, delta, energy, energy2
double precision :: dx, r(3), a(6), normalization, e_tmp, s2
integer :: i, k, l, j
double precision, external :: e_loc, psi
a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
dx = 10.d0/(size(x)-1)
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
do j=1,size(a)
! TODO
print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2
end do
end program variance_hydrogenTo compile and run:
gfortran hydrogen.F90 variance_hydrogen.F90 -o variance_hydrogen
./variance_hydrogen#include <stdio.h>
#include <math.h>
#include "hydrogen.h"
#define NPOINTS 50
#define NEXPO 6
int main() {
double x[NPOINTS], energy, dx, r[3], delta, normalization, w;
double a[NEXPO] = { 0.1, 0.2, 0.5, 1.0, 1.5, 2.0 };
double energy2, e_tmp, s2;
dx = 10.0/(NPOINTS-1);
for (int i = 0; i < NPOINTS; i++) {
x[i] = -5.0 + i*dx;
}
for (int j = 0; j < NEXPO; j++) {
// TODO
printf("a = %f E = %f s2 = %f\n", a[j], energy, s2);
}
}To compile and run:
gcc hydrogen.c variance_hydrogen.c -lm -o variance_hydrogen
./variance_hydrogen#!/usr/bin/env python3
import numpy as np
from hydrogen import e_loc, psi
interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3
r = np.array([0.,0.,0.])
for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
E = 0.
E2 = 0.
normalization = 0.
for x in interval:
r[0] = x
for y in interval:
r[1] = y
for z in interval:
r[2] = z
w = psi(a,r)
w = w * w * delta
e_tmp = e_loc(a,r)
E += w * e_tmp
E2 += w * e_tmp * e_tmp
normalization += w
E = E / normalization
E2 = E2 / normalization
s2 = E2 - E**2
print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")a = 0.1 E = -0.24518439 \sigma^2 = 0.02696522 a = 0.2 E = -0.26966058 \sigma^2 = 0.03719707 a = 0.5 E = -0.38563576 \sigma^2 = 0.05318597 a = 0.9 E = -0.49435710 \sigma^2 = 0.00577812 a = 1.0 E = -0.50000000 \sigma^2 = 0.00000000 a = 1.5 E = -0.39242967 \sigma^2 = 0.31449671 a = 2.0 E = -0.08086981 \sigma^2 = 1.80688143
program variance_hydrogen
implicit none
double precision :: x(50), w, delta, energy, energy2
double precision :: dx, r(3), a(6), normalization, e_tmp, s2
integer :: i, k, l, j
double precision, external :: e_loc, psi
a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
dx = 10.d0/(size(x)-1)
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
delta = dx**3
r(:) = 0.d0
do j=1,size(a)
energy = 0.d0
energy2 = 0.d0
normalization = 0.d0
do i=1,size(x)
r(1) = x(i)
do k=1,size(x)
r(2) = x(k)
do l=1,size(x)
r(3) = x(l)
w = psi(a(j),r)
w = w * w * delta
e_tmp = e_loc(a(j), r)
energy = energy + w * e_tmp
energy2 = energy2 + w * e_tmp * e_tmp
normalization = normalization + w
end do
end do
end do
energy = energy / normalization
energy2 = energy2 / normalization
s2 = energy2 - energy*energy
print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2
end do
end program variance_hydrogena = 0.10000000000000001 E = -0.24518438948809140 s2 = 2.6965218719722767E-002 a = 0.20000000000000001 E = -0.26966057967803236 s2 = 3.7197072370201284E-002 a = 0.50000000000000000 E = -0.38563576125173815 s2 = 5.3185967578480653E-002 a = 1.0000000000000000 E = -0.50000000000000000 s2 = 0.0000000000000000 a = 1.5000000000000000 E = -0.39242967082602065 s2 = 0.31449670909172917 a = 2.0000000000000000 E = -8.0869806678448772E-002 s2 = 1.8068814270846534
#include <stdio.h>
#include <math.h>
#include "hydrogen.h"
#define NPOINTS 50
#define NEXPO 6
int main() {
double x[NPOINTS], energy, dx, r[3], delta, normalization, w;
double a[NEXPO] = { 0.1, 0.2, 0.5, 1.0, 1.5, 2.0 };
double energy2, e_tmp, s2;
dx = 10.0/(NPOINTS-1);
for (int i = 0; i < NPOINTS; i++) {
x[i] = -5.0 + i*dx;
}
delta = dx*dx*dx;
for (int i = 0; i < 3; i++) {
r[i] = 0.0;
}
for (int j = 0; j < NEXPO; j++) {
energy = 0.0;
energy2 = 0.0;
normalization = 0.0;
for (int i = 0; i < NPOINTS; i++) {
r[0] = x[i];
for (int k = 0; k < NPOINTS; k++) {
r[1] = x[k];
for (int l = 0; l < NPOINTS; l++) {
r[2] = x[l];
w = psi(a[j], r);
w = w*w*delta;
e_tmp = e_loc(a[j], r);
energy += w * e_tmp;
energy2 += w * e_tmp * e_tmp;
normalization += w;
}
}
}
energy = energy/normalization;
energy2 = energy2/normalization;
s2 = energy2 - energy*energy;
printf("a = %f E = %f s2 = %f\n", a[j], energy, s2);
}
}gcc hydrogen.c variance_hydrogen.c -lm -o variance_hydrogen
./variance_hydrogena = 0.100000 E = -0.245184 s2 = 0.026965 a = 0.200000 E = -0.269661 s2 = 0.037197 a = 0.500000 E = -0.385636 s2 = 0.053186 a = 1.000000 E = -0.500000 s2 = 0.000000 a = 1.500000 E = -0.392430 s2 = 0.314497 a = 2.000000 E = -0.080870 s2 = 1.806881
Numerical integration with deterministic methods is very efficient in low dimensions. When the number of dimensions becomes large, instead of computing the average energy as a numerical integration on a grid, it is usually more efficient to use Monte Carlo sampling.
Moreover, Monte Carlo sampling will allow us to remove the bias due to the discretization of space, and compute a statistical confidence interval.
To compute the statistical error, you need to perform
The estimate of the energy is
The variance of the average energies can be computed as
And the confidence interval is given by
#!/usr/bin/env python3
from math import sqrt
def ave_error(arr):
#TODO
return (average, error)subroutine ave_error(x,n,ave,err)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: x(n)
double precision, intent(out) :: ave, err
! TODO
end subroutine ave_error#include <stdio.h>
#include <math.h>
#include <stddef.h> // for size_t
void ave_error(double* x, size_t n, double *ave, double *err) {
// TODO
}Hints if you are stuck
Write a subroutine called ave_error that calculates the average
and error of a given array of real numbers. The subroutine takes in
three arguments: an array x of real numbers, an integer n
representing the size of the array, and two output arguments ave
and err representing the average and error of the array,
respectively.
The subroutine starts by checking if the input integer n is less
than 1. If it is, the subroutine stops and prints an error message.
If n is equal to 1, the subroutine sets the average to the first
element of the array and the error to zero. If n is greater than
1, the subroutine calculates the average of the array by dividing
the sum of the elements by the number of elements in the
array. Then it calculates the variance of the array by taking the
sum of the square of the difference between each element and the
average and dividing by n-1. Finally, it calculates the error by
taking the square root of the variance divided by n.
#!/usr/bin/env python3
from math import sqrt
def ave_error(arr):
M = len(arr)
assert(M>0)
if M == 1:
average = arr[0]
error = 0.
else:
average = sum(arr)/M
variance = 1./(M-1) * sum( [ (x - average)**2 for x in arr ] )
error = sqrt(variance/M)
return (average, error)subroutine ave_error(x,n,ave,err)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: x(n)
double precision, intent(out) :: ave, err
double precision :: variance
if (n < 1) then
stop 'n<1 in ave_error'
else if (n == 1) then
ave = x(1)
err = 0.d0
else
ave = sum(x(:)) / dble(n)
variance = sum((x(:) - ave)**2) / dble(n-1)
err = dsqrt(variance/dble(n))
endif
end subroutine ave_error#include <stddef.h> // for size_t
void ave_error(double* x, size_t n, double *ave, double *err);#include <stdio.h>
#include <math.h>
#include <stddef.h> // for size_t
void ave_error(double* x, size_t n, double *ave, double *err) {
double variance;
if (n < 1) {
printf("n<1 in ave_error\n");
return;
} else if (n == 1) {
*ave = x[0];
*err = 0.0;
} else {
double sum = 0.0;
for (int i = 0; i < n; i++) {
sum += x[i];
}
*ave = sum / (double)n;
variance = 0.0;
for (int i = 0; i < n; i++) {
double x2 = x[i] - *ave;
variance += x2*x2;
}
variance = variance / (double)(n - 1);
*err = sqrt(variance / (double)n);
}
}We will now perform our first Monte Carlo calculation to compute the energy of the hydrogen atom.
Consider again the expression of the energy
\begin{eqnarray*} E & = & \frac{∫ E_L(\mathbf{r})|Ψ(\mathbf{r})|^2\,d\mathbf{r}}{∫ |Ψ(\mathbf{r}) |^2 d\mathbf{r}}\,. \end{eqnarray*}
Clearly, the square of the wave function is a good choice of probability density to sample but we will start with something simpler and rewrite the energy as
\begin{eqnarray*} E & = & \frac{∫ E_L(\mathbf{r})\frac{|Ψ(\mathbf{r})|^2}{P(\mathbf{r})}P(\mathbf{r})\, \,d\mathbf{r}}{∫ \frac{|Ψ(\mathbf{r})|^2 }{P(\mathbf{r})}P(\mathbf{r})d\mathbf{r}}\,. \end{eqnarray*}
Here, we will sample a uniform probability
and zero outside the cube.
One Monte Carlo run will consist of $N\rm MC$ Monte Carlo iterations. At every Monte Carlo iteration:
- Draw a random point
$\mathbf{r}_i$ in the box$(-5,-5,-5) ≤ (x,y,z) ≤ (5,5,5)$ - Compute
$|Ψ(\mathbf{r}_i)|^2$ and accumulate the result in a variablenormalization - Compute
$|Ψ(\mathbf{r}_i)|^2 × E_L(\mathbf{r}_i)$ , and accumulate the result in a variableenergy
Once all the iterations have been computed, the run returns the average energy
To compute the statistical error, perform
#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a, nmax):
# TODO
a = 1.2
nmax = 100000
#TODO
print(f"E = {E} +/- {deltaE}")subroutine uniform_montecarlo(a,nmax,energy)
implicit none
double precision, intent(in) :: a
integer , intent(in) :: nmax
double precision, intent(out) :: energy
integer :: istep
double precision :: normalization, r(3), w
double precision, external :: e_loc, psi
! TODO
end subroutine uniform_montecarlo
program qmc
implicit none
double precision, parameter :: a = 1.2d0
integer , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns)
double precision :: ave, err
!TODO
print *, 'E = ', ave, '+/-', err
end program qmcgfortran hydrogen.F90 qmc_stats.F90 qmc_uniform.F90 -o qmc_uniform
./qmc_uniform#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <stddef.h> // for size_t
#include <time.h>
#include "hydrogen.h"
#include "qmc_stats.h" // for ave_error
void uniform_montecarlo(double a, size_t nmax, double *energy) {
// TODO
}
int main(void) {
#define a 1.2
#define nmax 100000
#define nruns 30
srand48(time(NULL));
// TODO
printf("E = %f +/- %f\n", ave, err);
return 0;
}Hints if you are stuck
Write first a subroutine called uniform_montecarlo that
calculates the energy of the Hydrogen atom using the Monte Carlo
method with a uniform distribution. The subroutine takes in three
arguments: a real number a, an integer nmax representing the
number of Monte Carlo steps, and an output argument energy
representing the calculated energy.
The subroutine starts by initializing the energy and normalization
factor to 0 and defines some variables such as istep, normalization,
r and w. The subroutine also makes use of two external
functions: e_loc and psi which were defined in previous
examples.
The subroutine then enters a loop that iterates for nmax times. On
each iteration, the subroutine generates three random numbers
between 0 and 1, and then uses these random numbers to calculate a
random point in 3D space between -5 and 5. The subroutine then
calls the psi function to calculate the wave function at that
point and the e_loc function to calculate the local energy at
that point. The subroutine then accumulates the energy and
normalization factor using the generated point and the results of
the psi and e_loc functions.
At the end of the loop, the subroutine calculates the final energy by dividing the accumulated energy by the accumulated normalization factor.
Then, write a Fortran program called qmc that uses the
uniform_montecarlo subroutine to estimate the energy of the
Hydrogen atom using the Monte Carlo method. The program starts by
defining some parameters: a, nmax, and nruns.
The program then defines a variable irun which is used as a counter
in a loop, an array X of length nruns to store the energies
calculated by the uniform_montecarlo subroutine, and variables ave
and err to store the average and error of the energies,
respectively.
The program then enters a loop that iterates for nruns times. On
each iteration, the program calls the uniform_montecarlo
subroutine to calculate the energy of the Hydrogen atom and stores
the result in the X array.
After the loop, the program calls the ave_error subroutine to
calculate the average and error of the energies stored in the X
array and assigns the results to ave and err variables
respectively.
Finally, the program prints the average and error of the energies.
#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a, nmax):
energy = 0.
normalization = 0.
for istep in range(nmax):
r = np.random.uniform(-5., 5., (3))
f = psi(a,r)
w = f*f
energy += w * e_loc(a,r)
normalization += w
return energy / normalization
a = 1.2
nmax = 100000
X = [MonteCarlo(a,nmax) for i in range(30)]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")E = -0.4793311279357434 +/- 0.002563797463053474
subroutine uniform_montecarlo(a,nmax,energy)
implicit none
double precision, intent(in) :: a
integer*8 , intent(in) :: nmax
double precision, intent(out) :: energy
integer*8 :: istep
double precision :: normalization, r(3), w, f
double precision, external :: e_loc, psi
energy = 0.d0
normalization = 0.d0
do istep = 1,nmax
call random_number(r)
r(:) = -5.d0 + 10.d0*r(:)
f = psi(a,r)
w = f*f
energy = energy + w * e_loc(a,r)
normalization = normalization + w
end do
energy = energy / normalization
end subroutine uniform_montecarlo
program qmc
implicit none
double precision, parameter :: a = 1.2d0
integer*8 , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns)
double precision :: ave, err
do irun=1,nruns
call uniform_montecarlo(a, nmax, X(irun))
enddo
call ave_error(X, nruns, ave, err)
print *, 'E = ', ave, '+/-', err
end program qmcE = -0.47958062416368030 +/- 3.1460315707685389E-003
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <time.h>
#include <stddef.h> // for size_t
#include "hydrogen.h"
#include "qmc_stats.h" // for ave_error
void uniform_montecarlo(double a, size_t nmax, double *energy) {
size_t istep;
double normalization, r[3], w, f;
*energy = 0.0;
normalization = 0.0;
for (istep = 0; istep < nmax; istep++) {
for (int i = 0; i < 3; i++) {
r[i] = drand48();
}
r[0] = -5.0 + 10.0 * r[0];
r[1] = -5.0 + 10.0 * r[1];
r[2] = -5.0 + 10.0 * r[2];
f = psi(a, r);
w = f*f;
*energy += w * e_loc(a, r);
normalization += w;
}
*energy = *energy / normalization;
}
int main(void) {
#define a 1.2
#define nmax 100000
#define nruns 30
double X[nruns];
double ave, err;
srand48(time(NULL));
for (size_t irun = 0; irun < nruns; irun++) {
uniform_montecarlo(a, nmax, &X[irun]);
}
ave_error(X, nruns, &ave, &err);
printf("E = %f +/- %f\n", ave, err);
return 0;
}E = -0.479507 +/- 0.001972
We will now use the square of the wave function to sample random points distributed with the probability density \[ P(\mathbf{r}) = \frac{|Ψ(\mathbf{r})|^2}{∫ |Ψ(\mathbf{r})|^2 d\mathbf{r}}\,. \]
The expression of the average energy is now simplified as the average of the local energies, since the weights are taken care of by the sampling:
To sample a chosen probability density, an efficient method is the
Metropolis-Hastings sampling algorithm. Starting from a random
initial position
according to the following algorithm.
At every step, we propose a new move according to a transition probability $T(\mathbf{r}n→\mathbf{r}n+1)$ of our choice.
For simplicity, we will move the electron in a 3-dimensional box of side
where
After having moved the electron, we add the
accept/reject step that guarantees that the distribution of the
which, for our choice of transition probability, becomes
Also note that we do not need to compute the norm of the wave function!
The algorithm is summarized as follows:
- Evaluate the local energy at
$\mathbf{r}_n$ and add it to an accumulator - Compute a new position
$\mathbf{r’} = \mathbf{r}_n + δ L\, \mathbf{u}$ - Evaluate
$Ψ(\mathbf{r}’)$ at the new position - Compute the ratio $A = \frac{\left|Ψ(\mathbf{r’})\right|^2}{\left|Ψ(\mathbf{r}n)\right|^2}$
- Draw a uniform random number
$v ∈ [0,1]$ - if
$v ≤ A$ , accept the move : set $\mathbf{r}n+1 = \mathbf{r’}$ - else, reject the move : set $\mathbf{r}n+1 = \mathbf{r}_n$
If the box is infinitely small, the ratio will be very close to one and all the steps will be accepted. However, the moves will be very correlated and you will explore the configurational space very slowly.
On the other hand, if you propose too large moves, the number of accepted steps will decrease because the ratios might become small. If the number of accepted steps is close to zero, then the space is not well sampled either.
The size of the move should be adjusted so that it is as large as possible, keeping the number of accepted steps not too small. To achieve that, we define the acceptance rate as the number of accepted steps over the total number of steps. Adjusting the time step such that the acceptance rate is close to 0.5 is a good compromise for the current problem.
#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a,nmax,dt):
# TODO
return energy/nmax, N_accep/nmax
# Run simulation
a = 1.2
nmax = 100000
dt = #TODO
X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
# Energy
X = [ x for (x, _) in X0 ]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
print(f"A = {A} +/- {deltaA}")subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
implicit none
double precision, intent(in) :: a
integer*8 , intent(in) :: nmax
double precision, intent(in) :: dt
double precision, intent(out) :: energy
double precision, intent(out) :: accep
integer*8 :: istep
integer*8 :: n_accep
double precision :: r_old(3), r_new(3), psi_old, psi_new
double precision :: v, ratio
double precision, external :: e_loc, psi, gaussian
! TODO
end subroutine metropolis_montecarlo
program qmc
implicit none
double precision, parameter :: a = 1.2d0
double precision, parameter :: dt = ! TODO
integer*8 , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns), Y(nruns)
double precision :: ave, err
do irun=1,nruns
call metropolis_montecarlo(a,nmax,dt,X(irun),Y(irun))
enddo
call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err
call ave_error(Y,nruns,ave,err)
print *, 'A = ', ave, '+/-', err
end program qmc#include <stdio.h>
#include <stdlib.h>
#include <stddef.h> // for size_t
#include <time.h>
#include <math.h>
#include "hydrogen.h"
#include "qmc_stats.h"
void metropolis_montecarlo(double a, size_t nmax, double dt,
double *energy, double *accep)
{
// TODO
}
int main(void) {
#define a 1.2
#define nmax 100000
#define dt //TODO
#define nruns 30
double energy[nruns];
double accep[nruns];
double ave, err;
srand48(time(NULL));
for (size_t irun = 0; irun < nruns; irun++) {
metropolis_montecarlo(a, nmax, dt, energy, accep);
}
ave_error(energy, nruns, &ave, &err);
printf("E = %f +/- %f\n", ave, err);
ave_error(accep, nruns, &ave, &err);
printf("A = %f +/- %f\n", ave, err);
return 0;
}#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a,nmax,dt):
energy = 0.
N_accep = 0
r_old = np.random.uniform(-dt, dt, (3))
psi_old = psi(a,r_old)
for istep in range(nmax):
energy += e_loc(a,r_old)
r_new = r_old + np.random.uniform(-dt,dt,(3))
psi_new = psi(a,r_new)
ratio = (psi_new / psi_old)**2
if np.random.uniform() <= ratio:
N_accep += 1
r_old = r_new
psi_old = psi_new
return energy/nmax, N_accep/nmax
# Run simulation
a = 1.2
nmax = 100000
dt = 1.0
X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
# Energy
X = [ x for (x, _) in X0 ]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
print(f"A = {A} +/- {deltaA}")E = -0.4802595860693983 +/- 0.0005124420418289207 A = 0.5074913333333334 +/- 0.000350889422714878
subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
implicit none
double precision, intent(in) :: a
integer*8 , intent(in) :: nmax
double precision, intent(in) :: dt
double precision, intent(out) :: energy
double precision, intent(out) :: accep
double precision :: r_old(3), r_new(3), psi_old, psi_new
double precision :: v, ratio, u(3)
integer*8 :: n_accep
integer*8 :: istep
double precision, external :: e_loc, psi, gaussian
energy = 0.d0
n_accep = 0_8
call random_number(r_old)
r_old(:) = dt * (2.d0*r_old(:) - 1.d0)
psi_old = psi(a,r_old)
do istep = 1,nmax
energy = energy + e_loc(a,r_old)
call random_number(u)
r_new(:) = r_old(:) + dt*(2.d0*u(:) - 1.d0)
psi_new = psi(a,r_new)
ratio = (psi_new / psi_old)**2
call random_number(v)
if (v <= ratio) then
n_accep = n_accep + 1_8
r_old(:) = r_new(:)
psi_old = psi_new
endif
end do
energy = energy / dble(nmax)
accep = dble(n_accep) / dble(nmax)
end subroutine metropolis_montecarlo
program qmc
implicit none
double precision, parameter :: a = 1.2d0
double precision, parameter :: dt = 1.0d0
integer*8 , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns), Y(nruns)
double precision :: ave, err
do irun=1,nruns
call metropolis_montecarlo(a,nmax,dt,X(irun),Y(irun))
enddo
call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err
call ave_error(Y,nruns,ave,err)
print *, 'A = ', ave, '+/-', err
end program qmcE = -0.48031509056922989 +/- 4.6803231523088529E-004 A = 0.50765333333333340 +/- 3.4368267017377720E-004
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h> // for size_t
#include <math.h>
#include <time.h>
#include "hydrogen.h"
#include "qmc_stats.h"
void metropolis_montecarlo(double a, size_t nmax, double dt,
double *energy, double *accep)
{
double r_old[3], r_new[3], psi_old, psi_new, v, ratio;
size_t n_accep = 0;
*energy = 0.0;
for (int i = 0; i < 3; i++) {
r_old[i] = dt * (2.0*drand48() - 1.0);
}
psi_old = psi(a, r_old);
for (size_t istep = 0; istep < nmax; istep++) {
*energy += e_loc(a, r_old);
for (int i = 0; i < 3; i++) {
r_new[i] = r_old[i] + dt * (2.0*drand48() - 1.0);
}
psi_new = psi(a, r_new);
ratio = pow(psi_new / psi_old,2);
v = drand48();
if (v <= ratio) {
n_accep++;
for (int i = 0; i < 3; i++) {
r_old[i] = r_new[i];
}
psi_old = psi_new;
}
}
*energy = *energy / (double) nmax;
*accep = (double) n_accep / (double) nmax;
}
int main(void) {
#define a 1.2
#define nmax 100000
#define dt 1.0
#define nruns 30
double X[nruns];
double Y[nruns];
double ave, err;
srand48(time(NULL));
for (size_t irun = 0; irun < nruns; irun++) {
metropolis_montecarlo(a, nmax, dt, &X[irun], &Y[irun]);
}
ave_error(X, nruns, &ave, &err);
printf("E = %f +/- %f\n", ave, err);
ave_error(Y, nruns, &ave, &err);
printf("A = %f +/- %f\n", ave, err);
return 0;
}E = -0.479518 +/- 0.000466 A = 0.507560 +/- 0.000353
One can use more efficient numerical schemes to move the electrons by choosing a smarter expression for the transition probability.
The Metropolis acceptance step has to be adapted keeping in mind that
the detailed balance condition is satisfied. This means that the acceptance
probability
\[ P(\mathbf{r}n → \mathbf{r}n+1) = A(\mathbf{r}n → \mathbf{r}n+1) T(\mathbf{r}n → \mathbf{r}n+1) = A(\mathbf{r}n+1 → \mathbf{r}n) T(\mathbf{r}n+1 → \mathbf{r}n) \frac{P(\mathbf{r}n+1)}{P(\mathbf{r}n)} \]
where $T(\mathbf{r}_n → \mathbf{r}n+1)$ is the
probability of transition from
In the previous example, we were using uniform sampling in a box centered at the current position. Hence, the transition probability was symmetric
\[ T(\mathbf{r}n → \mathbf{r}n+1) = T(\mathbf{r}n+1 → \mathbf{r}n) = \text{constant}\,, \]
so the expression of
Now, if instead of drawing uniform random numbers, we
choose to draw Gaussian random numbers with zero mean and variance
\[ T(\mathbf{r}n → \mathbf{r}n+1) = \frac{1}{(2π\,δ t)3/2} exp \left[ - \frac{\left( \mathbf{r}n+1 - \mathbf{r}n \right)^2}{2δ t} \right]\,. \]
Furthermore, to sample the density even better, we can “push” the electrons into in the regions of high probability, and “pull” them away from the low-probability regions. This will increase the acceptance ratios and improve the sampling.
To do this, we can use the gradient of the probability density
\[ \frac{∇ [ Ψ^2 ]}{Ψ^2} = 2 \frac{∇ Ψ}{Ψ}\,, \]
and add the so-called drift vector,
\[ T(\mathbf{r}n → \mathbf{r}n+1) = \frac{1}{(2π\,δ t)3/2} exp \left[ - \frac{\left( \mathbf{r}n+1 - \mathbf{r}n - δ t\frac{∇ Ψ(\mathbf{r}_n)}{Ψ(\mathbf{r}_n)} \right)^2}{2\,δ t} \right]\,. \]
The corresponding move is proposed as
\[ \mathbf{r}n+1 = \mathbf{r}n + δ t\, \frac{∇ Ψ(\mathbf{r})}{Ψ(\mathbf{r})} + \sqrt{δ t}\, χ \,, \]
where
Here is an illustration of a trajectory:
Averaging all the trajectories converges to the density of the
hydrogen atom.
The algorithm of the previous exercise is only slightly modified as:
- Evaluate the local energy at $\mathbf{r}n$ and accumulate it
- Compute a new position
$\mathbf{r’} = \mathbf{r}_n + δ t\, \frac{∇ Ψ(\mathbf{r})}{Ψ(\mathbf{r})} + \sqrt{δ t}\,χ$ - Evaluate
$Ψ(\mathbf{r}’)$ and$\frac{∇ Ψ(\mathbf{r’})}{Ψ(\mathbf{r’})}$ at the new position - Compute the ratio $A = \frac{T(\mathbf{r}’ → \mathbf{r}n) P(\mathbf{r}’)}{T(\mathbf{r}n → \mathbf{r}’) P(\mathbf{r}n)}$
- Draw a uniform random number
$v ∈ [0,1]$ - if
$v ≤ A$ , accept the move : set $\mathbf{r}n+1 = \mathbf{r’}$ - else, reject the move : set $\mathbf{r}n+1 = \mathbf{r}_n$
To obtain Gaussian-distributed random numbers, you can apply the Box Muller transform to uniform random numbers:
\begin{eqnarray*}
z_1 &=& \sqrt{-2 ln u_1} cos(2 π u_2)
z_2 &=& \sqrt{-2 ln u_1} sin(2 π u_2)
\end{eqnarray*}
This gives Gaussian-distributed random numbers with zero mean and a variance of one.
Below is a Fortran and a C implementation returning a
Gaussian-distributed n-dimensional vector
subroutine random_gauss(z,n)
implicit none
integer, intent(in) :: n
double precision, intent(out) :: z(n)
double precision :: u(n+1)
double precision, parameter :: two_pi = 2.d0*dacos(-1.d0)
integer :: i
call random_number(u)
if (iand(n,1) == 0) then
! n is even
do i=1,n,2
z(i) = dsqrt(-2.d0*dlog(u(i)))
z(i+1) = z(i) * dsin( two_pi*u(i+1) )
z(i) = z(i) * dcos( two_pi*u(i+1) )
end do
else
! n is odd
do i=1,n-1,2
z(i) = dsqrt(-2.d0*dlog(u(i)))
z(i+1) = z(i) * dsin( two_pi*u(i+1) )
z(i) = z(i) * dcos( two_pi*u(i+1) )
end do
z(n) = dsqrt(-2.d0*dlog(u(n)))
z(n) = z(n) * dcos( two_pi*u(n+1) )
end if
end subroutine random_gaussvoid random_gauss(double* z, size_t n);#include <stdlib.h>
void random_gauss(double* z, size_t n) {
double u[n+1];
double two_pi = 2.0 * acos(-1.0);
size_t i;
//generate random numbers
for (i = 0; i <= n; i++) {
u[i] = drand48();
}
if (n % 2 == 0) {
// n is even
for (i = 0; i < n; i += 2) {
z[i] = sqrt(-2.0 * log(u[i]));
z[i+1] = z[i] * sin(two_pi * u[i+1]);
z[i] = z[i] * cos(two_pi * u[i+1]);
}
} else {
// n is odd
for (i = 0; i < n-1; i += 2) {
z[i] = sqrt(-2.0 * log(u[i]));
z[i+1] = z[i] * sin(two_pi * u[i+1]);
z[i] = z[i] * cos(two_pi * u[i+1]);
}
z[n-1] = sqrt(-2.0 * log(u[n-1]));
z[n-1] = z[n-1] * cos(two_pi * u[n]);
}
}def drift(a,r):
# TODOsubroutine drift(a,r,b)
implicit none
double precision, intent(in) :: a, r(3)
double precision, intent(out) :: b(3)
! TODO
end subroutine driftvoid drift(double a, double r[3], double b[3]) {
//TODO
}def drift(a,r):
ar_inv = -a/np.sqrt(np.dot(r,r))
return r * ar_invsubroutine drift(a,r,b)
implicit none
double precision, intent(in) :: a, r(3)
double precision, intent(out) :: b(3)
double precision :: ar_inv
ar_inv = -a / dsqrt(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
b(:) = r(:) * ar_inv
end subroutine driftvoid drift(double a, double r[3], double b[3]);void drift(double a, double r[3], double b[3]) {
double ar_inv = -a / sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
for (int i = 0; i < 3; i++) {
b[i] = r[i] * ar_inv;
}
}#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a,nmax,dt):
# TODO
# Run simulation
a = 1.2
nmax = 100000
dt = # TODO
X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
# Energy
X = [ x for (x, _) in X0 ]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
print(f"A = {A} +/- {deltaA}")subroutine variational_montecarlo(a,nmax,dt,energy,accep)
implicit none
double precision, intent(in) :: a, dt
integer*8 , intent(in) :: nmax
double precision, intent(out) :: energy, accep
integer*8 :: istep
integer*8 :: n_accep
double precision :: sq_dt, chi(3)
double precision :: psi_old, psi_new
double precision :: r_old(3), r_new(3)
double precision :: d_old(3), d_new(3)
double precision, external :: e_loc, psi
! TODO
end subroutine variational_montecarlo
program qmc
implicit none
double precision, parameter :: a = 1.2d0
double precision, parameter :: dt = ! TODO
integer*8 , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns), accep(nruns)
double precision :: ave, err
do irun=1,nruns
call variational_montecarlo(a,nmax,dt,X(irun),accep(irun))
enddo
call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err
call ave_error(accep,nruns,ave,err)
print *, 'A = ', ave, '+/-', err
end program qmcgfortran hydrogen.F90 qmc_stats.F90 vmc_metropolis.F90 -o vmc_metropolis
./vmc_metropolis#include <stdio.h>
#include <stdlib.h>
#include <stddef.h> // for size_t
#include <math.h>
#include <time.h>
#include "hydrogen.h"
#include "qmc_stats.h"
void variational_montecarlo(double a, size_t nmax, double dt,
double *energy, double *accep)
{
// TODO
}
int main(void) {
#define a 1.2
#define nmax 100000
#define dt // TODO
#define nruns 30
double X[nruns];
double Y[nruns];
double ave, err;
srand48(time(NULL));
for (size_t irun = 0; irun < nruns; irun++) {
variational_montecarlo(a, nmax, dt, &X[irun], &Y[irun]);
}
ave_error(X, nruns, &ave, &err);
printf("E = %f +/- %f\n", ave, err);
ave_error(Y, nruns, &ave, &err);
printf("A = %f +/- %f\n", ave, err);
return 0;
}#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a,nmax,dt):
sq_dt = np.sqrt(dt)
energy = 0.
N_accep = 0
r_old = np.random.normal(loc=0., scale=1.0, size=(3))
d_old = drift(a,r_old)
d2_old = np.dot(d_old,d_old)
psi_old = psi(a,r_old)
for istep in range(nmax):
chi = np.random.normal(loc=0., scale=1.0, size=(3))
energy += e_loc(a,r_old)
r_new = r_old + dt * d_old + sq_dt * chi
d_new = drift(a,r_new)
d2_new = np.dot(d_new,d_new)
psi_new = psi(a,r_new)
# Metropolis
prod = np.dot((d_new + d_old), (r_new - r_old))
argexpo = 0.5 * (d2_new - d2_old)*dt + prod
q = psi_new / psi_old
q = np.exp(-argexpo) * q*q
if np.random.uniform() <= q:
N_accep += 1
r_old = r_new
d_old = d_new
d2_old = d2_new
psi_old = psi_new
return energy/nmax, N_accep/nmax
# Run simulation
a = 1.2
nmax = 100000
dt = 1.0
X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
# Energy
X = [ x for (x, _) in X0 ]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
print(f"A = {A} +/- {deltaA}")E = -0.48034171558629885 +/- 0.0005286038561061781 A = 0.6210380000000001 +/- 0.0005457375900937905
subroutine variational_montecarlo(a,nmax,dt,energy,accep)
implicit none
double precision, intent(in) :: a, dt
integer*8 , intent(in) :: nmax
double precision, intent(out) :: energy, accep
integer*8 :: istep
integer*8 :: n_accep
double precision :: sq_dt, chi(3), d2_old, prod, u
double precision :: psi_old, psi_new, d2_new, argexpo, q
double precision :: r_old(3), r_new(3)
double precision :: d_old(3), d_new(3)
double precision, external :: e_loc, psi
sq_dt = dsqrt(dt)
! Initialization
energy = 0.d0
n_accep = 0_8
call random_gauss(r_old,3)
call drift(a,r_old,d_old)
d2_old = d_old(1)*d_old(1) + &
d_old(2)*d_old(2) + &
d_old(3)*d_old(3)
psi_old = psi(a,r_old)
do istep = 1,nmax
energy = energy + e_loc(a,r_old)
call random_gauss(chi,3)
r_new(:) = r_old(:) + dt*d_old(:) + chi(:)*sq_dt
call drift(a,r_new,d_new)
d2_new = d_new(1)*d_new(1) + &
d_new(2)*d_new(2) + &
d_new(3)*d_new(3)
psi_new = psi(a,r_new)
! Metropolis
prod = (d_new(1) + d_old(1))*(r_new(1) - r_old(1)) + &
(d_new(2) + d_old(2))*(r_new(2) - r_old(2)) + &
(d_new(3) + d_old(3))*(r_new(3) - r_old(3))
argexpo = 0.5d0 * (d2_new - d2_old)*dt + prod
q = psi_new / psi_old
q = dexp(-argexpo) * q*q
call random_number(u)
if (u <= q) then
n_accep = n_accep + 1_8
r_old(:) = r_new(:)
d_old(:) = d_new(:)
d2_old = d2_new
psi_old = psi_new
end if
end do
energy = energy / dble(nmax)
accep = dble(n_accep) / dble(nmax)
end subroutine variational_montecarlo
program qmc
implicit none
double precision, parameter :: a = 1.2d0
double precision, parameter :: dt = 1.0d0
integer*8 , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns), accep(nruns)
double precision :: ave, err
do irun=1,nruns
call variational_montecarlo(a,nmax,dt,X(irun),accep(irun))
enddo
call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err
call ave_error(accep,nruns,ave,err)
print *, 'A = ', ave, '+/-', err
end program qmcE = -0.47965766050774800 +/- 5.0599992102209337E-004 A = 0.62055466666666670 +/- 4.3100330247376571E-004
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h> // for size_t
#include <math.h>
#include <time.h>
#include "hydrogen.h"
#include "qmc_stats.h"
void variational_montecarlo(double a, size_t nmax, double dt,
double *energy, double *accep)
{
double chi[3], d2_old, prod, u;
double psi_old, psi_new, d2_new, argexpo, q;
double r_old[3], r_new[3];
double d_old[3], d_new[3];
size_t istep, n_accep = 0;
double sq_dt = sqrt(dt);
// Initialization
*energy = 0.0;
random_gauss(r_old, 3);
drift(a, r_old, d_old);
d2_old = d_old[0]*d_old[0] + d_old[1]*d_old[1] + d_old[2]*d_old[2];
psi_old = psi(a, r_old);
for (istep = 0; istep < nmax; istep++) {
*energy += e_loc(a, r_old);
random_gauss(chi, 3);
for (int i = 0; i < 3; i++) {
r_new[i] = r_old[i] + dt*d_old[i] + chi[i]*sq_dt;
}
drift(a, r_new, d_new);
d2_new = d_new[0]*d_new[0] + d_new[1]*d_new[1] + d_new[2]*d_new[2];
psi_new = psi(a, r_new);
// Metropolis
prod = (d_new[0] + d_old[0])*(r_new[0] - r_old[0]) +
(d_new[1] + d_old[1])*(r_new[1] - r_old[1]) +
(d_new[2] + d_old[2])*(r_new[2] - r_old[2]);
argexpo = 0.5 * (d2_new - d2_old)*dt + prod;
q = psi_new / psi_old;
q = exp(-argexpo) * q*q;
u = drand48();
if (u <= q) {
n_accep++;
for (int i = 0; i < 3; i++) {
r_old[i] = r_new[i];
d_old[i] = d_new[i];
}
d2_old = d2_new;
psi_old = psi_new;
}
}
*energy = *energy / (double) nmax;
*accep = (double) n_accep / (double) nmax;
}
int main(void) {
#define a 1.2
#define nmax 100000
#define dt 1.0
#define nruns 30
double X[nruns];
double Y[nruns];
double ave, err;
srand48(time(NULL));
for (size_t irun = 0; irun < nruns; irun++) {
variational_montecarlo(a, nmax, dt, &X[irun], &Y[irun]);
}
ave_error(X, nruns, &ave, &err);
printf("E = %f +/- %f\n", ave, err);
ave_error(Y, nruns, &ave, &err);
printf("A = %f +/- %f\n", ave, err);
return 0;
}E = -0.479814 +/- 0.000552 A = 0.621853 +/- 0.000401
As we have seen, Variational Monte Carlo is a powerful method to compute integrals in large dimensions. It is often used in cases where the expression of the wave function is such that the integrals can’t be evaluated (multi-centered Slater-type orbitals, correlation factors, etc).
Diffusion Monte Carlo is different. It goes beyond the computation of the integrals associated with an input wave function, and aims at finding a near-exact numerical solution to the Schrödinger equation.
Consider the time-dependent Schrödinger equation:
\[ i\frac{∂ Ψ(\mathbf{r},t)}{∂ t} = (\hat{H} -E\rm ref) Ψ(\mathbf{r},t)\,. \]
where we introduced a shift in the energy, $E\rm ref$, for reasons which will become apparent below.
We can expand a given starting wave function,
\[ Ψ(\mathbf{r},0) = ∑_k a_k\, Φ_k(\mathbf{r}). \]
The solution of the Schrödinger equation at time
\[ Ψ(\mathbf{r},t) = ∑_k a_k exp \left( -i\, (E_k-E\rm ref)\, t \right) Φ_k(\mathbf{r}). \]
Now, if we replace the time variable
\[ -\frac{∂ ψ(\mathbf{r}, τ)}{∂ τ} = (\hat{H} -E\rm ref) ψ(\mathbf{r}, τ) \]
where
\begin{eqnarray*}
ψ(\mathbf{r},τ) &=& ∑_k a_k exp( -(E_k-E\rm ref)\, τ) Φ_k(\mathbf{r})
&=& exp(-(E_0-E\rm ref)\, τ)∑_k a_k exp( -(E_k-E_0)\, τ) Φ_k(\mathbf{r})\,.
\end{eqnarray*}
For large positive values of
The diffusion equation of particles is given by
\[ \frac{∂ ψ(\mathbf{r},t)}{∂ t} = D\, Δ ψ(\mathbf{r},t) \]
where
\[ \frac{∂ ψ(\mathbf{r}, τ)}{∂ τ} = \left(\frac{1}{2}Δ - [V(\mathbf{r}) -E\rm ref]\right) ψ(\mathbf{r}, τ)\,, \]
it can be identified as the combination of:
- a diffusion equation (Laplacian)
- an equation whose solution is an exponential (potential)
The diffusion equation can be simulated by a Brownian motion:
\[ \mathbf{r}n+1 = \mathbf{r}n + \sqrt{δ t}\, χ \]
where
\[ ∏_i exp \left( - (V(\mathbf{r}_i) - E\text{ref}) δ t \right). \]
We note that the ground-state wave function of a Fermionic system is
antisymmetric and changes sign. Therefore, its interpretation as a probability
distribution is somewhat problematic. In fact, mathematically, since
the Bosonic ground state is lower in energy than the Fermionic one, for
large
For the systems you will study, this is not an issue:
- Hydrogen atom: You only have one electron!
- Two-electron system (
$H_2$ or He): The ground-wave function is antisymmetric in the spin variables but symmetric in the space ones.
Therefore, in both cases, you are dealing with a “Bosonic” ground state.
In a molecular system, the potential is far from being constant
and, in fact, diverges at the inter-particle coalescence points. Hence,
it results in very large fluctuations of the weight associated with
the potential, making the calculations impossible in practice.
Fortunately, if we multiply the Schrödinger equation by a chosen
trial wave function
\[ -\frac{∂ ψ(\mathbf{r},τ)}{∂ τ} Ψ_T(\mathbf{r}) = \left[ -\frac{1}{2} Δ ψ(\mathbf{r},τ) + V(\mathbf{r}) ψ(\mathbf{r},τ) \right] Ψ_T(\mathbf{r}) \]
Defining
\[ -\frac{∂ Π(\mathbf{r},τ)}{∂ τ} = -\frac{1}{2} Δ Π(\mathbf{r},τ) + ∇ \left[ Π(\mathbf{r},τ) \frac{∇ Ψ_T(\mathbf{r})}{Ψ_T(\mathbf{r})} \right] + (E_L(\mathbf{r})-E\rm ref)Π(\mathbf{r},τ) \]
The new “kinetic energy” can be simulated by the drift-diffusion
scheme presented in the previous section (VMC).
The new “potential” is the local energy, which has smaller fluctuations
when
This equation generates the N-electron density
To this aim, we use the mixed estimator of the energy:
\begin{eqnarray*}
E(τ) &=& \frac{\langle ψ(τ) | \hat{H} | Ψ_T \rangle}{\langle ψ(τ) | Ψ_T \rangle}
&=& \frac{∫ ψ(\mathbf{r},τ) \hat{H} Ψ_T(\mathbf{r}) d\mathbf{r}}
{∫ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) d\mathbf{r}} \
&=& \frac{∫ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) E_L(\mathbf{r}) d\mathbf{r}}
{∫ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) d\mathbf{r}} \,.
\end{eqnarray*}
For large
\[ Π(\mathbf{r},τ) =ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) → Φ_0(\mathbf{r}) Ψ_T(\mathbf{r})\,, \]
and, using that
\[ E(τ) = \frac{\langle ψ_τ | \hat{H} | Ψ_T \rangle} {\langle ψ_τ | Ψ_T \rangle} = \frac{\langle Ψ_T | \hat{H} | ψ_τ \rangle} {\langle Ψ_T | ψ_τ \rangle} → E_0 \frac{\langle Ψ_T | Φ_0 \rangle} {\langle Ψ_T | Φ_0 \rangle} = E_0 \]
Therefore, we can compute the energy within DMC by generating the
density
\[ -\frac{∂ ψ(\mathbf{r},τ)}{∂ τ} Ψ_T(\mathbf{r}) = \left[ -\frac{1}{2} Δ ψ(\mathbf{r},τ) + V(\mathbf{r}) ψ(\mathbf{r},τ) \right] Ψ_T(\mathbf{r}) \]
\[ -\frac{∂ \big[ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big]}{∂ τ} = -\frac{1}{2} \Big( Δ \big[ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big] - ψ(\mathbf{r},τ) Δ Ψ_T(\mathbf{r}) - 2 ∇ ψ(\mathbf{r},τ) ∇ Ψ_T(\mathbf{r}) \Big) + V(\mathbf{r}) ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \]
\[ -\frac{∂ \big[ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big]}{∂ τ} = -\frac{1}{2} Δ \big[ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big] + \frac{1}{2} ψ(\mathbf{r},τ) Δ Ψ_T(\mathbf{r}) + Ψ_T(\mathbf{r})∇ ψ(\mathbf{r},τ) \frac{∇ Ψ_T(\mathbf{r})}{Ψ_T(\mathbf{r})} + V(\mathbf{r}) ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \]
\[ -\frac{∂ \big[ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big]}{∂ τ} = -\frac{1}{2} Δ \big[ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big] + ψ(\mathbf{r},τ) Δ Ψ_T(\mathbf{r}) + Ψ_T(\mathbf{r})∇ ψ(\mathbf{r},τ) \frac{∇ Ψ_T(\mathbf{r})}{Ψ_T(\mathbf{r})} + E_L(\mathbf{r}) ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \] \[ -\frac{∂ \big[ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big]}{∂ τ} = -\frac{1}{2} Δ \big[ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \big] + ∇ \left[ ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \frac{∇ Ψ_T(\mathbf{r})}{Ψ_T(\mathbf{r})} \right] + E_L(\mathbf{r}) ψ(\mathbf{r},τ) Ψ_T(\mathbf{r}) \]
Defining
\[ -\frac{∂ Π(\mathbf{r},τ)}{∂ τ} = -\frac{1}{2} Δ Π(\mathbf{r},τ) + ∇ \left[ Π(\mathbf{r},τ) \frac{∇ Ψ_T(\mathbf{r})}{Ψ_T(\mathbf{r})} \right] + E_L(\mathbf{r}) Π(\mathbf{r},τ) \]
Instead of having a variable number of particles to simulate the branching process as in the Diffusion Monte Carlo (DMC) algorithm, we use variant called pure Diffusion Monte Carlo (PDMC) where the potential term is considered as a cumulative product of weights:
\begin{eqnarray*} W(\mathbf{r}_n, τ) = ∏i=1n exp \left( -δ t\, (E_L(\mathbf{r}_i) - E\text{ref}) \right) = ∏i=1n w(\mathbf{r}_i) \end{eqnarray*}
where
The PDMC algorithm is less stable than the DMC algorithm: it
requires to have a value of
- Start with
$W(\mathbf{r}_0)=1, τ_0 = 0$ - Evaluate the local energy at $\mathbf{r}n$
- Compute the contribution to the weight
$w(\mathbf{r}_n) = exp(-δ t(E_L(\mathbf{r}_n)-E_\text{ref}))$ - Update $W(\mathbf{r}n) = W(\mathbf{r}n-1) × w(\mathbf{r}_n)$
- Accumulate the weighted energy
$W(\mathbf{r}_n) × E_L(\mathbf{r}_n)$ , and the weight$W(\mathbf{r}_n)$ for the normalization - Update $τ_n = τn-1 + δ t$
- If $τn > τ_\text{max}$ (
$τ_\text{max}$ is an input parameter), the long projection time has been reached and we can start an new trajectory from the current position: reset$W(r_n) = 1$ and$τ_n = 0$ - Compute a new position
$\mathbf{r’} = \mathbf{r}_n + δ t\, \frac{∇ Ψ(\mathbf{r})}{Ψ(\mathbf{r})} + \sqrt{δ t}\, χ$ - Evaluate
$Ψ(\mathbf{r}’)$ and$\frac{∇ Ψ(\mathbf{r’})}{Ψ(\mathbf{r’})}$ at the new position - Compute the ratio $A = \frac{T(\mathbf{r}’ → \mathbf{r}n) P(\mathbf{r}’)}{T(\mathbf{r}n → \mathbf{r}’) P(\mathbf{r}n)}$
- Draw a uniform random number
$v ∈ [0,1]$ - if
$v ≤ A$ , accept the move : set $\mathbf{r}n+1 = \mathbf{r’}$ - else, reject the move : set $\mathbf{r}n+1 = \mathbf{r}_n$
Some comments are needed:
- You estimate the energy as
\begin{eqnarray*} E = \frac{∑k=1N_{\rm MC} E_L(\mathbf{r}_k) W(\mathbf{r}_k, kδ t)}{∑k=1N_{\rm MC} W(\mathbf{r}_k, kδ t)} \end{eqnarray*}
- The result will be affected by a time-step error
(the finite size of
$δ t$ ) due to the discretization of the integral, and one has in principle to extrapolate to the limit$δ t → 0$ . This amounts to fitting the energy computed for multiple values of$δ t$ .Here, you will be using a small enough time-step and you should not worry about the extrapolation.
- The accept/reject step (steps 9-12 in the algorithm) is in principle not needed for the correctness of the DMC algorithm. However, its use reduces significantly the time-step error.
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a, nmax, dt, Eref):
# TODO
# Run simulation
a = 1.2
nmax = 100000
dt = 0.05
tau = 100.
E_ref = -0.5
X0 = [ MonteCarlo(a, nmax, dt, E_ref) for i in range(30)]
# Energy
X = [ x for (x, _) in X0 ]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
print(f"A = {A} +/- {deltaA}")subroutine pdmc(a, nmax, dt, energy, accep, tau, E_ref)
implicit none
double precision, intent(in) :: a, dt, tau
integer*8 , intent(in) :: nmax
double precision, intent(out) :: energy, accep
double precision, intent(in) :: E_ref
integer*8 :: istep
integer*8 :: n_accep
double precision :: sq_dt, chi(3)
double precision :: psi_old, psi_new
double precision :: r_old(3), r_new(3)
double precision :: d_old(3), d_new(3)
double precision, external :: e_loc, psi
! TODO
end subroutine pdmc
program qmc
implicit none
double precision, parameter :: a = 1.2d0
double precision, parameter :: dt = 0.05d0
double precision, parameter :: E_ref = -0.5d0
double precision, parameter :: tau = 100.d0
integer*8 , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns), accep(nruns)
double precision :: ave, err
do irun=1,nruns
call pdmc(a, nmax, dt, X(irun), accep(irun), tau, E_ref)
enddo
call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err
call ave_error(accep,nruns,ave,err)
print *, 'A = ', ave, '+/-', err
end program qmcgfortran hydrogen.F90 qmc_stats.F90 pdmc.F90 -o pdmc
./pdmc#include <stdio.h>
#include <stdlib.h>
#include <stddef.h> // for size_t
#include <math.h>
#include <time.h>
#include "hydrogen.h"
#include "qmc_stats.h"
void pdmc(double a, size_t nmax, double dt, double *energy, double *accep,
double tau, double e_ref)
{
// TODO
}
int main(void) {
#define a 1.2
#define nmax 100000
#define dt 0.05
#define tau 100.0
#define e_ref -0.5
#define nruns 30
double X[nruns];
double Y[nruns];
double ave, err;
srand48(time(NULL));
for (size_t irun = 0; irun < nruns; irun++) {
pdmc(a, nmax, dt, &X[irun], &Y[irun], tau, e_ref);
}
ave_error(X, nruns, &ave, &err);
printf("E = %f +/- %f\n", ave, err);
ave_error(Y, nruns, &ave, &err);
printf("A = %f +/- %f\n", ave, err);
return 0;
}gcc hydrogen.c qmc_stats.c pdmc.c -lm -o pdmc
./pdmc#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a, nmax, dt, tau, Eref):
sq_dt = np.sqrt(dt)
energy = 0.
N_accep = 0
normalization = 0.
w = 1.
tau_current = 0.
r_old = np.random.normal(loc=0., scale=1.0, size=(3))
d_old = drift(a,r_old)
d2_old = np.dot(d_old,d_old)
psi_old = psi(a,r_old)
for istep in range(nmax):
el = e_loc(a,r_old)
w *= np.exp(-dt*(el - Eref))
normalization += w
energy += w * el
tau_current += dt
# Reset when tau is reached
if tau_current > tau:
w = 1.
tau_current = 0.
chi = np.random.normal(loc=0., scale=1.0, size=(3))
r_new = r_old + dt * d_old + sq_dt * chi
d_new = drift(a,r_new)
d2_new = np.dot(d_new,d_new)
psi_new = psi(a,r_new)
# Metropolis
prod = np.dot((d_new + d_old), (r_new - r_old))
argexpo = 0.5 * (d2_new - d2_old)*dt + prod
q = psi_new / psi_old
q = np.exp(-argexpo) * q*q
if np.random.uniform() <= q:
N_accep += 1
r_old = r_new
d_old = d_new
d2_old = d2_new
psi_old = psi_new
return energy/normalization, N_accep/nmax
# Run simulation
a = 1.2
nmax = 100000
dt = 0.05
tau = 100.
E_ref = -0.5
X0 = [ MonteCarlo(a, nmax, dt, tau, E_ref) for i in range(30)]
# Energy
X = [ x for (x, _) in X0 ]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
print(f"A = {A} +/- {deltaA}")subroutine pdmc(a, dt, nmax, energy, accep, tau, E_ref)
implicit none
double precision, intent(in) :: a, dt, tau
integer*8 , intent(in) :: nmax
double precision, intent(out) :: energy, accep
double precision, intent(in) :: E_ref
integer*8 :: istep
integer*8 :: n_accep
double precision :: sq_dt, chi(3), d2_old, prod, u
double precision :: psi_old, psi_new, d2_new, argexpo, q
double precision :: r_old(3), r_new(3)
double precision :: d_old(3), d_new(3)
double precision :: e, w, normalization, tau_current
double precision, external :: e_loc, psi
sq_dt = dsqrt(dt)
! Initialization
energy = 0.d0
n_accep = 0_8
normalization = 0.d0
w = 1.d0
tau_current = 0.d0
call random_gauss(r_old,3)
call drift(a,r_old,d_old)
d2_old = d_old(1)*d_old(1) + &
d_old(2)*d_old(2) + &
d_old(3)*d_old(3)
psi_old = psi(a,r_old)
do istep = 1,nmax
e = e_loc(a,r_old)
w = w * dexp(-dt*(e - E_ref))
normalization = normalization + w
energy = energy + w*e
tau_current = tau_current + dt
! Reset when tau is reached
if (tau_current > tau) then
w = 1.d0
tau_current = 0.d0
endif
call random_gauss(chi,3)
r_new(:) = r_old(:) + dt*d_old(:) + chi(:)*sq_dt
call drift(a,r_new,d_new)
d2_new = d_new(1)*d_new(1) + &
d_new(2)*d_new(2) + &
d_new(3)*d_new(3)
psi_new = psi(a,r_new)
! Metropolis
prod = (d_new(1) + d_old(1))*(r_new(1) - r_old(1)) + &
(d_new(2) + d_old(2))*(r_new(2) - r_old(2)) + &
(d_new(3) + d_old(3))*(r_new(3) - r_old(3))
argexpo = 0.5d0 * (d2_new - d2_old)*dt + prod
q = psi_new / psi_old
q = dexp(-argexpo) * q*q
call random_number(u)
if (u <= q) then
n_accep = n_accep + 1_8
r_old(:) = r_new(:)
d_old(:) = d_new(:)
d2_old = d2_new
psi_old = psi_new
end if
end do
energy = energy / normalization
accep = dble(n_accep) / dble(nmax)
end subroutine pdmc
program qmc
implicit none
double precision, parameter :: a = 1.2d0
double precision, parameter :: dt = 0.05d0
double precision, parameter :: E_ref = -0.5d0
double precision, parameter :: tau = 100.d0
integer*8 , parameter :: nmax = 100000
integer , parameter :: nruns = 30
integer :: irun
double precision :: X(nruns), accep(nruns)
double precision :: ave, err
do irun=1,nruns
call pdmc(a, dt, nmax, X(irun), accep(irun), tau, E_ref)
enddo
call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err
call ave_error(accep,nruns,ave,err)
print *, 'A = ', ave, '+/-', err
end program qmcgfortran hydrogen.F90 qmc_stats.F90 pdmc.F90 -o pdmc
./pdmc#include <stdio.h>
#include <stdlib.h>
#include <stddef.h> // for size_t
#include <math.h>
#include <time.h>
#include "hydrogen.h"
#include "qmc_stats.h"
void pdmc(double a, size_t nmax, double dt, double *energy, double *accep,
double tau, double e_ref)
{
double chi[3], d2_old, prod, u;
double psi_old, psi_new, d2_new, argexpo, q;
double r_old[3], r_new[3];
double d_old[3], d_new[3];
size_t istep, n_accep;
double e, w, tau_current, normalization;
double sq_dt = sqrt(dt);
// Initialization
*energy = 0.0;
n_accep = 0;
normalization = 0.0;
w = 1.0;
tau_current = 0.0;
random_gauss(r_old, 3);
drift(a, r_old, d_old);
d2_old = d_old[0]*d_old[0] +
d_old[1]*d_old[1] +
d_old[2]*d_old[2];
psi_old = psi(a, r_old);
for (istep = 0; istep < nmax; istep++) {
e = e_loc(a, r_old);
w *= exp(-dt*(e-e_ref));
normalization += w;
*energy += w*e;
tau_current += dt;
// Reset when tau is reached
if (tau_current > tau) {
w = 1.0;
tau_current = 0.0;
}
random_gauss(chi, 3);
for (int i = 0; i < 3; i++) {
r_new[i] = r_old[i] + dt*d_old[i] + chi[i]*sq_dt;
}
drift(a, r_new, d_new);
d2_new = d_new[0]*d_new[0] +
d_new[1]*d_new[1] +
d_new[2]*d_new[2];
psi_new = psi(a, r_new);
// Metropolis
prod = (d_new[0] + d_old[0])*(r_new[0] - r_old[0]) +
(d_new[1] + d_old[1])*(r_new[1] - r_old[1]) +
(d_new[2] + d_old[2])*(r_new[2] - r_old[2]);
argexpo = 0.5 * (d2_new - d2_old)*dt + prod;
q = psi_new / psi_old;
q = exp(-argexpo) * q*q;
u = drand48();
if (u <= q) {
n_accep++;
for (int i = 0; i < 3; i++) {
r_old[i] = r_new[i];
d_old[i] = d_new[i];
}
d2_old = d2_new;
psi_old = psi_new;
}
}
*energy = *energy / normalization;
*accep = (double) n_accep / (double) nmax;
}
int main(void) {
#define a 1.2
#define dt 0.05
#define e_ref -0.5
#define tau 100.0
#define nmax 100000
#define nruns 30
double X[nruns];
double Y[nruns];
double ave, err;
srand48(time(NULL));
for (size_t irun = 0; irun < nruns; irun++) {
pdmc(a, nmax, dt, &X[irun], &Y[irun], tau, e_ref);
}
ave_error(X, nruns, &ave, &err);
printf("E = %f +/- %f\n", ave, err);
ave_error(Y, nruns, &ave, &err);
printf("A = %f +/- %f\n", ave, err);
return 0;
}gcc hydrogen.c qmc_stats.c pdmc.c -lm -o pdmc
./pdmcChange your PDMC code for one of the following:
- the Helium atom, or
- the H_2 molecule at
$R$ =1.401 bohr.
And compute the ground state energy.
TREX : Targeting Real Chemical Accuracy at the Exascale project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content.