From cb6ceeb797a8db4af4cd26f8fd5cbca228c25112 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 11 Jan 2021 18:41:36 +0100 Subject: [PATCH] Fokker-Planck --- QMC.org | 346 ++++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 250 insertions(+), 96 deletions(-) diff --git a/QMC.org b/QMC.org index eb7da99..2d44712 100644 --- a/QMC.org +++ b/QMC.org @@ -1,8 +1,13 @@ #+TITLE: Quantum Monte Carlo #+AUTHOR: Anthony Scemama, Claudia Filippi -#+SETUPFILE: https://fniessen.github.io/org-html-themes/org/theme-readtheorg.setup +# SETUPFILE: https://fniessen.github.io/org-html-themes/org/theme-readtheorg.setup +# SETUPFILE: https://fniessen.github.io/org-html-themes/org/theme-bigblow.setup #+STARTUP: latexpreview +#+HTML_HEAD: +#+HTML_HEAD: +#+HTML_HEAD: + * Introduction @@ -14,10 +19,17 @@ computes a statistical estimate of the expectation value of the energy associated with a given wave function. Finally, we introduce the diffusion Monte Carlo (DMC) method which - gives the exact energy of the H$_2$ molecule. + gives the exact energy of the $H_2$ molecule. Code examples will be given in Python and Fortran. Whatever language can be chosen. + + We consider the stationary solution of the Schrödinger equation, so + the wave functions considered here are real: for an $N$ electron + system where the electrons move in the 3-dimensional space, + $\Psi : \mathbb{R}^{3N} \rightarrow \mathbb{R}$. In addition, $\Psi$ + is defined everywhere, continuous and infinitely differentiable. + ** Python ** Fortran @@ -45,7 +57,7 @@ $$ when $a=1$, by checking that $\hat{H}\Psi(\mathbf{r}) = E\Psi(\mathbf{r})$ for - all $\mathbf{r}$: we will check that the local energy, defined as + all $\mathbf{r}$. We will check that the local energy, defined as $$ E_L(\mathbf{r}) = \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}, @@ -54,6 +66,30 @@ is constant. + + The probabilistic /expected value/ of an arbitrary function $f(x)$ + with respect to a probability density function $p(x)$ is given by + + $$ \langle f \rangle_p = \int_{-\infty}^\infty p(x)\, f(x)\,dx $$. + + Recall that a probability density function $p(x)$ is non-negative + and integrates to one: + + $$ \int_{-\infty}^\infty p(x)\,dx = 1 $$. + + + The electronic energy of a system is the expectation value of the + local energy $E(\mathbf{r})$ with respect to the $3N$-dimensional + electron density given by the square of the wave function: + + \begin{eqnarray} + E & = & \frac{\langle \Psi| \hat{H} | \Psi\rangle}{\langle \Psi |\Psi \rangle} \\ + & = & \frac{\int \Psi(\mathbf{r})\, \hat{H} \Psi(\mathbf{r})\, d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} \\ + & = & \frac{\int \left[\Psi(\mathbf{r})\right]^2\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} \\ + & = & \frac{\int \left[\Psi(\mathbf{r})\right]^2\, E_L(\mathbf{r})\,d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} + = \langle E_L \rangle_{\Psi^2} + \end{eqnarray} + ** Local energy :PROPERTIES: :header-args:python: :tangle hydrogen.py @@ -63,9 +99,9 @@ The function accepts a 3-dimensional vector =r= as input arguments and returns the potential. - $\mathbf{r}=\sqrt{x^2 + y^2 + z^2})$, so + $\mathbf{r}=\sqrt{x^2 + y^2 + z^2}$, so $$ - V(x,y,z) = -\frac{1}{\sqrt{x^2 + y^2 + z^2})$ + V(x,y,z) = -\frac{1}{\sqrt{x^2 + y^2 + z^2}} $$ #+BEGIN_SRC python :results none @@ -175,7 +211,7 @@ double precision function e_loc(a,r) end function e_loc #+END_SRC -** Plot the local energy along the x axis +** Plot of the local energy along the $x$ axis :PROPERTIES: :header-args:python: :tangle plot_hydrogen.py :header-args:f90: :tangle plot_hydrogen.f90 @@ -184,7 +220,7 @@ end function e_loc For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the local energy along the $x$ axis. - #+begin_src python :results output + #+BEGIN_SRC python :results none import numpy as np import matplotlib.pyplot as plt @@ -270,72 +306,66 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \ #+RESULTS: [[file:plot.png]] -** Compute numerically the average energy +** Compute numerically the expectation value of the energy :PROPERTIES: :header-args:python: :tangle energy_hydrogen.py :header-args:f90: :tangle energy_hydrogen.f90 :END: - We want to compute - - \begin{eqnarray} - E & = & \frac{\langle \Psi| \hat{H} | \Psi\rangle}{\langle \Psi |\Psi \rangle} \\ - & = & \frac{\int \Psi(\mathbf{r})\, \hat{H} \Psi(\mathbf{r})\, d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} \\ - & = & \frac{\int \left[\Psi(\mathbf{r})\right]^2\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} - \end{eqnarray} - If the space is discretized in small volume elements $\delta - \mathbf{r}$, this last equation corresponds to a weighted average of - the local energy, where the weights are the values of the square of - the wave function at $\mathbf{r}$ multiplied by the volume element: + \mathbf{r}$, the expression of \langle E_L \rangle_{\Psi^2}$ becomes + a weighted average of the local energy, where the weights are the + values of the probability density at $\mathbf{r}$ multiplied + by the volume element: $$ - E \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\; + \langle E \rangle_{\Psi^2} \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\; w_i = \left[\Psi(\mathbf{r}_i)\right]^2 \delta \mathbf{r} $$ - We now compute an numerical estimate of the energy in a grid of - $50\times50\times50$ points in the range $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$. + In this section, we will compute a numerical estimate of the + energy in a grid of $50\times50\times50$ points in the range + $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$. Note: the energy is biased because: - - The energy is evaluated only inside the box - - The volume elements are not infinitely small + - The volume elements are not infinitely small (discretization error) + - The energy is evaluated only inside the box (incompleteness of the space) - #+BEGIN_SRC python :results output :exports both + #+BEGIN_SRC python :results none import numpy as np from hydrogen import e_loc, psi - interval = np.linspace(-5,5,num=50) - delta = (interval[1]-interval[0])**3 +interval = np.linspace(-5,5,num=50) +delta = (interval[1]-interval[0])**3 - r = np.array([0.,0.,0.]) +r = np.array([0.,0.,0.]) - for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: - E = 0. - norm = 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) - norm += w - E = E / norm - print(f"a = {a} \t E = {E}") +for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: + E = 0. + norm = 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) + norm += w + E = E / norm + print(f"a = {a} \t E = {E}") - #+end_src + #+end_src - #+RESULTS: - : 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 + #+RESULTS: + : 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 #+begin_src f90 @@ -367,7 +397,6 @@ program energy_hydrogen r(3) = x(l) w = psi(a(j),r) w = w * w * delta - energy = energy + w * e_loc(a(j), r) norm = norm + w end do @@ -380,7 +409,7 @@ program energy_hydrogen end program energy_hydrogen #+end_src - To compile and run: + To compile the Fortran and run it: #+begin_src sh :results output :exports both gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen @@ -401,20 +430,23 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen :header-args:f90: :tangle variance_hydrogen.f90 :END: - The variance of the local energy measures the magnitude of the - fluctuations of the local energy around the average. If the local - energy is constant (i.e. $\Psi$ is an eigenfunction of $\hat{H}$) - the variance is zero. + The variance of the local energy is a functional of $\Psi$ + which measures the magnitude of the fluctuations of the local + energy associated with $\Psi$ around the average: $$ \sigma^2(E_L) = \frac{\int \left[\Psi(\mathbf{r})\right]^2\, \left[ E_L(\mathbf{r}) - E \right]^2 \, d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} $$ + If the local energy is constant (i.e. $\Psi$ is an eigenfunction of + $\hat{H}$) the variance is zero, so the variance of the local + energy can be used as a measure of the quality of a wave function. + Compute a numerical estimate of the variance of the local energy in a grid of $50\times50\times50$ points in the range $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$. - #+BEGIN_SRC python :results output :exports both + #+begin_src python :results none import numpy as np from hydrogen import e_loc, psi @@ -451,7 +483,6 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: s2 += w * (El - E)**2 s2 = s2 / norm print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}") - #+end_src #+RESULTS: @@ -541,6 +572,8 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen * Variational Monte Carlo + Numerical integration with deterministic methods is very efficient + in low dimensions. When the number of dimensions becomes larger than Instead of computing the average energy as a numerical integration on a grid, we will do a Monte Carlo sampling, which is an extremely efficient method to compute integrals when the number of dimensions is @@ -582,7 +615,7 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen Write a function returning the average and statistical error of an input array. - #+BEGIN_SRC python + #+BEGIN_SRC python :results none from math import sqrt def ave_error(arr): M = len(arr) @@ -638,23 +671,23 @@ end subroutine ave_error from hydrogen import * from qmc_stats import * - def MonteCarlo(a, nmax): - E = 0. - N = 0. - for istep in range(nmax): - r = np.random.uniform(-5., 5., (3)) - w = psi(a,r) - w = w*w - N += w - E += w * e_loc(a,r) - return E/N +def MonteCarlo(a, nmax): + E = 0. + N = 0. + for istep in range(nmax): + r = np.random.uniform(-5., 5., (3)) + w = psi(a,r) + w = w*w + N += w + E += w * e_loc(a,r) + return E/N - a = 0.9 - nmax = 100000 - X = [MonteCarlo(a,nmax) for i in range(30)] - E, deltaE = ave_error(X) - print(f"E = {E} +/- {deltaE}") - #+END_SRC +a = 0.9 +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.4956255109300764 +/- 0.0007082875482711226 @@ -868,7 +901,12 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian #+RESULTS: : E = -0.49606057056767766 +/- 1.3918807547836872E-004 + ** Sampling with $\Psi^2$ + :PROPERTIES: + :header-args:python: :tangle vmc.py + :header-args:f90: :tangle vmc.f90 + :END: We will now use the square of the wave function to make the sampling: @@ -876,19 +914,79 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian P(\mathbf{r}) = \left[\Psi(\mathbf{r})\right]^2 \] - Now, the expression for the energy will be simplified to the - average of the local energies, each with a weight of 1. + The expression for the energy will be simplified to the average of + the local energies, each with a weight of 1. $$ - E \approx \frac{1}{M}\sum_{i=1}^M E_L(\mathbf{r}_i)} + E \approx \frac{1}{M}\sum_{i=1}^M E_L(\mathbf{r}_i) $$ - To generate the probability density $\Psi^2$, we can use a drifted - diffusion scheme: + 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) + \] + + following drifted diffusion scheme: \[ \mathbf{r}_{n+1} = \mathbf{r}_{n} + \tau \frac{\nabla - \Psi(r)}{\Psi(r)} + \eta \sqrt{\tau} + \Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \eta \sqrt{\tau} \] where $\eta$ is a normally-distributed Gaussian random number. @@ -902,7 +1000,66 @@ def drift(a,r): return r * ar_inv #+END_SRC - #+RESULTS: + #+BEGIN_SRC f90 +subroutine 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 drift + #+END_SRC + + + Now we can write the Monte Carlo sampling + #+BEGIN_SRC f90 +subroutine variational_montecarlo(a,nmax,energy) + implicit none + double precision, intent(in) :: a + integer , 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 variational_montecarlo + +program qmc + implicit none + double precision, parameter :: a = 0.9 + integer , 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 vmc.f90 -o vmc +./vmc + #+end_src #+BEGIN_SRC python def MonteCarlo(a,tau,nmax): @@ -932,11 +1089,8 @@ def MonteCarlo(a,tau,nmax): N += 1. E += e_loc(a,r_old) return E/N - #+END_SRC - #+RESULTS: - #+BEGIN_SRC python :results output nmax = 100000 tau = 0.1 X = [MonteCarlo(a,tau,nmax) for i in range(30)] @@ -950,15 +1104,15 @@ print(f"E = {E} +/- {deltaE}") * Diffusion Monte Carlo -We will now consider the H_2 molecule in a minimal basis composed of the -$1s$ orbitals of the hydrogen atoms: + 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. + $$ + \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.