scemama/qmc-lttc
1
0
Fork 1
Code Issues Pull Requests Projects Releases Wiki Activity

Added style

Browse Source
This commit is contained in:
Anthony Scemama 2021-01-11 20:54:40 +01:00
parent cb6ceeb797
commit 277e243545
2 changed files with 1110 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

287
QMC.org
View File

@ -4,9 +4,7 @@
# SETUPFILE: https://fniessen.github.io/org-html-themes/org/theme-bigblow.setup # SETUPFILE: https://fniessen.github.io/org-html-themes/org/theme-bigblow.setup
#+STARTUP: latexpreview #+STARTUP: latexpreview
#+HTML_HEAD: <link rel="stylesheet" title="Standard" href="https://orgmode.org/worg/style/worg.css" type="text/css" /> #+HTML_HEAD: <link rel="stylesheet" title="Standard" href="worg.css" type="text/css" />
#+HTML_HEAD: <link rel="alternate stylesheet" title="Zenburn" href="https://orgmode.org/worg/style/worg-zenburn.css" type="text/css" />
#+HTML_HEAD: <link rel="alternate stylesheet" title="Classic" href="https://orgmode.org/worg/style/worg-classic.css" type="text/css" />
* Introduction * Introduction
@ -30,15 +28,12 @@
$\Psi : \mathbb{R}^{3N} \rightarrow \mathbb{R}$. In addition, $\Psi$ $\Psi : \mathbb{R}^{3N} \rightarrow \mathbb{R}$. In addition, $\Psi$
is defined everywhere, continuous and infinitely differentiable. is defined everywhere, continuous and infinitely differentiable.
** Python *Note*
#+begin_important
** Fortran In Fortran, when you use a double precision constant, don't forget
to put d0 as a suffix (for example 2.0d0), or it will be
- 1.d0 interpreted as a single precision value
- external #+end_important
- r(:) = 0.d0
- a = (/ 0.1, 0.2 /)
- size(x)
* Numerical evaluation of the energy * Numerical evaluation of the energy
@ -63,9 +58,9 @@
E_L(\mathbf{r}) = \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}, E_L(\mathbf{r}) = \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})},
$$ $$
is constant. is constant. We will also see that when $a \ne 1$ the local energy
is not constant, so $\hat{H} \Psi \ne E \Psi$.
The probabilistic /expected value/ of an arbitrary function $f(x)$ The probabilistic /expected value/ of an arbitrary function $f(x)$
with respect to a probability density function $p(x)$ is given by with respect to a probability density function $p(x)$ is given by
@ -79,16 +74,16 @@
The electronic energy of a system is the expectation value of the The electronic energy of a system is the expectation value of the
local energy $E(\mathbf{r})$ with respect to the $3N$-dimensional local energy $E(\mathbf{r})$ with respect to the 3N-dimensional
electron density given by the square of the wave function: electron density given by the square of the wave function:
\begin{eqnarray} \begin{eqnarray*}
E & = & \frac{\langle \Psi| \hat{H} | \Psi\rangle}{\langle \Psi |\Psi \rangle} \\ 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 \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\, \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}} = \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} = \langle E_L \rangle_{\Psi^2}
\end{eqnarray} \end{eqnarray*}
** Local energy ** Local energy
:PROPERTIES: :PROPERTIES:
@ -99,67 +94,68 @@
The function accepts a 3-dimensional vector =r= as input arguments The function accepts a 3-dimensional vector =r= as input arguments
and returns the potential. and returns the potential.
$\mathbf{r}=\sqrt{x^2 + y^2 + z^2}$, so $\mathbf{r}=\left( \begin{array}{c} x \\ y\\ z\end{array} \right)$, so
$$ $$
V(x,y,z) = -\frac{1}{\sqrt{x^2 + y^2 + z^2}} V(\mathbf{r}) = -\frac{1}{\sqrt{x^2 + y^2 + z^2}}
$$ $$
#+BEGIN_SRC python :results none *Python*
#+BEGIN_SRC python :results none
import numpy as np import numpy as np
def potential(r): def potential(r):
return -1. / np.sqrt(np.dot(r,r)) return -1. / np.sqrt(np.dot(r,r))
#+END_SRC #+END_SRC
#+BEGIN_SRC f90 *Fortran*
#+BEGIN_SRC f90
double precision function potential(r) double precision function potential(r)
implicit none implicit none
double precision, intent(in) :: r(3) double precision, intent(in) :: r(3)
potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
end function potential end function potential
#+END_SRC #+END_SRC
*** Write a function which computes the wave function at $\mathbf{r}$ *** Write a function which computes the wave function at $\mathbf{r}$
The function accepts a scalar =a= and a 3-dimensional vector =r= as The function accepts a scalar =a= and a 3-dimensional vector =r= as
input arguments, and returns a scalar. input arguments, and returns a scalar.
#+BEGIN_SRC python :results none
*Python*
#+BEGIN_SRC python :results none
def psi(a, r): def psi(a, r):
return np.exp(-a*np.sqrt(np.dot(r,r))) return np.exp(-a*np.sqrt(np.dot(r,r)))
#+END_SRC #+END_SRC
#+BEGIN_SRC f90 *Fortran*
#+BEGIN_SRC f90
double precision function psi(a, r) double precision function psi(a, r)
implicit none implicit none
double precision, intent(in) :: a, r(3) double precision, intent(in) :: a, r(3)
psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )) psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
end function psi end function psi
#+END_SRC #+END_SRC
*** Write a function which computes the local kinetic energy at $\mathbf{r}$ *** Write a function which computes the local kinetic energy at $\mathbf{r}$
The function accepts =a= and =r= as input arguments and returns the The function accepts =a= and =r= as input arguments and returns the
local kinetic energy. local kinetic energy.
The local kinetic energy is defined as $$-\frac{1}{2}\frac{\Delta \Psi}{\Psi}$$. The local kinetic energy is defined as $$-\frac{1}{2}\frac{\Delta \Psi}{\Psi}.$$
$$
\Psi(x,y,z) = \exp(-a\,\sqrt{x^2 + y^2 + z^2}).
$$
We differentiate $\Psi$ with respect to $x$: We differentiate $\Psi$ with respect to $x$:
$$ \[\Psi(\mathbf{r}) = \exp(-a\,|\mathbf{r}|) \]
\frac{\partial \Psi}{\partial x} \[\frac{\partial \Psi}{\partial x}
= \frac{\partial \Psi}{\partial r} \frac{\partial r}{\partial x} = \frac{\partial \Psi}{\partial |\mathbf{r}|} \frac{\partial |\mathbf{r}|}{\partial x}
= - \frac{a\,x}{|\mathbf{r}|} \Psi(x,y,z) = - \frac{a\,x}{|\mathbf{r}|} \Psi(\mathbf{r}) \]
$$
and we differentiate a second time: and we differentiate a second time:
$$ $$
\frac{\partial^2 \Psi}{\partial x^2} = \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(x,y,z). \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} + The Laplacian operator $\Delta = \frac{\partial^2}{\partial x^2} +
@ -167,49 +163,54 @@ end function psi
applied to the wave function gives: applied to the wave function gives:
$$ $$
\Delta \Psi (x,y,z) = \left(a^2 - \frac{2a}{\mathbf{|r|}} \right) \Psi(x,y,z) \Delta \Psi (\mathbf{r}) = \left(a^2 - \frac{2a}{\mathbf{|r|}} \right) \Psi(\mathbf{r})
$$ $$
So the local kinetic energy is So the local kinetic energy is
$$ $$
-\frac{1}{2} \frac{\Delta \Psi}{\Psi} (x,y,z) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right)
$$ $$
#+BEGIN_SRC python :results none *Python*
#+BEGIN_SRC python :results none
def kinetic(a,r): def kinetic(a,r):
return -0.5 * (a**2 - (2.*a)/np.sqrt(np.dot(r,r))) return -0.5 * (a**2 - (2.*a)/np.sqrt(np.dot(r,r)))
#+END_SRC #+END_SRC
#+BEGIN_SRC f90 *Fortran*
#+BEGIN_SRC f90
double precision function kinetic(a,r) double precision function kinetic(a,r)
implicit none implicit none
double precision, intent(in) :: a, r(3) double precision, intent(in) :: a, r(3)
kinetic = -0.5d0 * (a*a - (2.d0*a) / & kinetic = -0.5d0 * (a*a - (2.d0*a) / &
dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) ) dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) )
end function kinetic end function kinetic
#+END_SRC #+END_SRC
*** Write a function which computes the local energy at $\mathbf{r}$ *** Write a function which computes the local energy at $\mathbf{r}$
The function accepts =x,y,z= as input arguments and returns the The function accepts =x,y,z= as input arguments and returns the
local energy. local energy.
$$ $$
E_L(x,y,z) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (x,y,z) + V(x,y,z) E_L(\mathbf{r}) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) + V(\mathbf{r})
$$ $$
#+BEGIN_SRC python :results none
*Python*
#+BEGIN_SRC python :results none
def e_loc(a,r): def e_loc(a,r):
return kinetic(a,r) + potential(r) return kinetic(a,r) + potential(r)
#+END_SRC #+END_SRC
#+BEGIN_SRC f90 *Fortran*
#+BEGIN_SRC f90
double precision function e_loc(a,r) double precision function e_loc(a,r)
implicit none implicit none
double precision, intent(in) :: a, r(3) double precision, intent(in) :: a, r(3)
double precision, external :: kinetic, potential double precision, external :: kinetic, potential
e_loc = kinetic(a,r) + potential(r) e_loc = kinetic(a,r) + potential(r)
end function e_loc end function e_loc
#+END_SRC #+END_SRC
** Plot of the local energy along the $x$ axis ** Plot of the local energy along the $x$ axis
:PROPERTIES: :PROPERTIES:
@ -220,6 +221,7 @@ end function e_loc
For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the
local energy along the $x$ axis. local energy along the $x$ axis.
*Python*
#+BEGIN_SRC python :results none #+BEGIN_SRC python :results none
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
@ -249,6 +251,8 @@ plt.savefig("plot_py.png")
[[./plot_py.png]] [[./plot_py.png]]
*Fortran*
#+begin_src f90 #+begin_src f90
program plot program plot
implicit none implicit none
@ -289,7 +293,7 @@ gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
#+RESULTS: #+RESULTS:
To plot the data using gnuplot" To plot the data using gnuplot:
#+begin_src gnuplot :file plot.png :exports both #+begin_src gnuplot :file plot.png :exports both
set grid set grid
@ -312,11 +316,11 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
:header-args:f90: :tangle energy_hydrogen.f90 :header-args:f90: :tangle energy_hydrogen.f90
:END: :END:
If the space is discretized in small volume elements $\delta If the space is discretized in small volume elements $\mathbf{r}_i$
\mathbf{r}$, the expression of \langle E_L \rangle_{\Psi^2}$ becomes of size $\delta \mathbf{r}$, the expression of $\langle E_L \rangle_{\Psi^2}$
a weighted average of the local energy, where the weights are the becomes a weighted average of the local energy, where the weights
values of the probability density at $\mathbf{r}$ multiplied are the values of the probability density at $\mathbf{r}_i$
by the volume element: 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}, \;\; \langle E \rangle_{\Psi^2} \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\;
@ -327,11 +331,14 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
energy in a grid of $50\times50\times50$ points in the range energy in a grid of $50\times50\times50$ points in the range
$(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$. $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$.
Note: the energy is biased because: #+begin_note
The energy is biased because:
- The volume elements are not infinitely small (discretization error) - The volume elements are not infinitely small (discretization error)
- The energy is evaluated only inside the box (incompleteness of the space) - The energy is evaluated only inside the box (incompleteness of the space)
#+end_note
#+BEGIN_SRC python :results none *Python*
#+BEGIN_SRC python :results none
import numpy as np import numpy as np
from hydrogen import e_loc, psi from hydrogen import e_loc, psi
@ -356,19 +363,19 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
E = E / norm E = E / norm
print(f"a = {a} \t E = {E}") print(f"a = {a} \t E = {E}")
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: a = 0.1 E = -0.24518438948809218 : a = 0.1 E = -0.24518438948809218
: a = 0.2 E = -0.26966057967803525 : a = 0.2 E = -0.26966057967803525
: a = 0.5 E = -0.3856357612517407 : a = 0.5 E = -0.3856357612517407
: a = 0.9 E = -0.49435709786716214 : a = 0.9 E = -0.49435709786716214
: a = 1.0 E = -0.5 : a = 1.0 E = -0.5
: a = 1.5 E = -0.39242967082602226 : a = 1.5 E = -0.39242967082602226
: a = 2.0 E = -0.08086980667844901 : a = 2.0 E = -0.08086980667844901
*Fortran*
#+begin_src f90 #+begin_src f90
program energy_hydrogen program energy_hydrogen
implicit none implicit none
double precision, external :: e_loc, psi double precision, external :: e_loc, psi
@ -407,22 +414,22 @@ program energy_hydrogen
end do end do
end program energy_hydrogen end program energy_hydrogen
#+end_src #+end_src
To compile the Fortran and run it: To compile the Fortran and run it:
#+begin_src sh :results output :exports both #+begin_src sh :results output :exports both
gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
./energy_hydrogen ./energy_hydrogen
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: a = 0.10000000000000001 E = -0.24518438948809140 : a = 0.10000000000000001 E = -0.24518438948809140
: a = 0.20000000000000001 E = -0.26966057967803236 : a = 0.20000000000000001 E = -0.26966057967803236
: a = 0.50000000000000000 E = -0.38563576125173815 : a = 0.50000000000000000 E = -0.38563576125173815
: a = 1.0000000000000000 E = -0.50000000000000000 : a = 1.0000000000000000 E = -0.50000000000000000
: a = 1.5000000000000000 E = -0.39242967082602065 : a = 1.5000000000000000 E = -0.39242967082602065
: a = 2.0000000000000000 E = -8.0869806678448772E-002 : a = 2.0000000000000000 E = -8.0869806678448772E-002
** Compute the variance of the local energy ** Compute the variance of the local energy
:PROPERTIES: :PROPERTIES:
@ -446,6 +453,7 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
Compute a numerical estimate of the variance of the local energy 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)$. in a grid of $50\times50\times50$ points in the range $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$.
*Python*
#+begin_src python :results none #+begin_src python :results none
import numpy as np import numpy as np
from hydrogen import e_loc, psi from hydrogen import e_loc, psi
@ -469,8 +477,8 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
El = e_loc(a, r) El = e_loc(a, r)
E += w * El E += w * El
norm += w norm += w
E = E / norm E = E / norm
s2 = 0. s2 = 0.
for x in interval: for x in interval:
r[0] = x r[0] = x
for y in interval: for y in interval:
@ -481,8 +489,8 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
w = w * w * delta w = w * w * delta
El = e_loc(a, r) El = e_loc(a, r)
s2 += w * (El - E)**2 s2 += w * (El - E)**2
s2 = s2 / norm s2 = s2 / norm
print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}") print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")
#+end_src #+end_src
#+RESULTS: #+RESULTS:
@ -494,6 +502,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 = 1.5 E = -0.39242967 \sigma^2 = 0.31449671
: a = 2.0 E = -0.08086981 \sigma^2 = 1.80688143 : a = 2.0 E = -0.08086981 \sigma^2 = 1.80688143
*Fortran*
#+begin_src f90 #+begin_src f90
program variance_hydrogen program variance_hydrogen
implicit none implicit none
@ -523,7 +532,6 @@ program variance_hydrogen
r(3) = x(l) r(3) = x(l)
w = psi(a(j),r) w = psi(a(j),r)
w = w * w * delta w = w * w * delta
energy = energy + w * e_loc(a(j), r) energy = energy + w * e_loc(a(j), r)
norm = norm + w norm = norm + w
end do end do
@ -541,7 +549,6 @@ program variance_hydrogen
r(3) = x(l) r(3) = x(l)
w = psi(a(j),r) w = psi(a(j),r)
w = w * w * delta w = w * w * delta
s2 = s2 + w * ( e_loc(a(j), r) - energy )**2 s2 = s2 + w * ( e_loc(a(j), r) - energy )**2
norm = norm + w norm = norm + w
end do end do
@ -573,11 +580,9 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
* Variational Monte Carlo * Variational Monte Carlo
Numerical integration with deterministic methods is very efficient Numerical integration with deterministic methods is very efficient
in low dimensions. When the number of dimensions becomes larger than in low dimensions. When the number of dimensions becomes large,
Instead of computing the average energy as a numerical integration instead of computing the average energy as a numerical integration
on a grid, we will do a Monte Carlo sampling, which is an extremely on a grid, it is usually more efficient to do a Monte Carlo sampling.
efficient method to compute integrals when the number of dimensions is
large.
Moreover, a Monte Carlo sampling will alow us to remove the bias due Moreover, a Monte Carlo sampling will alow us to remove the bias due
to the discretization of space, and compute a statistical confidence to the discretization of space, and compute a statistical confidence
@ -615,6 +620,7 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
Write a function returning the average and statistical error of an Write a function returning the average and statistical error of an
input array. input array.
*Python*
#+BEGIN_SRC python :results none #+BEGIN_SRC python :results none
from math import sqrt from math import sqrt
def ave_error(arr): def ave_error(arr):
@ -625,6 +631,7 @@ def ave_error(arr):
return (average, sqrt(variance/M)) return (average, sqrt(variance/M))
#+END_SRC #+END_SRC
*Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
subroutine ave_error(x,n,ave,err) subroutine ave_error(x,n,ave,err)
implicit none implicit none
@ -667,6 +674,7 @@ end subroutine ave_error
Compute the energy of the wave function with $a=0.9$. Compute the energy of the wave function with $a=0.9$.
*Python*
#+BEGIN_SRC python :results output #+BEGIN_SRC python :results output
from hydrogen import * from hydrogen import *
from qmc_stats import * from qmc_stats import *
@ -692,6 +700,7 @@ print(f"E = {E} +/- {deltaE}")
#+RESULTS: #+RESULTS:
: E = -0.4956255109300764 +/- 0.0007082875482711226 : E = -0.4956255109300764 +/- 0.0007082875482711226
*Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
subroutine uniform_montecarlo(a,nmax,energy) subroutine uniform_montecarlo(a,nmax,energy)
implicit none implicit none
@ -764,6 +773,7 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
z_2 &=& \sqrt{-2 \ln u_1} \sin(2 \pi u_2) z_2 &=& \sqrt{-2 \ln u_1} \sin(2 \pi u_2)
\end{eqnarray*} \end{eqnarray*}
*Fortran*
#+BEGIN_SRC f90 :tangle qmc_stats.f90 #+BEGIN_SRC f90 :tangle qmc_stats.f90
subroutine random_gauss(z,n) subroutine random_gauss(z,n)
implicit none implicit none
@ -813,6 +823,7 @@ end subroutine random_gauss
w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r}_i)} \delta \mathbf{r} w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r}_i)} \delta \mathbf{r}
$$ $$
*Python*
#+BEGIN_SRC python :results output #+BEGIN_SRC python :results output
from hydrogen import * from hydrogen import *
from qmc_stats import * from qmc_stats import *
@ -843,6 +854,7 @@ print(f"E = {E} +/- {deltaE}")
: E = -0.49507506093129827 +/- 0.00014164037765553668 : E = -0.49507506093129827 +/- 0.00014164037765553668
*Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
double precision function gaussian(r) double precision function gaussian(r)
implicit none implicit none
@ -994,12 +1006,14 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
First, write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$. First, write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$.
*Python*
#+BEGIN_SRC python #+BEGIN_SRC python
def drift(a,r): def drift(a,r):
ar_inv = -a/np.sqrt(np.dot(r,r)) ar_inv = -a/np.sqrt(np.dot(r,r))
return r * ar_inv return r * ar_inv
#+END_SRC #+END_SRC
*Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
subroutine drift(a,r,b) subroutine drift(a,r,b)
implicit none implicit none
@ -1012,7 +1026,50 @@ end subroutine drift
#+END_SRC #+END_SRC
Now we can write the Monte Carlo sampling Now we can write the Monte Carlo sampling:
*Python*
#+BEGIN_SRC python
def MonteCarlo(a,tau,nmax):
E = 0.
N = 0.
sq_tau = sqrt(tau)
r_old = np.random.normal(loc=0., scale=1.0, size=(3))
d_old = drift(a,r_old)
d2_old = np.dot(d_old,d_old)
psi_old = psi(a,r_old)
for istep in range(nmax):
eta = np.random.normal(loc=0., scale=1.0, size=(3))
r_new = r_old + tau * d_old + sq_tau * eta
d_new = drift(a,r_new)
d2_new = np.dot(d_new,d_new)
psi_new = psi(a,r_new)
# Metropolis
prod = np.dot((d_new + d_old), (r_new - r_old))
argexpo = 0.5 * (d2_new - d2_old)*tau + prod
q = psi_new / psi_old
q = np.exp(-argexpo) * q*q
if np.random.uniform() < q:
r_old = r_new
d_old = d_new
d2_old = d2_new
psi_old = psi_new
N += 1.
E += e_loc(a,r_old)
return E/N
nmax = 100000
tau = 0.1
X = [MonteCarlo(a,tau,nmax) for i in range(30)]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
#+END_SRC
#+RESULTS:
: E = -0.4951783346213532 +/- 0.00022067316984271938
*Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
subroutine variational_montecarlo(a,nmax,energy) subroutine variational_montecarlo(a,nmax,energy)
implicit none implicit none
@ -1061,46 +1118,6 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
./vmc ./vmc
#+end_src #+end_src
#+BEGIN_SRC python
def MonteCarlo(a,tau,nmax):
E = 0.
N = 0.
sq_tau = sqrt(tau)
r_old = np.random.normal(loc=0., scale=1.0, size=(3))
d_old = drift(a,r_old)
d2_old = np.dot(d_old,d_old)
psi_old = psi(a,r_old)
for istep in range(nmax):
eta = np.random.normal(loc=0., scale=1.0, size=(3))
r_new = r_old + tau * d_old + sq_tau * eta
d_new = drift(a,r_new)
d2_new = np.dot(d_new,d_new)
psi_new = psi(a,r_new)
# Metropolis
prod = np.dot((d_new + d_old), (r_new - r_old))
argexpo = 0.5 * (d2_new - d2_old)*tau + prod
q = psi_new / psi_old
q = np.exp(-argexpo) * q*q
if np.random.uniform() < q:
r_old = r_new
d_old = d_new
d2_old = d2_new
psi_old = psi_new
N += 1.
E += e_loc(a,r_old)
return E/N
nmax = 100000
tau = 0.1
X = [MonteCarlo(a,tau,nmax) for i in range(30)]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
#+END_SRC
#+RESULTS:
: E = -0.4951783346213532 +/- 0.00022067316984271938
* Diffusion Monte Carlo * Diffusion Monte Carlo

958
worg.css Normal file
View File

@ -0,0 +1,958 @@
@import url(https://fonts.googleapis.com/css?family=Droid+Sans|Droid+Sans+Mono|Droid+Serif);
@media all
{
html {
margin: 0;
font: .9em/1.6em "Droid Serif", Cambria, Georgia, "DejaVu Serif", serif;
background-image: url(/img/org-mode-unicorn-logo-worg.png);
background-attachment: fixed;
background-position: right bottom;
background-repeat: no-repeat;
background-color: white;
}
body {
font-size: 14pt;
line-height: 22pt;
color: black;
margin-top: 0;
}
body #content {
padding-top: 2em;
margin: auto;
max-width: 70%;
background-color: white;
}
body #support {
position: fixed;
top:0;
display:block;
font-size: 12pt;
right:0pt;
text-align: right;
padding: .2em 1em;
background: #EEE;
border-radius: 10px;
}
body .title {
margin-left: 0px;
font-size: 22pt;
}
#org-div-home-and-up{
position: fixed;
right: 0.5em;
margin-top: 70px;
font-family:sans-serif;
}
/* TOC inspired by http://jashkenas.github.com/coffee-script */
#table-of-contents {
margin-top: 105px;
font-size: 10pt;
font-family:sans-serif;
position: fixed;
right: 0em;
top: 0em;
background: white;
line-height: 12pt;
text-align: right;
box-shadow: 0 0 1em #777777;
-webkit-box-shadow: 0 0 1em #777777;
-moz-box-shadow: 0 0 1em #777777;
-webkit-border-bottom-left-radius: 5px;
-moz-border-radius-bottomleft: 5px;
/* ensure doesn't flow off the screen when expanded */
max-height: 80%;
overflow: auto; }
#table-of-contents h2 {
font-size: 13pt;
max-width: 9em;
border: 0;
font-weight: normal;
padding-left: 0.5em;
padding-right: 0.5em;
padding-top: 0.05em;
padding-bottom: 0.05em; }
#table-of-contents #text-table-of-contents {
display: none;
text-align: left; }
#table-of-contents:hover #text-table-of-contents {
display: block;
padding: 0.5em;
margin-top: -1.5em; }
#license {
background-color: #eeeeee;
}
h1 {
font-size:2.1em;
padding:0 0 30px 0;
margin-top: 10px;
margin-bottom: 10px;
margin-right: 7%;
color: grey;
}
h2 {
font-family:sans-serif;
font-size:1.45em;
padding:10px 0 10px 0;
color: black;
border-bottom: 1px solid #ddd;
padding-top: 1.5em;
}
.outline-text-2 {
margin-left: 0.1em
}
h3 {
font-family:sans-serif;
font-size:1.3em;
color: grey;
margin-left: 0.6em;
padding-top: 1.5em;
}
/* #A34D32;*/
.outline-text-3 {
margin-left: 0.9em;
}
h4 {
font-family:sans-serif;
font-size:1.2em;
margin-left: 1.2em;
color: #A5573E;
padding-top: 1.5em;
}
.outline-text-4 {
margin-left: 1.45em;
}
a {text-decoration: none; font-weight: 400;}
a:visited {text-decoration: none; font-weight: 400;}
a:hover {text-decoration: underline;}
.todo {
color: #CA0000;
}
.done {
color: #006666;
}
.timestamp-kwd {
color: #444;
}
.tag {
}
li {
margin: .4em;
}
table {
border: 0;
}
thead {
border: 0;
}
tbody {
border: 0;
}
tr {
border: 0;
}
td {
border-left: 0px;
border-right: 0px;
border-top: 0px;
border-bottom: 0px;
}
th {
border-left: 0px;
border-right: 0px;
border-top: 1px solid grey;
border-bottom: 1px solid grey;
}
code {
font-size: 100%;
color: black;
padding: 0px 0.2em;
}
img {
border: 0;
}
.share img {
opacity: .4;
-moz-opacity: .4;
filter: alpha(opacity=40);
}
.share img:hover {
opacity: 1;
-moz-opacity: 1;
filter: alpha(opacity=100);
}
pre {
font-family: Droid Sans Mono, Monaco, Consolas, "Lucida Console", monospace;
color: black;
font-size: 90%;
padding: 0.5em;
overflow: auto;
border: none;
background-color: #f2f2f2;
border-radius: 5px;
}
.org-info-box {
clear:both;
margin-left:auto;
margin-right:auto;
padding:0.7em;
}
.org-info-box img {
float:left;
margin:0em 0.5em 0em 0em;
}
.org-info-box p {
margin:0em;
padding:0em;
}
.builtin {
/* font-lock-builtin-face */
color: #f4a460;
}
.comment {
/* font-lock-comment-face */
color: #737373;
}
.comment-delimiter {
/* font-lock-comment-delimiter-face */
color: #666666;
}
.constant {
/* font-lock-constant-face */
color: #db7093;
}
.doc {
/* font-lock-doc-face */
color: #b3b3b3;
}
.function-name {
/* font-lock-function-name-face */
color: #5f9ea0;
}
.headline {
/* headline-face */
color: #ffffff;
background-color: #000000;
font-weight: bold;
}
.keyword {
/* font-lock-keyword-face */
color: #4682b4;
}
.negation-char {
}
.regexp-grouping-backslash {
}
.regexp-grouping-construct {
}
.string {
/* font-lock-string-face */
color: #ccc79a;
}
.todo-comment {
/* todo-comment-face */
color: #ffffff;
background-color: #000000;
font-weight: bold;
}
.variable-name {
/* font-lock-variable-name-face */
color: #ff6a6a;
}
.warning {
/* font-lock-warning-face */
color: #ffffff;
background-color: #cd5c5c;
font-weight: bold;
}
.important {
/* font-lock-warning-face */
background-color: #e3e3f7;
}
.note {
/* font-lock-warning-face */
background-color: #f7f7d9;
}
pre.a {
color: inherit;
background-color: inherit;
font: inherit;
text-decoration: inherit;
}
pre.a:hover {
text-decoration: underline;
}
/* Styles for org-info.js */
.org-info-js_info-navigation
{
border-style:none;
}
#org-info-js_console-label
{
font-size:10px;
font-weight:bold;
white-space:nowrap;
}
.org-info-js_search-highlight
{
background-color:#ffff00;
color:#000000;
font-weight:bold;
}
#org-info-js-window
{
border-bottom:1px solid black;
padding-bottom:10px;
margin-bottom:10px;
}
.org-info-search-highlight
{
background-color:#adefef; /* same color as emacs default */
color:#000000;
font-weight:bold;
}
.org-bbdb-company {
/* bbdb-company */
font-style: italic;
}
.org-bbdb-field-name {
}
.org-bbdb-field-value {
}
.org-bbdb-name {
/* bbdb-name */
text-decoration: underline;
}
.org-bold {
/* bold */
font-weight: bold;
}
.org-bold-italic {
/* bold-italic */
font-weight: bold;
font-style: italic;
}
.org-border {
/* border */
background-color: #000000;
}
.org-buffer-menu-buffer {
/* buffer-menu-buffer */
font-weight: bold;
}
.org-builtin {
/* font-lock-builtin-face */
color: #da70d6;
}
.org-button {
/* button */
text-decoration: underline;
}
.org-c-nonbreakable-space {
/* c-nonbreakable-space-face */
background-color: #ff0000;
font-weight: bold;
}
.org-calendar-today {
/* calendar-today */
text-decoration: underline;
}
.org-comment {
/* font-lock-comment-face */
color: #b22222;
}
.org-comment-delimiter {
/* font-lock-comment-delimiter-face */
color: #b22222;
}
.org-constant {
/* font-lock-constant-face */
color: #5f9ea0;
}
.org-cursor {
/* cursor */
background-color: #000000;
}
.org-default {
/* default */
color: #000000;
background-color: #ffffff;
}
.org-diary {
/* diary */
color: #ff0000;
}
.org-doc {
/* font-lock-doc-face */
color: #bc8f8f;
}
.org-escape-glyph {
/* escape-glyph */
color: #a52a2a;
}
.org-file-name-shadow {
/* file-name-shadow */
color: #7f7f7f;
}
.org-fixed-pitch {
}
.org-fringe {
/* fringe */
background-color: #f2f2f2;
}
.org-function-name {
/* font-lock-function-name-face */
color: #0000ff;
}
.org-header-line {
/* header-line */
color: #333333;
background-color: #e5e5e5;
}
.org-help-argument-name {
/* help-argument-name */
font-style: italic;
}
.org-highlight {
/* highlight */
background-color: #b4eeb4;
}
.org-holiday {
/* holiday */
background-color: #ffc0cb;
}
.org-info-header-node {
/* info-header-node */
color: #a52a2a;
font-weight: bold;
font-style: italic;
}
.org-info-header-xref {
/* info-header-xref */
color: #0000ff;
text-decoration: underline;
}
.org-info-menu-header {
/* info-menu-header */
font-weight: bold;
}
.org-info-menu-star {
/* info-menu-star */
color: #ff0000;
}
.org-info-node {
/* info-node */
color: #a52a2a;
font-weight: bold;
font-style: italic;
}
.org-info-title-1 {
/* info-title-1 */
font-size: 172%;
font-weight: bold;
}
.org-info-title-2 {
/* info-title-2 */
font-size: 144%;
font-weight: bold;
}
.org-info-title-3 {
/* info-title-3 */
font-size: 120%;
font-weight: bold;
}
.org-info-title-4 {
/* info-title-4 */
font-weight: bold;
}
.org-info-xref {
/* info-xref */
color: #0000ff;
text-decoration: underline;
}
.org-isearch {
/* isearch */
color: #b0e2ff;
background-color: #cd00cd;
}
.org-italic {
/* italic */
font-style: italic;
}
.org-keyword {
/* font-lock-keyword-face */
color: #a020f0;
}
.org-lazy-highlight {
/* lazy-highlight */
background-color: #afeeee;
}
.org-link {
/* link */
color: #0000ff;
text-decoration: underline;
}
.org-link-visited {
/* link-visited */
color: #8b008b;
text-decoration: underline;
}
.org-match {
/* match */
background-color: #ffff00;
}
.org-menu {
}
.org-message-cited-text {
/* message-cited-text */
color: #ff0000;
}
.org-message-header-cc {
/* message-header-cc */
color: #191970;
}
.org-message-header-name {
/* message-header-name */
color: #6495ed;
}
.org-message-header-newsgroups {
/* message-header-newsgroups */
color: #00008b;
font-weight: bold;
font-style: italic;
}
.org-message-header-other {
/* message-header-other */
color: #4682b4;
}
.org-message-header-subject {
/* message-header-subject */
color: #000080;
font-weight: bold;
}
.org-message-header-to {
/* message-header-to */
color: #191970;
font-weight: bold;
}
.org-message-header-xheader {
/* message-header-xheader */
color: #0000ff;
}
.org-message-mml {
/* message-mml */
color: #228b22;
}
.org-message-separator {
/* message-separator */
color: #a52a2a;
}
.org-minibuffer-prompt {
/* minibuffer-prompt */
color: #0000cd;
}
.org-mm-uu-extract {
/* mm-uu-extract */
color: #006400;
background-color: #ffffe0;
}
.org-mode-line {
/* mode-line */
color: #000000;
background-color: #bfbfbf;
}
.org-mode-line-buffer-id {
/* mode-line-buffer-id */
font-weight: bold;
}
.org-mode-line-highlight {
}
.org-mode-line-inactive {
/* mode-line-inactive */
color: #333333;
background-color: #e5e5e5;
}
.org-mouse {
/* mouse */
background-color: #000000;
}
.org-negation-char {
}
.org-next-error {
/* next-error */
background-color: #eedc82;
}
.org-nobreak-space {
/* nobreak-space */
color: #a52a2a;
text-decoration: underline;
}
.org-org-agenda-date {
/* org-agenda-date */
color: #0000ff;
}
.org-org-agenda-date-weekend {
/* org-agenda-date-weekend */
color: #0000ff;
font-weight: bold;
}
.org-org-agenda-restriction-lock {
/* org-agenda-restriction-lock */
background-color: #ffff00;
}
.org-org-agenda-structure {
/* org-agenda-structure */
color: #0000ff;
}
.org-org-archived {
/* org-archived */
color: #7f7f7f;
}
.org-org-code {
/* org-code */
color: #7f7f7f;
}
.org-org-column {
/* org-column */
background-color: #e5e5e5;
}
.org-org-column-title {
/* org-column-title */
background-color: #e5e5e5;
font-weight: bold;
text-decoration: underline;
}
.org-org-date {
/* org-date */
color: #a020f0;
text-decoration: underline;
}
.org-org-done {
/* org-done */
color: #228b22;
font-weight: bold;
}
.org-org-drawer {
/* org-drawer */
color: #0000ff;
}
.org-org-ellipsis {
/* org-ellipsis */
color: #b8860b;
text-decoration: underline;
}
.org-org-formula {
/* org-formula */
color: #b22222;
}
.org-org-headline-done {
/* org-headline-done */
color: #bc8f8f;
}
.org-org-hide {
/* org-hide */
color: #e5e5e5;
}
.org-org-latex-and-export-specials {
/* org-latex-and-export-specials */
color: #8b4513;
}
.org-org-level-1 {
/* org-level-1 */
color: #0000ff;
}
.org-org-level-2 {
/* org-level-2 */
color: #b8860b;
}
.org-org-level-3 {
/* org-level-3 */
color: #a020f0;
}
.org-org-level-4 {
/* org-level-4 */
color: #b22222;
}
.org-org-level-5 {
/* org-level-5 */
color: #228b22;
}
.org-org-level-6 {
/* org-level-6 */
color: #5f9ea0;
}
.org-org-level-7 {
/* org-level-7 */
color: #da70d6;
}
.org-org-level-8 {
/* org-level-8 */
color: #bc8f8f;
}
.org-org-link {
/* org-link */
color: #a020f0;
text-decoration: underline;
}
.org-org-property-value {
}
.org-org-scheduled-previously {
/* org-scheduled-previously */
color: #b22222;
}
.org-org-scheduled-today {
/* org-scheduled-today */
color: #006400;
}
.org-org-sexp-date {
/* org-sexp-date */
color: #a020f0;
}
.org-org-special-keyword {
/* org-special-keyword */
color: #bc8f8f;
}
.org-org-table {
/* org-table */
color: #0000ff;
}
.org-org-tag {
/* org-tag */
font-weight: bold;
}
.org-org-target {
/* org-target */
text-decoration: underline;
}
.org-org-time-grid {
/* org-time-grid */
color: #b8860b;
}
.org-org-todo {
/* org-todo */
color: #ff0000;
}
.org-org-upcoming-deadline {
/* org-upcoming-deadline */
color: #b22222;
}
.org-org-verbatim {
/* org-verbatim */
color: #7f7f7f;
text-decoration: underline;
}
.org-org-warning {
/* org-warning */
color: #ff0000;
font-weight: bold;
}
.org-outline-1 {
/* outline-1 */
color: #0000ff;
}
.org-outline-2 {
/* outline-2 */
color: #b8860b;
}
.org-outline-3 {
/* outline-3 */
color: #a020f0;
}
.org-outline-4 {
/* outline-4 */
color: #b22222;
}
.org-outline-5 {
/* outline-5 */
color: #228b22;
}
.org-outline-6 {
/* outline-6 */
color: #5f9ea0;
}
.org-outline-7 {
/* outline-7 */
color: #da70d6;
}
.org-outline-8 {
/* outline-8 */
color: #bc8f8f;
}
.outline-text-1, .outline-text-2, .outline-text-3, .outline-text-4, .outline-text-5, .outline-text-6 {
/* Add more spacing between section. Padding, so that folding with org-info.js works as expected. */
}
.org-preprocessor {
/* font-lock-preprocessor-face */
color: #da70d6;
}
.org-query-replace {
/* query-replace */
color: #b0e2ff;
background-color: #cd00cd;
}
.org-regexp-grouping-backslash {
/* font-lock-regexp-grouping-backslash */
font-weight: bold;
}
.org-regexp-grouping-construct {
/* font-lock-regexp-grouping-construct */
font-weight: bold;
}
.org-region {
/* region */
background-color: #eedc82;
}
.org-rmail-highlight {
}
.org-scroll-bar {
/* scroll-bar */
background-color: #bfbfbf;
}
.org-secondary-selection {
/* secondary-selection */
background-color: #ffff00;
}
.org-shadow {
/* shadow */
color: #7f7f7f;
}
.org-show-paren-match {
/* show-paren-match */
background-color: #40e0d0;
}
.org-show-paren-mismatch {
/* show-paren-mismatch */
color: #ffffff;
background-color: #a020f0;
}
.org-string {
/* font-lock-string-face */
color: #bc8f8f;
}
.org-texinfo-heading {
/* texinfo-heading */
color: #0000ff;
}
.org-tool-bar {
/* tool-bar */
color: #000000;
background-color: #bfbfbf;
}
.org-tooltip {
/* tooltip */
color: #000000;
background-color: #ffffe0;
}
.org-trailing-whitespace {
/* trailing-whitespace */
background-color: #ff0000;
}
.org-type {
/* font-lock-type-face */
color: #228b22;
}
.org-underline {
/* underline */
text-decoration: underline;
}
.org-variable-name {
/* font-lock-variable-name-face */
color: #b8860b;
}
.org-variable-pitch {
}
.org-vertical-border {
}
.org-warning {
/* font-lock-warning-face */
color: #ff0000;
font-weight: bold;
}
.rss_box {}
.rss_title, rss_title a {}
.rss_items {}
.rss_item a:link, .rss_item a:visited, .rss_item a:active {}
.rss_item a:hover {}
.rss_date {}
label.org-src-name {
font-size: 80%;
font-style: italic;
}
#show_source {margin: 0; padding: 0;}
#postamble {
font-size: 75%;
min-width: 700px;
max-width: 80%;
line-height: 14pt;
margin-left: 20px;
margin-top: 10px;
padding: .2em;
background-color: #ffffff;
z-index: -1000;
}
} /* END OF @media all */
@media screen
{
#table-of-contents {
position: fixed;
margin-top: 105px;
float: right;
border: 1px solid #red;
max-width: 50%;
overflow: auto;
}
} /* END OF @media screen */