qmc-lttc/QMC.org

#+TITLE: Quantum Monte Carlo
#+AUTHOR: Anthony Scemama, Claudia Filippi
#+LANGUAGE:  en
#+INFOJS_OPT: toc:t mouse:underline path:org-info.js
#+STARTUP: latexpreview
#+LATEX_CLASS: report
#+LATEX_HEADER_EXTRA: \usepackage{minted}
#+HTML_HEAD: <link rel="stylesheet" title="Standard" href="worg.css" type="text/css" />

#+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: solution

  #+BEGIN_SRC elisp :output none :exports none
(setq org-latex-listings 'minted
      org-latex-packages-alist '(("" "minted"))
      org-latex-pdf-process
      '("pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f"
        "pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f"
        "pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f"))
(setq org-latex-minted-options '(("breaklines" "true")
                                 ("breakanywhere" "true")))
(setq org-latex-minted-options
      '(("frame" "lines")
        ("fontsize" "\\scriptsize")
        ("linenos" "")))
(org-beamer-export-to-pdf)
                            
  #+END_SRC   

  #+RESULTS:
  : /home/scemama/TREX/qmc-lttc/QMC.pdf

* 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
  to introduce the /local energy/.
  Then we introduce the variational Monte Carlo (VMC) method which
  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 hydrogen atom and of the H_2 molecule. 

  Code examples will be given in Python and Fortran. You can use
  whatever language you prefer to write the program.

  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.

  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
  wave function:

  $$
  \Psi(\mathbf{r}) = \exp(-a |\mathbf{r}|)
  $$

  We will first verify that, for a given value of $a$, $\Psi$ is an
  eigenfunction of the Hamiltonian

  $$
  \hat{H} = \hat{T} + \hat{V} = - \frac{1}{2} \Delta - \frac{1}{|\mathbf{r}|}
  $$

  To do that, we will check if the local energy, defined as

  $$
  E_L(\mathbf{r}) = \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})},
  $$

  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
   :header-args:f90: :tangle hydrogen.f90
   :END:

   Write all the functions of this section in a single file :
 ~hydrogen.py~ if you use Python, or ~hydrogen.f90~ is you use
   Fortran.
   
   #+begin_note
   - When computing a square root in $\mathbb{R}$, *always* make sure
     that the argument of the square root is non-negative.
   - When you divide, *always* make sure that you will not divide by zero

   If a /floating-point exception/ can occur, you should make a test
   to catch the error.
   #+end_note
   
*** Exercise 1

    #+begin_exercise
    Write a function which computes the potential at $\mathbf{r}$.
    The function accepts a 3-dimensional vector =r= as input arguments
    and returns the potential.
    #+end_exercise

    $\mathbf{r}=\left( \begin{array}{c} x \\ y\\ z\end{array} \right)$, so
    $$
    V(\mathbf{r}) = -\frac{1}{\sqrt{x^2 + y^2 + z^2}}
    $$

 *Python*
     #+BEGIN_SRC python :results none :tangle none
import numpy as np

def potential(r):
    # TODO
     #+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

**** Solution                                                      :solution:
 *Python*
     #+BEGIN_SRC python :results none
import numpy as np

def potential(r):
    distance = np.sqrt(np.dot(r,r))
    assert (distance > 0)
    return -1. / distance
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90 
double precision function potential(r)
  implicit none
  double precision, intent(in) :: r(3)

  double precision             :: distance

  distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )

  if (distance > 0.d0) then
     potential = -1.d0 / distance
  else
     stop 'potential at r=0.d0 diverges'
  end if

end function potential
     #+END_SRC

*** 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*
     #+BEGIN_SRC python :results none  :tangle none
def psi(a, r):
    # TODO
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90  :tangle none
double precision function psi(a, r)
  implicit none
  double precision, intent(in) :: a, r(3)

  ! TODO

end function psi
     #+END_SRC
     
**** 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
  double precision, intent(in) :: a, r(3)

  psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
end function psi
     #+END_SRC
     
*** 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
    local kinetic energy.
    #+end_exercise

    The local kinetic energy is defined as $$-\frac{1}{2}\frac{\Delta \Psi}{\Psi}.$$
     
    We differentiate $\Psi$ with respect to $x$:
     
    \[\Psi(\mathbf{r})  =  \exp(-a\,|\mathbf{r}|) \]
    \[\frac{\partial \Psi}{\partial x}
      = \frac{\partial \Psi}{\partial |\mathbf{r}|} \frac{\partial |\mathbf{r}|}{\partial x}   
      =  - \frac{a\,x}{|\mathbf{r}|} \Psi(\mathbf{r}) \]

    and we differentiate a second time:

    $$
    \frac{\partial^2 \Psi}{\partial x^2} =
    \left( \frac{a^2\,x^2}{|\mathbf{r}|^2}  -
    \frac{a(y^2+z^2)}{|\mathbf{r}|^{3}} \right) \Psi(\mathbf{r}).
    $$

    The Laplacian operator $\Delta = \frac{\partial^2}{\partial x^2} +
    \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$
    applied to the wave function gives:

    $$
    \Delta \Psi (\mathbf{r}) = \left(a^2 - \frac{2a}{\mathbf{|r|}} \right) \Psi(\mathbf{r})
    $$

    So the local kinetic energy is
    $$
    -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) 
    $$
     
 *Python*
     #+BEGIN_SRC python :results none :tangle none
def kinetic(a,r):
    # TODO
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90  :tangle none
double precision function kinetic(a,r)
  implicit none
  double precision, intent(in) :: a, r(3)

  ! TODO

end function kinetic
     #+END_SRC

**** Solution                                                      :solution:
 *Python*
     #+BEGIN_SRC python :results none
def kinetic(a,r):
    distance = np.sqrt(np.dot(r,r))
    assert (distance > 0.)

    return a * (1./distance - 0.5 * a)
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90 
double precision function kinetic(a,r)
  implicit none
  double precision, intent(in) :: a, r(3)

  double precision             :: distance

  distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) 

  if (distance > 0.d0) then

     kinetic =  a * (1.d0 / distance - 0.5d0 * a)

  else
     stop 'kinetic energy diverges at r=0'
  end if

end function kinetic
     #+END_SRC

*** Exercise 4
    #+begin_exercise
    Write a function which computes the local energy at $\mathbf{r}$,
    using the previously defined functions.
    The function accepts =a= and =r= as input arguments and returns the
    local kinetic energy.
    #+end_exercise
   
    $$
    E_L(\mathbf{r}) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) + V(\mathbf{r})
    $$

    
 *Python*
     #+BEGIN_SRC python :results none :tangle none
def e_loc(a,r):
    #TODO
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90 :tangle none
double precision function e_loc(a,r)
  implicit none
  double precision, intent(in) :: a, r(3)

  ! TODO

end function e_loc
     #+END_SRC
   
**** 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
  double precision, intent(in) :: a, r(3)

  double precision, external   :: kinetic, potential

  e_loc = kinetic(a,r) + potential(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
   :header-args:f90: :tangle plot_hydrogen.f90
   :END:
   
   #+begin_note
   The potential and the kinetic energy both diverge at $r=0$, so we
   choose a grid which does not contain the origin.
   #+end_note

*** Exercise
    #+begin_exercise
    For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the
    local energy along the $x$ axis. In Python, you can use matplotlib
    for example. In Fortran, it is convenient to write in a text file
    the values of $x$ and $E_L(\mathbf{r})$ for each point, and use
    Gnuplot to plot the files.
    #+end_exercise

 *Python*
     #+BEGIN_SRC python :results none :tangle 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))

# TODO

plt.tight_layout()
plt.legend()
plt.savefig("plot_py.png")
     #+end_src

 *Fortran*
     #+begin_src f90  :tangle none
program plot
  implicit none
  double precision, external :: e_loc

  double precision :: x(50), dx
  integer :: i, j

  dx = 10.d0/(size(x)-1)
  do i=1,size(x)
     x(i) = -5.d0 + (i-1)*dx
  end do

  ! TODO

end program plot
     #+end_src

     To compile and run:

     #+begin_src sh :exports both
gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
./plot_hydrogen > data
     #+end_src

     To plot the data using Gnuplot:

     #+begin_src gnuplot :file plot.png :exports code
set grid
set xrange [-5:5]
set yrange [-2:1]
plot './data' index 0 using 1:2 with lines title 'a=0.1', \
     './data' index 1 using 1:2 with lines title 'a=0.2', \
     './data' index 2 using 1:2 with lines title 'a=0.5', \
     './data' index 3 using 1:2 with lines title 'a=1.0', \
     './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

**** 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
  double precision, external :: e_loc

  double precision :: x(50), energy, dx, r(3), a(6)
  integer :: i, 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

  r(:) = 0.d0

  do j=1,size(a)
     print *, '# a=', a(j)
     do i=1,size(x)
        r(1) = x(i)
        energy = e_loc( a(j), r )
        print *, x(i), energy
     end do
     print *, ''
     print *, ''
  end do

end program plot
     #+end_src

     #+begin_src sh :exports none
gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
./plot_hydrogen > data
     #+end_src

     #+begin_src gnuplot :file plot.png :exports results
set grid
set xrange [-5:5]
set yrange [-2:1]
plot './data' index 0 using 1:2 with lines title 'a=0.1', \
     './data' index 1 using 1:2 with lines title 'a=0.2', \
     './data' index 2 using 1:2 with lines title 'a=0.5', \
     './data' index 3 using 1:2 with lines title 'a=1.0', \
     './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]]

** Numerical estimation of the energy
   :PROPERTIES:
   :header-args:python: :tangle energy_hydrogen.py
   :header-args:f90: :tangle energy_hydrogen.f90
   :END:

   If the space is discretized in small volume elements $\mathbf{r}_i$
   of size $\delta \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}_i$
   multiplied by the volume element:
     
   $$
   \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}
   $$
     
   #+begin_note
   The energy is biased because:
   - The volume elements are not infinitely small (discretization error)
   - The energy is evaluated only inside the box (incompleteness of the space)
   #+end_note

   
*** Exercise
     #+begin_exercise
    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)$.
     #+end_exercise

 *Python*
     #+BEGIN_SRC python :results none :tangle none
import numpy as np
from hydrogen import e_loc, psi

interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3

r = np.array([0.,0.,0.])

for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
    # TODO
    print(f"a = {a} \t E = {E}")                

     #+end_src

 *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

interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3

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}")                

     #+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

 *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

  delta = dx**3

  r(:) = 0.d0

  do j=1,size(a)
     energy = 0.d0
     norm   = 0.d0
     
     do i=1,size(x)
        r(1) = x(i)

        do k=1,size(x)
           r(2) = x(k)

           do l=1,size(x)
              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

        end do

     end do

     energy = energy / norm
     print *, 'a = ', a(j), '    E = ', energy
  end do

end program energy_hydrogen
     #+end_src

     #+begin_src sh :results output :exports results
gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
./energy_hydrogen 
     #+end_src

     #+RESULTS:
     :  a =   0.10000000000000001          E =  -0.24518438948809140     
     :  a =   0.20000000000000001          E =  -0.26966057967803236     
     :  a =   0.50000000000000000          E =  -0.38563576125173815     
     :  a =    1.0000000000000000          E =  -0.50000000000000000     
     :  a =    1.5000000000000000          E =  -0.39242967082602065     
     :  a =    2.0000000000000000          E =   -8.0869806678448772E-002

** Variance of the local energy
   :PROPERTIES:
   :header-args:python: :tangle variance_hydrogen.py
   :header-args:f90: :tangle variance_hydrogen.f90
   :END:

   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 its 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}}
   $$
   which can be simplified as
   
   $$ \sigma^2(E_L) = \langle E_L^2 \rangle_{\Psi^2} - \langle E_L \rangle_{\Psi^2}^2.$$

   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.

*** Exercise (optional)
   #+begin_exercise
   Prove that :
   $$\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
   
**** Solution                                                      :solution:

   $\bar{E} = \langle E \rangle$ is a constant, so $\langle \bar{E}
   \rangle = \bar{E}$ .
   
   \begin{eqnarray*}
   \langle E - \bar{E} \rangle^2 & = & 
   \langle E^2 - 2 E \bar{E} + \bar{E}^2 \rangle \\
   &=& \langle E^2 \rangle - 2 \langle E \bar{E} \rangle + \langle \bar{E}^2 \rangle \\
   &=& \langle E^2 \rangle - 2 \langle E \rangle \bar{E}  + \bar{E}^2 \\
   &=& \langle E^2 \rangle - 2 \bar{E}^2  + \bar{E}^2 \\
   &=& \langle E^2 \rangle - \bar{E}^2 \\
   &=& \langle E^2 \rangle - \langle E \rangle^2 \\
   \end{eqnarray*}
*** Exercise
   #+begin_exercise
   Add the calculation of the variance to the previous code, and
   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)$ for different values of $a$.
   #+end_exercise
     
 *Python*
     #+begin_src python :results none :tangle none
import numpy as np from hydrogen import e_loc, psi

interval = np.linspace(-5,5,num=50)

delta = (interval[1]-interval[0])**3

r = np.array([0.,0.,0.])

for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:

    # TODO

    print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")
    #+end_src

 *Fortran*
     #+begin_src f90 :tangle none
program variance_hydrogen
  implicit none

  double precision :: x(50), w, delta, energy, energy2
  double precision :: dx, r(3), a(6), norm, e_tmp, s2
  integer          :: i, k, l, j

  double precision, external :: e_loc, psi

  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 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*