From c10c7697b9db9dd276d794592ca707768659d3ba Mon Sep 17 00:00:00 2001
From: Anthony Scemama
Date: Tue, 2 Feb 2021 23:41:07 +0100
Subject: [PATCH] Update PDMC

QMC.org  101 ++++++++++++++++++++++++++++++++++
1 file changed, 62 insertions(+), 39 deletions()
diff git a/QMC.org b/QMC.org
index c0d6901..f7e4eee 100644
 a/QMC.org
+++ b/QMC.org
@@ 2235,10 +2235,11 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
where $\chi$ is a Gaussian random variable, and the potential term
can be simulated by creating or destroying particles over time (a
socalled branching process) or by simply considering it as a
 cumulative multiplicative weight along the diffusion trajectory:
+ cumulative multiplicative weight along the diffusion trajectory
+ (pure Diffusion Monte Carlo):
\[
 \exp \left( \int_0^\tau  (E_L(\mathbf{r}_t)  E_{\text{ref}}) dt \right)
+ \exp \left( \int_0^\tau  (V(\mathbf{r}_t)  E_{\text{ref}}) dt \right).
\]
@@ 2260,7 +2261,7 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
In a molecular system, the potential is far from being constant
and, in fact, diverges at the interparticle coalescence points. Hence,
 it results in very large fluctuations of the term associated with
+ it results in very large fluctuations of the erm weight associated with
the potental, making the calculations impossible in practice.
Fortunately, if we multiply the SchrÃ¶dinger equation by a chosen
/trial wave function/ $\Psi_T(\mathbf{r})$ (HartreeFock, KohnSham
@@ 2284,17 +2285,19 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
scheme presented in the previous section (VMC).
The new "potential" is the local energy, which has smaller fluctuations
when $\Psi_T$ gets closer to the exact wave function.
 This term can be simulated by t particles according to $\exp\left[ \delta t\,
 \left(E_L(\mathbf{r})  E_{\rm ref}\right)\right]$
+ This term can be simulated by
+ \[
+ \exp \left( \int_0^\tau  (E_L(\mathbf{r}_t)  E_{\text{ref}}) dt \right).
+ \]
where $E_{\rm ref}$ is the constant we had introduced above, which is adjusted to
 the running average energy to keep the number of particles
 reasonably constant.
+ an estimate of the average energy to keep the weights close to one.
This equation generates the /N/electron density $\Pi$, which is the
 product of the ground state with the trial wave function. You may then ask: how
 can we compute the total energy of the system?
+ product of the ground state solution with the trial wave
+ function. You may then ask: how can we compute the total energy of
+ the system?
 To this aim, we use the mixed estimator of the energy:
+ To this aim, we use the /mixed estimator/ of the energy:
\begin{eqnarray*}
E(\tau) &=& \frac{\langle \psi(\tau)  \hat{H}  \Psi_T \rangle}{\langle \psi(\tau)  \Psi_T \rangle}\\
@@ 2310,7 +2313,8 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
\Pi(\mathbf{r},\tau) =\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \rightarrow \Phi_0(\mathbf{r}) \Psi_T(\mathbf{r})\,,
\]
 and, using that $\hat{H}$ is Hermitian and that $\Phi_0$ is an eigenstate of the Hamiltonian, we obtain for large $\tau$
+ and, using that $\hat{H}$ is Hermitian and that $\Phi_0$ is an
+ eigenstate of the Hamiltonian, we obtain for large $\tau$
\[
E(\tau) = \frac{\langle \psi_\tau  \hat{H}  \Psi_T \rangle}
@@ 2372,12 +2376,12 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
\right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau)
\]
** Pure Diffusion Monte Carlo (PDMC)
+** Pure Diffusion Monte Carlo
Instead of having a variable number of particles to simulate the
 branching process, one can consider the term
 $\exp \left( \delta t\,( E_L(\mathbf{r})  E_{\rm ref}) \right)$ as a
 cumulative product of weights:
+ 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, \tau)
@@ 2387,21 +2391,39 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
\prod_{i=1}^{n} w(\mathbf{r}_i)
\end{eqnarray*}
 where $\mathbf{r}_i$ are the coordinates along the trajectory and we introduced a timestep $\delta t$.

 The algorithm can be rather easily built on top of your VMC code:
+ where $\mathbf{r}_i$ are the coordinates along the trajectory and
+ we introduced a /timestep variable/ $\delta t$ to discretize the
+ integral.
 0) Start with $W=1$
 1) Evaluate the local energy at $\mathbf{r}_{n}$ and accumulate it
 2) Compute the weight $w(\mathbf{r}_n)$
 3) Update $W$
 4) Compute a new position $\mathbf{r'} = \mathbf{r}_n +
+ The PDMC algorithm is less stable than the DMC algorithm: it
+ requires to have a value of $E_\text{ref}$ which is close to the
+ fixednode energy, and a good trial wave function. Moreover, we
+ can't let $\tau$ become too large because the weight whether
+ explode or vanish: we need to have a fixed value of $\tau$
+ (projection time).
+ The big advantage of PDMC is that it is rather simple to implement
+ starting from a VMC code:
+
+ 0) Start with $W(\mathbf{r}_0)=1, \tau_0 = 0$
+ 1) Evaluate the local energy at $\mathbf{r}_{n}$
+ 2) Compute the contribution to the weight $w(\mathbf{r}_n) =
+ \exp(\delta t(E_L(\mathbf{r}_n)E_\text{ref}))$
+ 3) Update $W(\mathbf{r}_{n}) = W(\mathbf{r}_{n1}) \times w(\mathbf{r}_n)$
+ 4) Accumulate the weighted energy $W(\mathbf{r}_n) \times
+ E_L(\mathbf{r}_n)$,
+ and the weight $W(\mathbf{r}_n)$ for the normalization
+ 5) Update $\tau_n = \tau_{n1} + \delta t$
+ 6) If $\tau_{n} > \tau_\text{max}$, the long projection time has
+ been reached and we can start an new trajectory from the current
+ position: reset $W(r_n) = 1$ and $\tau_n
+ = 0$
+ 7) Compute a new position $\mathbf{r'} = \mathbf{r}_n +
\delta t\, \frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \chi$
 5) Evaluate $\Psi(\mathbf{r}')$ and $\frac{\nabla \Psi(\mathbf{r'})}{\Psi(\mathbf{r'})}$ at the new position
 6) Compute the ratio $A = \frac{T(\mathbf{r}' \rightarrow \mathbf{r}_{n}) P(\mathbf{r}')}{T(\mathbf{r}_{n} \rightarrow \mathbf{r}') P(\mathbf{r}_{n})}$
 7) Draw a uniform random number $v \in [0,1]$
 8) if $v \le A$, accept the move : set $\mathbf{r}_{n+1} = \mathbf{r'}$
 9) else, reject the move : set $\mathbf{r}_{n+1} = \mathbf{r}_n$
+ 8) Evaluate $\Psi(\mathbf{r}')$ and $\frac{\nabla \Psi(\mathbf{r'})}{\Psi(\mathbf{r'})}$ at the new position
+ 9) Compute the ratio $A = \frac{T(\mathbf{r}' \rightarrow \mathbf{r}_{n}) P(\mathbf{r}')}{T(\mathbf{r}_{n} \rightarrow \mathbf{r}') P(\mathbf{r}_{n})}$
+ 10) Draw a uniform random number $v \in [0,1]$
+ 11) if $v \le A$, accept the move : set $\mathbf{r}_{n+1} = \mathbf{r'}$
+ 12) else, reject the move : set $\mathbf{r}_{n+1} = \mathbf{r}_n$
Some comments are needed:
@@ 2412,19 +2434,16 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
E = \frac{\sum_{k=1}^{N_{\rm MC}} E_L(\mathbf{r}_k) W(\mathbf{r}_k, k\delta t)}{\sum_{k=1}^{N_{\rm MC}} W(\mathbf{r}_k, k\delta t)}
\end{eqnarray*}
  The result will be affected by a timestep error (the finite size of $\delta t$) and one
 has in principle to extrapolate to the limit $\delta t \rightarrow 0$. This amounts to fitting
 the energy computed for multiple values of $\delta t$.
+  The result will be affected by a timestep error
+ (the finite size of $\delta t$) due to the discretization of the
+ integral, and one has in principle to extrapolate to the limit
+ $\delta t \rightarrow 0$. This amounts to fitting the energy
+ computed for multiple values of $\delta t$.
Here, you will be using a small enough timestep and you should not worry about the extrapolation.
  The accept/reject step (steps 25 in the algorithm) is in principle not needed for the correctness of
+  The accept/reject step (steps 912 in the algorithm) is in principle not needed for the correctness of
the DMC algorithm. However, its use reduces significantly the timestep error.
 The PDMC algorithm is less stable than the branching algorithm: it
 requires to have a value of $E_\text{ref}$ which is close to the
 fixednode energy, and a good trial wave function. Its big
 advantage is that it is very easy to program starting from a VMC
 code, so this is what we will do in the next section.
** Hydrogen atom
:PROPERTIES:
@@ 2435,9 +2454,10 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 o vmc_metropolis
*** Exercise
#+begin_exercise
 Modify the Metropolis VMC program to introduce the PDMC weight.
+ Modify the Metropolis VMC program into a PDMC program.
In the limit $\delta t \rightarrow 0$, you should recover the exact
 energy of H for any value of $a$.
+ energy of H for any value of $a$, as long as the simulation is stable.
+ We choose here a fixed projection time $\tau=10$ a.u.
#+end_exercise
*Python*
@@ 2452,6 +2472,7 @@ def MonteCarlo(a, nmax, dt, Eref):
a = 1.2
nmax = 100000
dt = 0.01
+tau = 10.
E_ref = 0.5
X0 = [ MonteCarlo(a, nmax, dt, E_ref) for i in range(30)]
@@ 2467,6 +2488,8 @@ A, deltaA = ave_error(X)
print(f"A = {A} +/ {deltaA}")
#+END_SRC
+ #+RESULTS:
+
*Fortran*
#+BEGIN_SRC f90 :tangle none
subroutine pdmc(a, dt, nmax, energy, accep, tau, E_ref)