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 so-called 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 inter-particle 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})$ (Hartree-Fock, Kohn-Sham @@ -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 time-step $\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 /time-step 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 + fixed-node 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}_{n-1}) \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_{n-1} + \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 time-step 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 time-step 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 time-step and you should not worry about the extrapolation. - - The accept/reject step (steps 2-5 in the algorithm) is in principle not needed for the correctness of + - 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. - 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 - fixed-node 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)