diff --git a/QMC.org b/QMC.org index c53902e..a2e7e1e 100644 --- a/QMC.org +++ b/QMC.org @@ -9,7 +9,7 @@ #+LATEX_HEADER_EXTRA: \usepackage{minted} #+HTML_HEAD: -#+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t +#+OPTIONS: H:4 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t #+OPTIONS: TeX:t LaTeX:t skip:nil d:nil todo:t pri:nil tags:not-in-toc # EXCLUDE_TAGS: Python solution # EXCLUDE_TAGS: Fortran solution @@ -36,6 +36,9 @@ * Introduction + This web site is the QMC tutorial of the LTTC winter school + [[https://www.irsamc.ups-tlse.fr/lttc/Luchon][Tutorials in Theoretical Chemistry]]. + We propose different exercises to understand quantum Monte Carlo (QMC) methods. In the first section, we propose to compute the energy of a hydrogen atom using numerical integration. The goal of this section is @@ -58,7 +61,6 @@ All the quantities are expressed in /atomic units/ (energies, coordinates, etc). - * Numerical evaluation of the energy In this section we consider the Hydrogen atom with the following @@ -128,12 +130,6 @@ *** Exercise 1 - #+begin_exercise - Find the theoretical value of $a$ for which $\Psi$ is an eigenfunction of $\hat{H}$. - #+end_exercise - -*** Exercise 2 - #+begin_exercise Write a function which computes the potential at $\mathbf{r}$. The function accepts a 3-dimensional vector =r= as input arguments @@ -145,7 +141,7 @@ V(\mathbf{r}) = -\frac{1}{\sqrt{x^2 + y^2 + z^2}} $$ -**** Python + *Python* #+BEGIN_SRC python :results none :tangle none import numpy as np @@ -153,7 +149,17 @@ def potential(r): # TODO #+END_SRC -**** Python :solution: + *Fortran* + #+BEGIN_SRC f90 :tangle none +double precision function potential(r) + implicit none + double precision, intent(in) :: r(3) + ! TODO +end function potential + #+END_SRC + +**** Solution :solution: + *Python* #+BEGIN_SRC python :results none import numpy as np @@ -163,16 +169,7 @@ def potential(r): return -1. / distance #+END_SRC -**** Fortran - #+BEGIN_SRC f90 :tangle none -double precision function potential(r) - implicit none - double precision, intent(in) :: r(3) - ! TODO -end function potential - #+END_SRC - -**** Fortran :solution: + *Fortran* #+BEGIN_SRC f90 double precision function potential(r) implicit none @@ -187,26 +184,20 @@ double precision function potential(r) end function potential #+END_SRC -*** Exercise 3 +*** Exercise 2 #+begin_exercise Write a function which computes the wave function at $\mathbf{r}$. The function accepts a scalar =a= and a 3-dimensional vector =r= as input arguments, and returns a scalar. #+end_exercise -**** Python + *Python* #+BEGIN_SRC python :results none :tangle none def psi(a, r): # TODO #+END_SRC -**** Python :solution: - #+BEGIN_SRC python :results none -def psi(a, r): - return np.exp(-a*np.sqrt(np.dot(r,r))) - #+END_SRC - -**** Fortran + *Fortran* #+BEGIN_SRC f90 :tangle none double precision function psi(a, r) implicit none @@ -215,7 +206,14 @@ double precision function psi(a, r) end function psi #+END_SRC -**** Fortran :solution: +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results none +def psi(a, r): + return np.exp(-a*np.sqrt(np.dot(r,r))) + #+END_SRC + + *Fortran* #+BEGIN_SRC f90 double precision function psi(a, r) implicit none @@ -224,7 +222,7 @@ double precision function psi(a, r) end function psi #+END_SRC -*** Exercise 4 +*** Exercise 3 #+begin_exercise Write a function which computes the local kinetic energy at $\mathbf{r}$. The function accepts =a= and =r= as input arguments and returns the @@ -261,21 +259,13 @@ end function psi -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) $$ -**** Python + *Python* #+BEGIN_SRC python :results none :tangle none def kinetic(a,r): # TODO #+END_SRC -**** Python :solution: - #+BEGIN_SRC python :results none -def kinetic(a,r): - distance = np.sqrt(np.dot(r,r)) - assert (distance > 0.) - return -0.5 * (a**2 - (2.*a)/distance) - #+END_SRC - -**** Fortran + *Fortran* #+BEGIN_SRC f90 :tangle none double precision function kinetic(a,r) implicit none @@ -284,7 +274,16 @@ double precision function kinetic(a,r) end function kinetic #+END_SRC -**** Fortran :solution: +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results none +def kinetic(a,r): + distance = np.sqrt(np.dot(r,r)) + assert (distance > 0.) + return -0.5 * (a**2 - (2.*a)/distance) + #+END_SRC + + *Fortran* #+BEGIN_SRC f90 double precision function kinetic(a,r) implicit none @@ -299,7 +298,7 @@ double precision function kinetic(a,r) end function kinetic #+END_SRC -*** Exercise 5 +*** Exercise 4 #+begin_exercise Write a function which computes the local energy at $\mathbf{r}$, using the previously defined functions. @@ -312,19 +311,13 @@ end function kinetic $$ -**** Python + *Python* #+BEGIN_SRC python :results none :tangle none def e_loc(a,r): #TODO #+END_SRC -**** Python :solution: - #+BEGIN_SRC python :results none -def e_loc(a,r): - return kinetic(a,r) + potential(r) - #+END_SRC - -**** Fortran + *Fortran* #+BEGIN_SRC f90 :tangle none double precision function e_loc(a,r) implicit none @@ -333,7 +326,14 @@ double precision function e_loc(a,r) end function e_loc #+END_SRC -**** Fortran :solution: +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results none +def e_loc(a,r): + return kinetic(a,r) + potential(r) + #+END_SRC + + *Fortran* #+BEGIN_SRC f90 double precision function e_loc(a,r) implicit none @@ -343,6 +343,26 @@ double precision function e_loc(a,r) end function e_loc #+END_SRC +*** Exercise 5 + + #+begin_exercise + Find the theoretical value of $a$ for which $\Psi$ is an eigenfunction of $\hat{H}$. + #+end_exercise + +**** Solution :solution: + + \begin{eqnarray*} + E &=& \frac{\hat{H} \Psi}{\Psi} = - \frac{1}{2} \frac{\Delta \Psi}{\Psi} - + \frac{1}{|\mathbf{r}|} \\ + &=& -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) - + \frac{1}{|\mathbf{r}|} \\ + &=& + -\frac{1}{2} a^2 + \frac{a-1}{\mathbf{|r|}} + \end{eqnarray*} + + $a=1$ cancels the $1/|r|$ term, and makes the energy constant, + equal to -0.5 atomic units. + ** Plot of the local energy along the $x$ axis :PROPERTIES: :header-args:python: :tangle plot_hydrogen.py @@ -363,7 +383,7 @@ end function e_loc Gnuplot to plot the files. #+end_exercise -**** Python + *Python* #+BEGIN_SRC python :results none :tangle none import numpy as np import matplotlib.pyplot as plt @@ -380,30 +400,7 @@ plt.legend() plt.savefig("plot_py.png") #+end_src -**** Python :solution: - #+BEGIN_SRC python :results none -import numpy as np -import matplotlib.pyplot as plt - -from hydrogen import e_loc - -x=np.linspace(-5,5) -plt.figure(figsize=(10,5)) - -for a in [0.1, 0.2, 0.5, 1., 1.5, 2.]: - y=np.array([ e_loc(a, np.array([t,0.,0.]) ) for t in x]) - plt.plot(x,y,label=f"a={a}") - -plt.tight_layout() -plt.legend() -plt.savefig("plot_py.png") - #+end_src - - #+RESULTS: - - [[./plot_py.png]] - -**** Fortran + *Fortran* #+begin_src f90 :tangle none program plot implicit none @@ -431,7 +428,7 @@ gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen To plot the data using Gnuplot: - #+begin_src gnuplot :file plot.png :exports both + #+begin_src gnuplot :file plot.png :exports code set grid set xrange [-5:5] set yrange [-2:1] @@ -443,7 +440,31 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \ './data' index 5 using 1:2 with lines title 'a=2.0' #+end_src -**** Fortran :solution: +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results none +import numpy as np +import matplotlib.pyplot as plt + +from hydrogen import e_loc + +x=np.linspace(-5,5) +plt.figure(figsize=(10,5)) + +for a in [0.1, 0.2, 0.5, 1., 1.5, 2.]: + y=np.array([ e_loc(a, np.array([t,0.,0.]) ) for t in x]) + plt.plot(x,y,label=f"a={a}") + +plt.tight_layout() +plt.legend() +plt.savefig("plot_py.png") + #+end_src + + #+RESULTS: + + [[./plot_py.png]] + + *Fortran* #+begin_src f90 program plot implicit none @@ -475,18 +496,12 @@ program plot end program plot #+end_src - To compile and run: - - #+begin_src sh :exports both + #+begin_src sh :exports none gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen ./plot_hydrogen > data #+end_src - #+RESULTS: - - To plot the data using Gnuplot: - - #+begin_src gnuplot :file plot.png :exports both + #+begin_src gnuplot :file plot.png :exports results set grid set xrange [-5:5] set yrange [-2:1] @@ -497,7 +512,6 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \ './data' index 4 using 1:2 with lines title 'a=1.5', \ './data' index 5 using 1:2 with lines title 'a=2.0' #+end_src - #+RESULTS: [[file:plot.png]] @@ -532,7 +546,7 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \ \mathbf{r} \le (5,5,5)$. #+end_exercise -**** Python + *Python* #+BEGIN_SRC python :results none :tangle none import numpy as np from hydrogen import e_loc, psi @@ -548,8 +562,40 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: #+end_src -**** Python :solution: - #+BEGIN_SRC python :results none + *Fortran* + #+begin_src f90 +program energy_hydrogen + implicit none + double precision, external :: e_loc, psi + double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm + integer :: i, k, l, j + + a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) + + dx = 10.d0/(size(x)-1) + do i=1,size(x) + x(i) = -5.d0 + (i-1)*dx + end do + + do j=1,size(a) + ! TODO + print *, 'a = ', a(j), ' E = ', energy + end do + +end program energy_hydrogen + #+end_src + + To compile the Fortran and run it: + + #+begin_src sh :results output :exports code + +gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen +./energy_hydrogen + #+end_src + +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results none :exports both import numpy as np from hydrogen import e_loc, psi @@ -585,37 +631,7 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: : a = 1.5 E = -0.39242967082602226 : a = 2.0 E = -0.08086980667844901 -**** Fortran - #+begin_src f90 -program energy_hydrogen - implicit none - double precision, external :: e_loc, psi - double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm - integer :: i, k, l, j - - a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) - - dx = 10.d0/(size(x)-1) - do i=1,size(x) - x(i) = -5.d0 + (i-1)*dx - end do - - do j=1,size(a) - ! TODO - print *, 'a = ', a(j), ' E = ', energy - end do - -end program energy_hydrogen - #+end_src - - To compile the Fortran and run it: - - #+begin_src sh :results output :exports both -gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen -./energy_hydrogen - #+end_src - -**** Fortran :solution: + *Fortran* #+begin_src f90 program energy_hydrogen implicit none @@ -657,9 +673,7 @@ program energy_hydrogen end program energy_hydrogen #+end_src - To compile the Fortran and run it: - - #+begin_src sh :results output :exports both + #+begin_src sh :results output :exports results gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen ./energy_hydrogen #+end_src @@ -700,6 +714,7 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen $$\left( \langle E - \langle E \rangle_{\Psi^2} \rangle_{\Psi^2} \right)^2 = \langle E^2 \rangle_{\Psi^2} - \langle E \rangle_{\Psi^2}^2 $$ #+end_exercise +**** TODO Solution :solution: *** Exercise #+begin_exercise Add the calculation of the variance to the previous code, and @@ -709,7 +724,7 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen \le \mathbf{r} \le (5,5,5)$ for different values of $a$. #+end_exercise -**** Python + *Python* #+begin_src python :results none :tangle none import numpy as np from hydrogen import e_loc, psi @@ -724,8 +739,44 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}") #+end_src -**** Python :solution: - #+begin_src python :results none + *Fortran* + #+begin_src f90 :tangle none +program variance_hydrogen + implicit none + double precision, external :: e_loc, psi + double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm, s2 + double precision :: e, energy2 + integer :: i, k, l, j + + a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) + + dx = 10.d0/(size(x)-1) + do i=1,size(x) + x(i) = -5.d0 + (i-1)*dx + end do + + delta = dx**3 + + r(:) = 0.d0 + + do j=1,size(a) + ! TODO + print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2 + end do + +end program variance_hydrogen + #+end_src + + To compile and run: + + #+begin_src sh :results output :exports both +gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen +./variance_hydrogen + #+end_src + +**** Solution :solution: + *Python* + #+begin_src python :results none :exports both import numpy as np from hydrogen import e_loc, psi @@ -765,42 +816,7 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: : a = 1.5 E = -0.39242967 \sigma^2 = 0.31449671 : a = 2.0 E = -0.08086981 \sigma^2 = 1.80688143 -**** Fortran - #+begin_src f90 :tangle none -program variance_hydrogen - implicit none - double precision, external :: e_loc, psi - double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm, s2 - double precision :: e, energy2 - integer :: i, k, l, j - - a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) - - dx = 10.d0/(size(x)-1) - do i=1,size(x) - x(i) = -5.d0 + (i-1)*dx - end do - - delta = dx**3 - - r(:) = 0.d0 - - do j=1,size(a) - ! TODO - print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2 - end do - -end program variance_hydrogen - #+end_src - - To compile and run: - - #+begin_src sh :results output :exports both -gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen -./variance_hydrogen - #+end_src - -**** Fortran :solution: + *Fortran* #+begin_src f90 program variance_hydrogen implicit none @@ -848,9 +864,7 @@ program variance_hydrogen end program variance_hydrogen #+end_src - To compile and run: - - #+begin_src sh :results output :exports both + #+begin_src sh :results output :exports results gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen ./variance_hydrogen #+end_src @@ -863,7 +877,6 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen : a = 1.5000000000000000 E = -0.39242967082602065 s2 = 0.31449670909172917 : a = 2.0000000000000000 E = -8.0869806678448772E-002 s2 = 1.8068814270846534 - * Variational Monte Carlo Numerical integration with deterministic methods is very efficient @@ -910,7 +923,7 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen input array. #+end_exercise -**** Python + *Python* #+BEGIN_SRC python :results none :tangle none from math import sqrt def ave_error(arr): @@ -918,8 +931,20 @@ def ave_error(arr): return (average, error) #+END_SRC -**** Python :solution: - #+BEGIN_SRC python :results none + *Fortran* + #+BEGIN_SRC f90 +subroutine ave_error(x,n,ave,err) + implicit none + integer, intent(in) :: n + double precision, intent(in) :: x(n) + double precision, intent(out) :: ave, err + ! TODO +end subroutine ave_error + #+END_SRC + +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results none :exports code from math import sqrt def ave_error(arr): M = len(arr) @@ -933,19 +958,8 @@ def ave_error(arr): return (average, error) #+END_SRC -**** Fortran - #+BEGIN_SRC f90 -subroutine ave_error(x,n,ave,err) - implicit none - integer, intent(in) :: n - double precision, intent(in) :: x(n) - double precision, intent(out) :: ave, err - ! TODO -end subroutine ave_error - #+END_SRC - -**** Fortran :solution: - #+BEGIN_SRC f90 + *Fortran* + #+BEGIN_SRC f90 :exports both subroutine ave_error(x,n,ave,err) implicit none integer, intent(in) :: n @@ -1001,14 +1015,13 @@ end subroutine ave_error compute the average energy and the associated error bar. #+end_exercise -**** Python - + *Python* #+begin_note To draw a uniform random number in Python, you can use the [[https://numpy.org/doc/stable/reference/random/generated/numpy.random.uniform.html][~random.uniform~]] function of Numpy. #+end_note - #+BEGIN_SRC python :results output :tangle none + #+BEGIN_SRC python :tangle none :exports code from hydrogen import * from qmc_stats import * @@ -1021,36 +1034,7 @@ nmax = 100000 print(f"E = {E} +/- {deltaE}") #+END_SRC - #+RESULTS: - : E = -0.4956255109300764 +/- 0.0007082875482711226 - -**** Python :solution: - #+BEGIN_SRC python :results output -from hydrogen import * -from qmc_stats import * - -def MonteCarlo(a, nmax): - energy = 0. - normalization = 0. - for istep in range(nmax): - r = np.random.uniform(-5., 5., (3)) - w = psi(a,r) - w = w*w - normalization += w - energy += w * e_loc(a,r) - return energy/normalization - -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 - -**** Fortran + *Fortran* #+begin_note To draw a uniform random number in Fortran, you can use the [[https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fNUMBER.html][~RANDOM_NUMBER~]] subroutine. @@ -1094,20 +1078,47 @@ program qmc end program qmc #+END_SRC - #+begin_src sh :results output :exports both + #+begin_src sh :results output :exports code gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform ./qmc_uniform #+end_src -**** Fortran :solution: - #+begin_note - When running Monte Carlo calculations, the number of steps is - usually very large. We expect =nmax= to be possibly larger than 2 - billion, so we use 8-byte integers (=integer*8=) to represent it, as - well as the index of the current step. - #+end_note +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results output :exports both +from hydrogen import * +from qmc_stats import * - #+BEGIN_SRC f90 +def MonteCarlo(a, nmax): + energy = 0. + normalization = 0. + for istep in range(nmax): + r = np.random.uniform(-5., 5., (3)) + w = psi(a,r) + w = w*w + normalization += w + energy += w * e_loc(a,r) + return energy/normalization + +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 + + *Fortran* + #+begin_note + When running Monte Carlo calculations, the number of steps is + usually very large. We expect =nmax= to be possibly larger than 2 + billion, so we use 8-byte integers (=integer*8=) to represent it, as + well as the index of the current step. + #+end_note + + #+BEGIN_SRC f90 :exports code subroutine uniform_montecarlo(a,nmax,energy) implicit none double precision, intent(in) :: a @@ -1149,15 +1160,15 @@ program qmc call ave_error(X,nruns,ave,err) print *, 'E = ', ave, '+/-', err end program qmc - #+END_SRC + #+END_SRC - #+begin_src sh :results output :exports both + #+begin_src sh :results output :exports results gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform ./qmc_uniform - #+end_src + #+end_src - #+RESULTS: - : E = -0.49588321986667677 +/- 7.1758863546737969E-004 + #+RESULTS: + : E = -0.49588321986667677 +/- 7.1758863546737969E-004 ** Metropolis sampling with $\Psi^2$ :PROPERTIES: @@ -1237,7 +1248,7 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform Can you observe a reduction in the statistical error? #+end_exercise -**** Python + *Python* #+BEGIN_SRC python :results output :tangle none from hydrogen import * from qmc_stats import * @@ -1263,54 +1274,7 @@ A, deltaA = ave_error(X) print(f"A = {A} +/- {deltaA}") #+END_SRC - #+RESULTS: - : E = -0.4950720838131573 +/- 0.00019089638602238043 - : A = 0.5172960000000001 +/- 0.0003443446549306529 - -**** Python :solution: - #+BEGIN_SRC python :results output -from hydrogen import * -from qmc_stats import * - -def MonteCarlo(a,nmax,tau): - energy = 0. - N_accep = 0 - r_old = np.random.uniform(-tau, tau, (3)) - psi_old = psi(a,r_old) - for istep in range(nmax): - r_new = r_old + np.random.uniform(-tau,tau,(3)) - psi_new = psi(a,r_new) - ratio = (psi_new / psi_old)**2 - v = np.random.uniform() - if v <= ratio: - N_accep += 1 - r_old = r_new - psi_old = psi_new - energy += e_loc(a,r_old) - return energy/nmax, N_accep/nmax - -# Run simulation -a = 0.9 -nmax = 100000 -tau = 1.3 -X0 = [ MonteCarlo(a,nmax,tau) for i in range(30)] - -# Energy -X = [ x for (x, _) in X0 ] -E, deltaE = ave_error(X) -print(f"E = {E} +/- {deltaE}") - -# Acceptance rate -X = [ x for (_, x) in X0 ] -A, deltaA = ave_error(X) -print(f"A = {A} +/- {deltaA}") - #+END_SRC - - #+RESULTS: - : E = -0.4950720838131573 +/- 0.00019089638602238043 - : A = 0.5172960000000001 +/- 0.0003443446549306529 - -**** Fortran + *Fortran* #+BEGIN_SRC f90 :tangle none subroutine metropolis_montecarlo(a,nmax,tau,energy,accep) implicit none @@ -1358,8 +1322,52 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_metropolis.f90 -o qmc_metropolis ./qmc_metropolis #+end_src -**** Fortran :solution: - #+BEGIN_SRC f90 +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results output :exports both +from hydrogen import * +from qmc_stats import * + +def MonteCarlo(a,nmax,tau): + energy = 0. + N_accep = 0 + r_old = np.random.uniform(-tau, tau, (3)) + psi_old = psi(a,r_old) + for istep in range(nmax): + r_new = r_old + np.random.uniform(-tau,tau,(3)) + psi_new = psi(a,r_new) + ratio = (psi_new / psi_old)**2 + v = np.random.uniform() + if v <= ratio: + N_accep += 1 + r_old = r_new + psi_old = psi_new + energy += e_loc(a,r_old) + return energy/nmax, N_accep/nmax + +# Run simulation +a = 0.9 +nmax = 100000 +tau = 1.3 +X0 = [ MonteCarlo(a,nmax,tau) for i in range(30)] + +# Energy +X = [ x for (x, _) in X0 ] +E, deltaE = ave_error(X) +print(f"E = {E} +/- {deltaE}") + +# Acceptance rate +X = [ x for (_, x) in X0 ] +A, deltaA = ave_error(X) +print(f"A = {A} +/- {deltaA}") + #+END_SRC + + #+RESULTS: + : E = -0.4950720838131573 +/- 0.00019089638602238043 + : A = 0.5172960000000001 +/- 0.0003443446549306529 + + *Fortran* + #+BEGIN_SRC f90 :exports code subroutine metropolis_montecarlo(a,nmax,tau,energy,accep) implicit none double precision, intent(in) :: a @@ -1420,7 +1428,7 @@ program qmc end program qmc #+END_SRC - #+begin_src sh :results output :exports both + #+begin_src sh :results output :exports results gfortran hydrogen.f90 qmc_stats.f90 qmc_metropolis.f90 -o qmc_metropolis ./qmc_metropolis #+end_src @@ -1554,20 +1562,13 @@ end subroutine random_gauss Write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$. #+end_exercise -**** Python + *Python* #+BEGIN_SRC python :tangle hydrogen.py :tangle none def drift(a,r): # TODO #+END_SRC -**** Python :solution: - #+BEGIN_SRC python :tangle hydrogen.py -def drift(a,r): - ar_inv = -a/np.sqrt(np.dot(r,r)) - return r * ar_inv - #+END_SRC - -**** Fortran + *Fortran* #+BEGIN_SRC f90 :tangle hydrogen.f90 :tangle none subroutine drift(a,r,b) implicit none @@ -1577,7 +1578,15 @@ subroutine drift(a,r,b) end subroutine drift #+END_SRC -**** Fortran :solution: +**** Solution :solution: + *Python* + #+BEGIN_SRC python :tangle hydrogen.py +def drift(a,r): + ar_inv = -a/np.sqrt(np.dot(r,r)) + return r * ar_inv + #+END_SRC + + *Fortran* #+BEGIN_SRC f90 :tangle hydrogen.f90 subroutine drift(a,r,b) implicit none @@ -1596,7 +1605,7 @@ end subroutine drift (This is a necessary step for the next section). #+end_exercise -**** Python + *Python* #+BEGIN_SRC python :results output :tangle none from hydrogen import * from qmc_stats import * @@ -1621,8 +1630,54 @@ A, deltaA = ave_error(X) print(f"A = {A} +/- {deltaA}") #+END_SRC -**** Python :solution: - #+BEGIN_SRC python :results output + *Fortran* + #+BEGIN_SRC f90 :tangle none +subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate) + implicit none + double precision, intent(in) :: a, tau + integer*8 , intent(in) :: nmax + double precision, intent(out) :: energy, accep_rate + + integer*8 :: istep + double precision :: sq_tau, chi(3) + double precision :: psi_old, psi_new + double precision :: r_old(3), r_new(3) + double precision :: d_old(3), d_new(3) + double precision, external :: e_loc, psi + + ! TODO + +end subroutine variational_montecarlo + +program qmc + implicit none + double precision, parameter :: a = 0.9 + double precision, parameter :: tau = 1.0 + integer*8 , parameter :: nmax = 100000 + integer , parameter :: nruns = 30 + + integer :: irun + double precision :: X(nruns), accep(nruns) + double precision :: ave, err + + do irun=1,nruns + call variational_montecarlo(a,tau,nmax,X(irun),accep(irun)) + enddo + call ave_error(X,nruns,ave,err) + print *, 'E = ', ave, '+/-', err + call ave_error(accep,nruns,ave,err) + print *, 'A = ', ave, '+/-', err +end program qmc + #+END_SRC + + #+begin_src sh :results output :exports code +gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis +./vmc_metropolis + #+end_src + +**** Solution :solution: + *Python* + #+BEGIN_SRC python :results output :exports both from hydrogen import * from qmc_stats import * @@ -1676,52 +1731,7 @@ print(f"A = {A} +/- {deltaA}") : E = -0.4951317910667116 +/- 0.00014045774335059988 : A = 0.7200673333333333 +/- 0.00045942791345632793 -**** Fortran - #+BEGIN_SRC f90 :tangle none -subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate) - implicit none - double precision, intent(in) :: a, tau - integer*8 , intent(in) :: nmax - double precision, intent(out) :: energy, accep_rate - - integer*8 :: istep - double precision :: sq_tau, chi(3) - double precision :: psi_old, psi_new - double precision :: r_old(3), r_new(3) - double precision :: d_old(3), d_new(3) - double precision, external :: e_loc, psi - - ! TODO - -end subroutine variational_montecarlo - -program qmc - implicit none - double precision, parameter :: a = 0.9 - double precision, parameter :: tau = 1.0 - integer*8 , parameter :: nmax = 100000 - integer , parameter :: nruns = 30 - - integer :: irun - double precision :: X(nruns), accep(nruns) - double precision :: ave, err - - do irun=1,nruns - call variational_montecarlo(a,tau,nmax,X(irun),accep(irun)) - enddo - call ave_error(X,nruns,ave,err) - print *, 'E = ', ave, '+/-', err - call ave_error(accep,nruns,ave,err) - print *, 'A = ', ave, '+/-', err -end program qmc - #+END_SRC - - #+begin_src sh :results output :exports both -gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis -./vmc_metropolis - #+end_src - -**** Fortran :solution: + *Fortran* #+BEGIN_SRC f90 subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate) implicit none @@ -1794,7 +1804,7 @@ program qmc end program qmc #+END_SRC - #+begin_src sh :results output :exports both + #+begin_src sh :results output :exports results gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis ./vmc_metropolis #+end_src @@ -1810,135 +1820,158 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis :END: - Consider the time-dependent Schrödinger equation: - - \[ - i\frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \hat{H} \Psi(\mathbf{r},t) - \] - - We can expand $\Psi(\mathbf{r},0)$, in the basis of the eigenstates - of the time-independent Hamiltonian: - - \[ - \Psi(\mathbf{r},0) = \sum_k a_k\, \Phi_k(\mathbf{r}). - \] - - The solution of the Schrödinger equation at time $t$ is - - \[ - \Psi(\mathbf{r},t) = \sum_k a_k \exp \left( -i\, E_k\, t \right) \Phi_k(\mathbf{r}). - \] - - Now, let's replace the time variable $t$ by an imaginary time variable - $\tau=i\,t$, we obtain - - \[ - -\frac{\partial \psi(\mathbf{r}, \tau)}{\partial \tau} = \hat{H} \psi(\mathbf{r}, \tau) - \] - - where $\psi(\mathbf{r},\tau) = \Psi(\mathbf{r},-i\tau) = \Psi(\mathbf{r},t)$ - and - \[ - \psi(\mathbf{r},\tau) = \sum_k a_k \exp( -E_k\, \tau) \phi_k(\mathbf{r}). - \] - For large positive values of $\tau$, $\psi$ is dominated by the - $k=0$ term, namely the lowest eigenstate. - So we can expect that simulating the differetial equation in - imaginary time will converge to the exact ground state of the - system. - - - The [[https://en.wikipedia.org/wiki/Diffusion_equation][diffusion equation]] of particles is given by - - \[ - \frac{\partial \phi(\mathbf{r},t)}{\partial t} = D\, \Delta \phi(\mathbf{r},t). - \] - - The [[https://en.wikipedia.org/wiki/Reaction_rate][rate of reaction]] $v$ is the speed at which a chemical reaction - takes place. In a solution, the rate is given as a function of the - concentration $[A]$ by - - \[ - v = \frac{d[A]}{dt}, - \] - - where the concentration $[A]$ is proportional to the number of particles. - - These two equations allow us to interpret the Schrödinger equation - in imaginary time as the combination of: - - a diffusion equation with a diffusion coefficient $D=1/2$ for the - kinetic energy, and - - a rate equation for the potential. - - The diffusion equation can be simulated by a Brownian motion: - \[ \mathbf{r}_{n+1} = \mathbf{r}_{n} + \sqrt{\tau}\, \chi \] - where $\chi$ is a Gaussian random variable, and the rate equation - can be simulated by creating or destroying particles over time. - - /Diffusion Monte Carlo/ (DMC) consists in obtaining the ground state of a - system by simulating the Schrödinger equation in imaginary time, by - the combination of a diffusion process and a rate process. - In a molecular system, the potential is far from being constant, - and diverges at inter-particle coalescence points. Hence, when the - rate equation is simulated, it results in very large fluctuations - in the numbers of particles, 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 - determinant, CI wave function, /etc/), one obtains +** Schrödinger equation in imaginary time + + Consider the time-dependent Schrödinger equation: + \[ + i\frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \hat{H} \Psi(\mathbf{r},t) + \] - \begin{eqnarray*} - \end{eqnarray*} -\[ - -\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) = - \left[ -\frac{1}{2} \Delta \psi(\mathbf{r},\tau) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \right] \Psi_T(\mathbf{r}) -\] + We can expand $\Psi(\mathbf{r},0)$, in the basis of the eigenstates + of the time-independent Hamiltonian: -\[ --\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} -= -\frac{1}{2} \Big( \Delta \big[ -\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] - -\psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) - 2 -\nabla \psi(\mathbf{r},\tau) \nabla \Psi_T(\mathbf{r}) \Big) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) -\] + \[ + \Psi(\mathbf{r},0) = \sum_k a_k\, \Phi_k(\mathbf{r}). + \] -\[ --\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} -= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] + - \frac{1}{2} \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) + -\Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) -\] + The solution of the Schrödinger equation at time $t$ is -\[ --\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} -= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] + - \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) + -\Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) -\] -\[ --\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} -= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] + -\nabla \left[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) -\frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} -\right] + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) -\] + \[ + \Psi(\mathbf{r},t) = \sum_k a_k \exp \left( -i\, E_k\, t \right) \Phi_k(\mathbf{r}). + \] -Defining $\Pi(\mathbf{r},t) = \psi(\mathbf{r},\tau) -\Psi_T(\mathbf{r})$, + Now, let's replace the time variable $t$ by an imaginary time variable + $\tau=i\,t$, we obtain -\[ --\frac{\partial \Pi(\mathbf{r},\tau)}{\partial \tau} -= -\frac{1}{2} \Delta \Pi(\mathbf{r},\tau) + -\nabla \left[ \Pi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} -\right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau) -\] + \[ + -\frac{\partial \psi(\mathbf{r}, \tau)}{\partial \tau} = \hat{H} \psi(\mathbf{r}, \tau) + \] -The new "potential" is the local energy, which has smaller fluctuations -as $\Psi_T$ tends to the exact wave function. The new "kinetic energy" -can be simulated by the drifted diffusion scheme presented in the -previous section (VMC). + where $\psi(\mathbf{r},\tau) = \Psi(\mathbf{r},-i\tau) = \Psi(\mathbf{r},t)$ + and + \[ + \psi(\mathbf{r},\tau) = \sum_k a_k \exp( -E_k\, \tau) \phi_k(\mathbf{r}). + \] + For large positive values of $\tau$, $\psi$ is dominated by the + $k=0$ term, namely the lowest eigenstate. + So we can expect that simulating the differetial equation in + imaginary time will converge to the exact ground state of the + system. + + +** Diffusion and branching + + The [[https://en.wikipedia.org/wiki/Diffusion_equation][diffusion equation]] of particles is given by + + \[ + \frac{\partial \phi(\mathbf{r},t)}{\partial t} = D\, \Delta \phi(\mathbf{r},t). + \] + + The [[https://en.wikipedia.org/wiki/Reaction_rate][rate of reaction]] $v$ is the speed at which a chemical reaction + takes place. In a solution, the rate is given as a function of the + concentration $[A]$ by + + \[ + v = \frac{d[A]}{dt}, + \] + + where the concentration $[A]$ is proportional to the number of particles. + + These two equations allow us to interpret the Schrödinger equation + in imaginary time as the combination of: + - a diffusion equation with a diffusion coefficient $D=1/2$ for the + kinetic energy, and + - a rate equation for the potential. + + The diffusion equation can be simulated by a Brownian motion: + \[ \mathbf{r}_{n+1} = \mathbf{r}_{n} + \sqrt{\tau}\, \chi \] + where $\chi$ is a Gaussian random variable, and the rate equation + can be simulated by creating or destroying particles over time (a + so-called branching process). + + /Diffusion Monte Carlo/ (DMC) consists in obtaining the ground state of a + system by simulating the Schrödinger equation in imaginary time, by + the combination of a diffusion process and a branching process. + +** Importance sampling + + In a molecular system, the potential is far from being constant, + and diverges at inter-particle coalescence points. Hence, when the + rate equation is simulated, it results in very large fluctuations + in the numbers of particles, 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 + determinant, CI wave function, /etc/), one obtains + + \[ + -\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) = + \left[ -\frac{1}{2} \Delta \psi(\mathbf{r},\tau) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \right] \Psi_T(\mathbf{r}) + \] + + Defining $\Pi(\mathbf{r},t) = \psi(\mathbf{r},\tau) + \Psi_T(\mathbf{r})$, + + \[ + -\frac{\partial \Pi(\mathbf{r},\tau)}{\partial \tau} + = -\frac{1}{2} \Delta \Pi(\mathbf{r},\tau) + + \nabla \left[ \Pi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + \right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau) + \] + + The new "potential" is the local energy, which has smaller fluctuations + as $\Psi_T$ tends to the exact wave function. The new "kinetic energy" + can be simulated by the drifted diffusion scheme presented in the + previous section (VMC). + + +*** Appendix : Details of the Derivation + + \[ + -\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) = + \left[ -\frac{1}{2} \Delta \psi(\mathbf{r},\tau) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \right] \Psi_T(\mathbf{r}) + \] + + \[ + -\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} + = -\frac{1}{2} \Big( \Delta \big[ + \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] - + \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) - 2 + \nabla \psi(\mathbf{r},\tau) \nabla \Psi_T(\mathbf{r}) \Big) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) + \] + + \[ + -\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} + = -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] + + \frac{1}{2} \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) + + \Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) + \] + + \[ + -\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} + = -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] + + \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) + + \Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) + \] + \[ + -\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau} + = -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] + + \nabla \left[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) + \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + \right] + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) + \] + + Defining $\Pi(\mathbf{r},t) = \psi(\mathbf{r},\tau) + \Psi_T(\mathbf{r})$, + + \[ + -\frac{\partial \Pi(\mathbf{r},\tau)}{\partial \tau} + = -\frac{1}{2} \Delta \Pi(\mathbf{r},\tau) + + \nabla \left[ \Pi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + \right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau) + \] ** TODO Hydrogen atom diff --git a/docs/worg.css b/docs/worg.css index c9725e1..f7fd150 100644 --- a/docs/worg.css +++ b/docs/worg.css @@ -156,7 +156,8 @@ } .tag { - + background-color: #ffff; + color: #ffff; } li {