From 707338aa8af92563b60b74e8532f330575a281bc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Feb 2021 00:05:04 +0100 Subject: [PATCH] Cleaning --- QMC.org | 343 -------------------------------------------------------- 1 file changed, 343 deletions(-) diff --git a/QMC.org b/QMC.org index f2ffb8a..cb23a00 100644 --- a/QMC.org +++ b/QMC.org @@ -2758,350 +2758,7 @@ gfortran hydrogen.f90 qmc_stats.f90 pdmc.f90 -o pdmc : 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 Change your PDMC code for one of the following: