TREX/qmc-lttc
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmc-lttc.git synced 2024-12-21 11:53:58 +01:00
Code Issues Projects Releases Wiki Activity

Introduced fixed-node error

Browse Source
This commit is contained in:
Anthony Scemama 2021-01-27 15:21:44 +01:00
parent d49102cf2e
commit f9354e4401
1 changed files with 67 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

163
QMC.org
View File

@ -714,98 +714,83 @@ 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 $$ $$\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 #+end_exercise
**** TODO Solution :solution: **** 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 *** Exercise
#+begin_exercise #+begin_exercise
Add the calculation of the variance to the previous code, and Add the calculation of the variance to the previous code, and
compute a numerical estimate of the variance of the local energy compute a numerical estimate of the variance of the local energy in
in a grid of $50\times50\times50$ points in the range a grid of $50\times50\times50$ points in the range $(-5,-5,-5) \le
$(-5,-5,-5) \mathbf{r} \le (5,5,5)$ for different values of $a$.
\le \mathbf{r} \le (5,5,5)$ for different values of $a$.
#+end_exercise #+end_exercise
*Python* *Python*
#+begin_src python :results none :tangle none #+begin_src python :results none :tangle none
import numpy as np import numpy as np from hydrogen import e_loc, psi
from hydrogen import e_loc, psi
interval = np.linspace(-5,5,num=50) interval = np.linspace(-5,5,num=50) delta =
delta = (interval[1]-interval[0])**3 (interval[1]-interval[0])**3
r = np.array([0.,0.,0.]) r = np.array([0.,0.,0.])
for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
# TODO # TODO
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
*Fortran* *Fortran*
#+begin_src f90 :tangle none #+begin_src f90 :tangle none
program variance_hydrogen program variance_hydrogen implicit none double precision, external ::
implicit none e_loc, psi double precision :: x(50), w, delta, energy, dx, r(3),
double precision, external :: e_loc, psi a(6), norm, s2 double precision :: e, energy2 integer :: i, k, l, j
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 /) a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
dx = 10.d0/(size(x)-1) dx = 10.d0/(size(x)-1) do i=1,size(x) x(i) = -5.d0 + (i-1)*dx end do
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
delta = dx**3 delta = dx**3
r(:) = 0.d0 r(:) = 0.d0
do j=1,size(a) do j=1,size(a) ! TODO print *, 'a = ', a(j), ' E = ', energy, ' s2 =
! TODO ', s2 end do
print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2
end do
end program variance_hydrogen end program variance_hydrogen #+end_src
#+end_src
To compile and run: To compile and run:
#+begin_src sh :results output :exports both #+begin_src sh :results output :exports both
gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
./variance_hydrogen ./variance_hydrogen #+end_src
#+end_src
**** Solution :solution: **** Solution :solution:
*Python* *Python*
#+begin_src python :results none :exports both #+begin_src python :results none :exports both
import numpy as np import numpy as np from hydrogen import e_loc, psi
from hydrogen import e_loc, psi
interval = np.linspace(-5,5,num=50) interval = np.linspace(-5,5,num=50) delta =
delta = (interval[1]-interval[0])**3 (interval[1]-interval[0])**3
r = np.array([0.,0.,0.]) r = np.array([0.,0.,0.])
for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: E = 0. E2 = 0. norm = 0.
E = 0. for x in interval: r[0] = x for y in interval: r[1] = y for z in
E2 = 0. interval: r[2] = z w = psi(a, r) w = w * w * delta El = e_loc(a,
norm = 0. r) E += w * El E2 += w * El*El norm += w E = E / norm E2 = E2 /
for x in interval: norm s2 = E2 - E*E print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 =
r[0] = x {s2:10.8f}") #+end_src
for y in interval:
r[1] = y
for z in interval:
r[2] = z
w = psi(a, r)
w = w * w * delta
El = e_loc(a, r)
E += w * El
E2 += w * El*El
norm += w
E = E / norm
E2 = E2 / norm
s2 = E2 - E*E
print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")
#+end_src
#+RESULTS: #+RESULTS:
: a = 0.1 E = -0.24518439 \sigma^2 = 0.02696522 : a = 0.1 E = -0.24518439 \sigma^2 = 0.02696522
@ -818,56 +803,31 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
*Fortran* *Fortran*
#+begin_src f90 #+begin_src f90
program variance_hydrogen program variance_hydrogen implicit none double precision, external ::
implicit none e_loc, psi double precision :: x(50), w, delta, energy, dx, r(3),
double precision, external :: e_loc, psi a(6), norm, s2 double precision :: e, energy2 integer :: i, k, l, j
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 /) a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
dx = 10.d0/(size(x)-1) dx = 10.d0/(size(x)-1) do i=1,size(x) x(i) = -5.d0 + (i-1)*dx end do
do i=1,size(x)
x(i) = -5.d0 + (i-1)*dx
end do
delta = dx**3 delta = dx**3
r(:) = 0.d0 r(:) = 0.d0
do j=1,size(a) do j=1,size(a) energy = 0.d0 energy2 = 0.d0 norm = 0.d0 do
energy = 0.d0 i=1,size(x) r(1) = x(i) do k=1,size(x) r(2) = x(k) do l=1,size(x)
energy2 = 0.d0 r(3) = x(l) w = psi(a(j),r) w = w * w * delta e = e_loc(a(j), r)
norm = 0.d0 energy = energy + w * e energy2 = energy2 + w * e * e norm =
do i=1,size(x) norm + w end do end do end do energy = energy / norm energy2 =
r(1) = x(i) energy2 / norm s2 = energy2 - energy*energy print *, 'a = ',
do k=1,size(x) a(j), ' E = ', energy, ' s2 = ', s2 end do
r(2) = x(k)
do l=1,size(x)
r(3) = x(l)
w = psi(a(j),r)
w = w * w * delta
e = e_loc(a(j), r)
energy = energy + w * e
energy2 = energy2 + w * e * e
norm = norm + w
end do
end do
end do
energy = energy / norm
energy2 = energy2 / norm
s2 = energy2 - energy*energy
print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2
end do
end program variance_hydrogen end program variance_hydrogen #+end_src
#+end_src
#+begin_src sh :results output :exports results #+begin_src sh :results output :exports results
gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
./variance_hydrogen ./variance_hydrogen #+end_src
#+end_src
#+RESULTS: #+RESULTS:
: a = 0.10000000000000001 E = -0.24518438948809140 s2 = 2.6965218719722767E-002 : a = 0.10000000000000001 E = -0.24518438948809140 s2 = 2.6965218719722767E-002
@ -1860,7 +1820,6 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
imaginary time will converge to the exact ground state of the imaginary time will converge to the exact ground state of the
system. system.
** Diffusion and branching ** Diffusion and branching
The [[https://en.wikipedia.org/wiki/Diffusion_equation][diffusion equation]] of particles is given by The [[https://en.wikipedia.org/wiki/Diffusion_equation][diffusion equation]] of particles is given by
@ -1904,7 +1863,8 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
practice. practice.
Fortunately, if we multiply the Schrödinger equation by a chosen Fortunately, if we multiply the Schrödinger equation by a chosen
/trial wave function/ $\Psi_T(\mathbf{r})$ (Hartree-Fock, Kohn-Sham /trial wave function/ $\Psi_T(\mathbf{r})$ (Hartree-Fock, Kohn-Sham
determinant, CI wave function, /etc/), one obtains determinant, CI wave function, /etc/), one obtains (see appendix
for details)
\[ \[
-\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) = -\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) =
@ -1926,7 +1886,16 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
can be simulated by the drifted diffusion scheme presented in the can be simulated by the drifted diffusion scheme presented in the
previous section (VMC). previous section (VMC).
This equation generates the /N/-electron density $\Pi$, which is the
product of the ground state with the trial wave function. It
introduces the constraint that $\Pi(\mathbf{r},\tau)=0$ where
$\Psi_T(\mathbf{r})=0$. If the wave function has the same sign
everywhere, as in the hydrogen atom or the H_2 molecule, this
constraint is harmless and we can obtain the exact energy. But for
systems where the wave function has nodes (systems with 3 or more
electrons), this scheme introduces an error known as the /fixed node
error/.
*** Appendix : Details of the Derivation *** Appendix : Details of the Derivation
\[ \[
@ -1973,6 +1942,8 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
\right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau) \right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau)
\] \]
** Pure Diffusion Monte Carlo
** TODO Hydrogen atom ** TODO Hydrogen atom
*** Exercise *** Exercise