mirror of
https://github.com/TREX-CoE/qmc-lttc.git
synced 2024-12-21 11:53:58 +01:00
Cleaning
This commit is contained in:
parent
0140aff11c
commit
707338aa8a
343
QMC.org
343
QMC.org
@ -2758,350 +2758,7 @@ gfortran hydrogen.f90 qmc_stats.f90 pdmc.f90 -o pdmc
|
|||||||
: A = 0.98963533333333342 +/- 6.3052128284666221E-005
|
: A = 0.98963533333333342 +/- 6.3052128284666221E-005
|
||||||
|
|
||||||
|
|
||||||
** H_2
|
|
||||||
|
|
||||||
We will now consider the H_2 molecule in a minimal basis composed of the
|
|
||||||
$1s$ orbitals of the hydrogen atoms:
|
|
||||||
|
|
||||||
$$
|
|
||||||
\Psi(\mathbf{r}_1, \mathbf{r}_2) =
|
|
||||||
\exp(-(\mathbf{r}_1 - \mathbf{R}_A)) +
|
|
||||||
$$
|
|
||||||
where $\mathbf{r}_1$ and $\mathbf{r}_2$ denote the electron
|
|
||||||
coordinates and $\mathbf{R}_A$ and $\mathbf{R}_B$ the coordinates of
|
|
||||||
the nuclei.
|
|
||||||
|
|
||||||
|
|
||||||
* Old sections to be removed :noexport:
|
|
||||||
|
|
||||||
:PROPERTIES:
|
|
||||||
:header-args:python: :tangle none
|
|
||||||
:header-args:f90: :tangle none
|
|
||||||
:END:
|
|
||||||
|
|
||||||
** Gaussian sampling :noexport:
|
|
||||||
:PROPERTIES:
|
|
||||||
:header-args:python: :tangle qmc_gaussian.py
|
|
||||||
:header-args:f90: :tangle qmc_gaussian.f90
|
|
||||||
:END:
|
|
||||||
|
|
||||||
We will now improve the sampling and allow to sample in the whole
|
|
||||||
3D space, correcting the bias related to the sampling in the box.
|
|
||||||
|
|
||||||
Instead of drawing uniform random numbers, we will draw Gaussian
|
|
||||||
random numbers centered on 0 and with a variance of 1.
|
|
||||||
|
|
||||||
Now the sampling probability can be inserted into the equation of the energy:
|
|
||||||
|
|
||||||
\[
|
|
||||||
E = \frac{\int P(\mathbf{r})
|
|
||||||
\frac{\left|\Psi(\mathbf{r})\right|^2}{P(\mathbf{r})}\,
|
|
||||||
\frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int P(\mathbf{r}) \frac{\left|\Psi(\mathbf{r}) \right|^2}{P(\mathbf{r})} d\mathbf{r}}
|
|
||||||
\]
|
|
||||||
|
|
||||||
with the Gaussian probability
|
|
||||||
|
|
||||||
\[
|
|
||||||
P(\mathbf{r}) = \frac{1}{(2 \pi)^{3/2}}\exp\left( -\frac{\mathbf{r}^2}{2} \right).
|
|
||||||
\]
|
|
||||||
|
|
||||||
As the coordinates are drawn with probability $P(\mathbf{r})$, the
|
|
||||||
average energy can be computed as
|
|
||||||
|
|
||||||
$$
|
|
||||||
E \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\;
|
|
||||||
w_i = \frac{\left|\Psi(\mathbf{r}_i)\right|^2}{P(\mathbf{r}_i)} \delta \mathbf{r}
|
|
||||||
$$
|
|
||||||
|
|
||||||
|
|
||||||
*** Exercise
|
|
||||||
|
|
||||||
#+begin_exercise
|
|
||||||
Modify the program of the previous section to sample with
|
|
||||||
Gaussian-distributed random numbers. Can you see an reduction in
|
|
||||||
the statistical error?
|
|
||||||
#+end_exercise
|
|
||||||
|
|
||||||
**** Python
|
|
||||||
#+BEGIN_SRC python :results output
|
|
||||||
#!/usr/bin/env python3
|
|
||||||
|
|
||||||
from hydrogen import *
|
|
||||||
from qmc_stats import *
|
|
||||||
|
|
||||||
norm_gauss = 1./(2.*np.pi)**(1.5)
|
|
||||||
def gaussian(r):
|
|
||||||
return norm_gauss * np.exp(-np.dot(r,r)*0.5)
|
|
||||||
|
|
||||||
def MonteCarlo(a,nmax):
|
|
||||||
E = 0.
|
|
||||||
N = 0.
|
|
||||||
for istep in range(nmax):
|
|
||||||
r = np.random.normal(loc=0., scale=1.0, size=(3))
|
|
||||||
w = psi(a,r)
|
|
||||||
w = w*w / gaussian(r)
|
|
||||||
N += w
|
|
||||||
E += w * e_loc(a,r)
|
|
||||||
return E/N
|
|
||||||
|
|
||||||
a = 1.2
|
|
||||||
nmax = 100000
|
|
||||||
X = [MonteCarlo(a,nmax) for i in range(30)]
|
|
||||||
E, deltaE = ave_error(X)
|
|
||||||
print(f"E = {E} +/- {deltaE}")
|
|
||||||
#+END_SRC
|
|
||||||
|
|
||||||
#+RESULTS:
|
|
||||||
: E = -0.49511014287471955 +/- 0.00012402022172236656
|
|
||||||
|
|
||||||
**** Fortran
|
|
||||||
#+BEGIN_SRC f90
|
|
||||||
double precision function gaussian(r)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(in) :: r(3)
|
|
||||||
double precision, parameter :: norm_gauss = 1.d0/(2.d0*dacos(-1.d0))**(1.5d0)
|
|
||||||
gaussian = norm_gauss * dexp( -0.5d0 * (r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
|
|
||||||
end function gaussian
|
|
||||||
|
|
||||||
|
|
||||||
subroutine gaussian_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 :: norm, r(3), w
|
|
||||||
|
|
||||||
double precision, external :: e_loc, psi, gaussian
|
|
||||||
|
|
||||||
energy = 0.d0
|
|
||||||
norm = 0.d0
|
|
||||||
do istep = 1,nmax
|
|
||||||
call random_gauss(r,3)
|
|
||||||
w = psi(a,r)
|
|
||||||
w = w*w / gaussian(r)
|
|
||||||
norm = norm + w
|
|
||||||
energy = energy + w * e_loc(a,r)
|
|
||||||
end do
|
|
||||||
energy = energy / norm
|
|
||||||
end subroutine gaussian_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 gaussian_montecarlo(a,nmax,X(irun))
|
|
||||||
enddo
|
|
||||||
call ave_error(X,nruns,ave,err)
|
|
||||||
print *, 'E = ', ave, '+/-', err
|
|
||||||
end program qmc
|
|
||||||
#+END_SRC
|
|
||||||
|
|
||||||
#+begin_src sh :results output :exports both
|
|
||||||
gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
|
|
||||||
./qmc_gaussian
|
|
||||||
#+end_src
|
|
||||||
|
|
||||||
#+RESULTS:
|
|
||||||
: E = -0.49517104619091717 +/- 1.0685523607878961E-004
|
|
||||||
|
|
||||||
** Improved sampling with $\Psi^2$ :noexport:
|
|
||||||
|
|
||||||
*** Importance sampling
|
|
||||||
:PROPERTIES:
|
|
||||||
:header-args:python: :tangle vmc.py
|
|
||||||
:header-args:f90: :tangle vmc.f90
|
|
||||||
:END:
|
|
||||||
|
|
||||||
To generate the probability density $\Psi^2$, we consider a
|
|
||||||
diffusion process characterized by a time-dependent density
|
|
||||||
$|\Psi(\mathbf{r},t)|^2$, which obeys the Fokker-Planck equation:
|
|
||||||
|
|
||||||
\[
|
|
||||||
\frac{\partial \Psi^2}{\partial t} = \sum_i D
|
|
||||||
\frac{\partial}{\partial \mathbf{r}_i} \left(
|
|
||||||
\frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right) |\Psi(\mathbf{r},t)|^2.
|
|
||||||
\]
|
|
||||||
|
|
||||||
$D$ is the diffusion constant and $F_i$ is the i-th component of a
|
|
||||||
drift velocity caused by an external potential. For a stationary
|
|
||||||
density, \( \frac{\partial \Psi^2}{\partial t} = 0 \), so
|
|
||||||
|
|
||||||
\begin{eqnarray*}
|
|
||||||
0 & = & \sum_i D
|
|
||||||
\frac{\partial}{\partial \mathbf{r}_i} \left(
|
|
||||||
\frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right) |\Psi(\mathbf{r})|^2 \\
|
|
||||||
0 & = & \sum_i D
|
|
||||||
\frac{\partial}{\partial \mathbf{r}_i} \left(
|
|
||||||
\frac{\partial |\Psi(\mathbf{r})|^2}{\partial \mathbf{r}_i} -
|
|
||||||
F_i(\mathbf{r})\,|\Psi(\mathbf{r})|^2 \right) \\
|
|
||||||
0 & = &
|
|
||||||
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} -
|
|
||||||
\frac{\partial F_i }{\partial \mathbf{r}_i}|\Psi(\mathbf{r})|^2 -
|
|
||||||
\frac{\partial \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
|
|
||||||
\end{eqnarray*}
|
|
||||||
|
|
||||||
we search for a drift function which satisfies
|
|
||||||
|
|
||||||
\[
|
|
||||||
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} = |\Psi(\mathbf{r})|^2 \frac{\partial F_i }{\partial \mathbf{r}_i} +
|
|
||||||
\frac{\partial \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
|
|
||||||
\]
|
|
||||||
|
|
||||||
to obtain a second derivative on the left, we need the drift to be
|
|
||||||
of the form
|
|
||||||
\[
|
|
||||||
F_i(\mathbf{r}) = g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
|
|
||||||
\]
|
|
||||||
|
|
||||||
\[
|
|
||||||
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2}
|
|
||||||
= |\Psi(\mathbf{r})|^2 \frac{\partial g(\mathbf{r})}{\partial
|
|
||||||
\mathbf{r}_i}\frac{\partial \Psi^2}{\partial
|
|
||||||
\mathbf{r}_i} + |\Psi(\mathbf{r})|^2 g(\mathbf{r})
|
|
||||||
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} + \frac{\partial
|
|
||||||
\Psi^2}{\partial \mathbf{r}_i} g(\mathbf{r}) \frac{\partial
|
|
||||||
\Psi^2}{\partial \mathbf{r}_i}
|
|
||||||
\]
|
|
||||||
|
|
||||||
$g = 1 / \Psi^2$ satisfies this equation, so
|
|
||||||
|
|
||||||
\[
|
|
||||||
F(\mathbf{r}) = \frac{\nabla |\Psi(\mathbf{r})|^2}{|\Psi(\mathbf{r})|^2} = 2 \frac{\nabla
|
|
||||||
\Psi(\mathbf{r})}{\Psi(\mathbf{r})} = 2 \nabla \left( \log \Psi(\mathbf{r}) \right)
|
|
||||||
\]
|
|
||||||
|
|
||||||
In statistical mechanics, Fokker-Planck trajectories are generated
|
|
||||||
by a Langevin equation:
|
|
||||||
|
|
||||||
\[
|
|
||||||
\frac{\partial \mathbf{r}(t)}{\partial t} = 2D \frac{\nabla
|
|
||||||
\Psi(\mathbf{r}(t))}{\Psi} + \eta
|
|
||||||
\]
|
|
||||||
|
|
||||||
where $\eta$ is a normally-distributed fluctuating random force.
|
|
||||||
|
|
||||||
Discretizing this differential equation gives the following drifted
|
|
||||||
diffusion scheme:
|
|
||||||
|
|
||||||
\[
|
|
||||||
\mathbf{r}_{n+1} = \mathbf{r}_{n} + \delta t\, 2D \frac{\nabla
|
|
||||||
\Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \chi
|
|
||||||
\]
|
|
||||||
where $\chi$ is a Gaussian random variable with zero mean and
|
|
||||||
variance $\delta t\,2D$.
|
|
||||||
|
|
||||||
**** Exercise 2
|
|
||||||
|
|
||||||
#+begin_exercise
|
|
||||||
Sample $\Psi^2$ approximately using the drifted diffusion scheme,
|
|
||||||
with a diffusion constant $D=1/2$. You can use a time step of
|
|
||||||
0.001 a.u.
|
|
||||||
#+end_exercise
|
|
||||||
|
|
||||||
*Python*
|
|
||||||
#+BEGIN_SRC python :results output
|
|
||||||
#!/usr/bin/env python3
|
|
||||||
|
|
||||||
from hydrogen import *
|
|
||||||
from qmc_stats import *
|
|
||||||
|
|
||||||
def MonteCarlo(a,dt,nmax):
|
|
||||||
sq_dt = np.sqrt(dt)
|
|
||||||
|
|
||||||
# Initialization
|
|
||||||
E = 0.
|
|
||||||
N = 0.
|
|
||||||
r_old = np.random.normal(loc=0., scale=1.0, size=(3))
|
|
||||||
|
|
||||||
for istep in range(nmax):
|
|
||||||
d_old = drift(a,r_old)
|
|
||||||
chi = np.random.normal(loc=0., scale=1.0, size=(3))
|
|
||||||
r_new = r_old + dt * d_old + chi*sq_dt
|
|
||||||
N += 1.
|
|
||||||
E += e_loc(a,r_new)
|
|
||||||
r_old = r_new
|
|
||||||
return E/N
|
|
||||||
|
|
||||||
|
|
||||||
a = 1.2
|
|
||||||
nmax = 100000
|
|
||||||
dt = 0.2
|
|
||||||
X = [MonteCarlo(a,dt,nmax) for i in range(30)]
|
|
||||||
E, deltaE = ave_error(X)
|
|
||||||
print(f"E = {E} +/- {deltaE}")
|
|
||||||
#+END_SRC
|
|
||||||
|
|
||||||
#+RESULTS:
|
|
||||||
: E = -0.4858534479298907 +/- 0.00010203236131158794
|
|
||||||
|
|
||||||
*Fortran*
|
|
||||||
#+BEGIN_SRC f90
|
|
||||||
subroutine variational_montecarlo(a,dt,nmax,energy)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(in) :: a, dt
|
|
||||||
integer*8 , intent(in) :: nmax
|
|
||||||
double precision, intent(out) :: energy
|
|
||||||
|
|
||||||
integer*8 :: istep
|
|
||||||
double precision :: norm, r_old(3), r_new(3), d_old(3), sq_dt, chi(3)
|
|
||||||
double precision, external :: e_loc
|
|
||||||
|
|
||||||
sq_dt = dsqrt(dt)
|
|
||||||
|
|
||||||
! Initialization
|
|
||||||
energy = 0.d0
|
|
||||||
norm = 0.d0
|
|
||||||
call random_gauss(r_old,3)
|
|
||||||
|
|
||||||
do istep = 1,nmax
|
|
||||||
call drift(a,r_old,d_old)
|
|
||||||
call random_gauss(chi,3)
|
|
||||||
r_new(:) = r_old(:) + dt * d_old(:) + chi(:)*sq_dt
|
|
||||||
norm = norm + 1.d0
|
|
||||||
energy = energy + e_loc(a,r_new)
|
|
||||||
r_old(:) = r_new(:)
|
|
||||||
end do
|
|
||||||
energy = energy / norm
|
|
||||||
end subroutine variational_montecarlo
|
|
||||||
|
|
||||||
program qmc
|
|
||||||
implicit none
|
|
||||||
double precision, parameter :: a = 1.2d0
|
|
||||||
double precision, parameter :: dt = 0.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 variational_montecarlo(a,dt,nmax,X(irun))
|
|
||||||
enddo
|
|
||||||
call ave_error(X,nruns,ave,err)
|
|
||||||
print *, 'E = ', ave, '+/-', err
|
|
||||||
end program qmc
|
|
||||||
#+END_SRC
|
|
||||||
|
|
||||||
#+begin_src sh :results output :exports both
|
|
||||||
gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
|
|
||||||
./vmc
|
|
||||||
#+end_src
|
|
||||||
|
|
||||||
#+RESULTS:
|
|
||||||
: E = -0.48584030499187431 +/- 1.0411743995438257E-004
|
|
||||||
|
|
||||||
|
|
||||||
* Project
|
* Project
|
||||||
|
|
||||||
Change your PDMC code for one of the following:
|
Change your PDMC code for one of the following:
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user