From 9c937da3086151a5604445754bec684ad12a8703 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 2 Feb 2021 13:57:53 +0100 Subject: [PATCH] Changing a=0.9 to a=1.2 to avoid population explosion in DMC --- QMC.org | 115 ++++++++++++++++++++++++------------------- energy_hydrogen.py | 2 + hydrogen.f90 | 3 +- hydrogen.py | 1 + plot_hydrogen.py | 2 + qmc_gaussian.f90 | 2 +- qmc_gaussian.py | 4 +- qmc_metropolis.f90 | 4 +- qmc_metropolis.py | 6 ++- qmc_stats.py | 2 + qmc_uniform.f90 | 2 +- qmc_uniform.py | 4 +- variance_hydrogen.py | 2 + vmc.f90 | 4 +- vmc.py | 4 +- vmc_metropolis.f90 | 4 +- vmc_metropolis.py | 4 +- 17 files changed, 99 insertions(+), 66 deletions(-) diff --git a/QMC.org b/QMC.org index a6afabb..4fc1bb7 100644 --- a/QMC.org +++ b/QMC.org @@ -1185,7 +1185,7 @@ end subroutine ave_error *** Exercise #+begin_exercise - Parameterize the wave function with $a=0.9$. Perform 30 + Parameterize the wave function with $a=1.2$. Perform 30 independent Monte Carlo runs, each with 100 000 Monte Carlo steps. Store the final energies of each run and use this array to compute the average energy and the associated error bar. @@ -1206,7 +1206,7 @@ from qmc_stats import * def MonteCarlo(a, nmax): # TODO -a = 0.9 +a = 1.2 nmax = 100000 #TODO @@ -1244,7 +1244,7 @@ end subroutine uniform_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -1287,7 +1287,7 @@ def MonteCarlo(a, nmax): return energy / normalization -a = 0.9 +a = 1.2 nmax = 100000 X = [MonteCarlo(a,nmax) for i in range(30)] @@ -1297,7 +1297,7 @@ print(f"E = {E} +/- {deltaE}") #+END_SRC #+RESULTS: - : E = -0.4956255109300764 +/- 0.0007082875482711226 + : E = -0.4773221805255284 +/- 0.0022489426510479975 *Fortran* #+begin_note @@ -1341,7 +1341,7 @@ end subroutine uniform_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -1365,7 +1365,7 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform #+end_src #+RESULTS: - : E = -0.49518773675598715 +/- 5.2391494923686175E-004 + : E = -0.48084122147238995 +/- 2.4983775878329355E-003 ** Metropolis sampling with $\Psi^2$ :PROPERTIES: @@ -1420,10 +1420,14 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform which, for our choice of transition probability, becomes $$ - A(\mathbf{r}_{n}\rightarrow\mathbf{r}_{n+1}) = \min\left(1,\frac{P(\mathbf{r}_{n+1})}{P(\mathbf{r}_{n})}\right)= \min\left(1,\frac{\Psi(\mathbf{r}_{n+1})^2}{\Psi(\mathbf{r}_{n})^2}\right)\,. + A(\mathbf{r}_{n}\rightarrow\mathbf{r}_{n+1}) = \min\left(1,\frac{P(\mathbf{r}_{n+1})}{P(\mathbf{r}_{n})}\right)= \min\left(1,\frac{|\Psi(\mathbf{r}_{n+1})|^2}{|\Psi(\mathbf{r}_{n})|^2}\right)\,. $$ - Explain why the transition probability cancels out in the expression of $A$. Also note that we do not need to compute the norm of the wave function! + #+begin_exercise + Explain why the transition probability cancels out in the + expression of $A$. + #+end_exercise + Also note that we do not need to compute the norm of the wave function! The algorithm is summarized as follows: @@ -1439,27 +1443,32 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform A common error is to remove the rejected samples from the calculation of the average. *Don't do it!* - All samples should be kept, from both accepted and rejected moves. + All samples should be kept, from both accepted /and/ rejected moves. #+end_note - If the box is infinitely small, the ratio will be very close - to one and all the steps will be accepted. However, the moves will be - very correlated and you will visit the configurational space very slowly. - - On the other hand, if you propose too large moves, the number of - accepted steps will decrease because the ratios might become - small. If the number of accepted steps is close to zero, then the - space is not well sampled either. - - The size of the move should be adjusted so that it is as large as - possible, keeping the number of accepted steps not too small. To - achieve that, we define the acceptance rate as the number of - accepted steps over the total number of steps. Adjusting the time - step such that the acceptance rate is close to 0.5 is a good - compromise for the current problem. - NOTE: below, we use the symbol $\delta t$ to denote $\delta L$ since we will use - the same variable later on to store a time step. +*** Optimal step size + + If the box is infinitely small, the ratio will be very close + to one and all the steps will be accepted. However, the moves will be + very correlated and you will visit the configurational space very slowly. + + On the other hand, if you propose too large moves, the number of + accepted steps will decrease because the ratios might become + small. If the number of accepted steps is close to zero, then the + space is not well sampled either. + + The size of the move should be adjusted so that it is as large as + possible, keeping the number of accepted steps not too small. To + achieve that, we define the acceptance rate as the number of + accepted steps over the total number of steps. Adjusting the time + step such that the acceptance rate is close to 0.5 is a good + compromise for the current problem. + + #+begin_note + Below, we use the symbol $\delta t$ to denote $\delta L$ since we will use + the same variable later on to store a time step. + #+end_note *** Exercise @@ -1489,7 +1498,7 @@ def MonteCarlo(a,nmax,dt): # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 dt = #TODO @@ -1506,6 +1515,8 @@ A, deltaA = ave_error(X) print(f"A = {A} +/- {deltaA}") #+END_SRC + #+RESULTS: + *Fortran* #+BEGIN_SRC f90 :tangle none subroutine metropolis_montecarlo(a,nmax,dt,energy,accep) @@ -1529,7 +1540,7 @@ end subroutine metropolis_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9d0 + double precision, parameter :: a = 1.2d0 double precision, parameter :: dt = ! TODO integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -1588,9 +1599,9 @@ def MonteCarlo(a,nmax,dt): return energy/nmax, N_accep/nmax # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 -dt = 1.3 +dt = 1.0 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)] @@ -1606,8 +1617,8 @@ print(f"A = {A} +/- {deltaA}") #+END_SRC #+RESULTS: - : E = -0.4950720838131573 +/- 0.00019089638602238043 - : A = 0.5172960000000001 +/- 0.0003443446549306529 + : E = -0.4802595860693983 +/- 0.0005124420418289207 + : A = 0.5074913333333334 +/- 0.000350889422714878 *Fortran* #+BEGIN_SRC f90 :exports code @@ -1662,8 +1673,8 @@ end subroutine metropolis_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9d0 - double precision, parameter :: dt = 1.3d0 + double precision, parameter :: a = 1.2d0 + double precision, parameter :: dt = 1.0d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -1689,8 +1700,8 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_metropolis.f90 -o qmc_metropolis ./qmc_metropolis #+end_src #+RESULTS: - : E = -0.49503130891988767 +/- 1.7160104275040037E-004 - : A = 0.51695266666666673 +/- 4.0445505648997396E-004 + : E = -0.47948142754168033 +/- 4.8410177863919307E-004 + : A = 0.50762633333333318 +/- 3.4601756760043725E-004 ** Gaussian random number generator @@ -1899,7 +1910,7 @@ def MonteCarlo(a,nmax,dt): # TODO # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 dt = # TODO @@ -1939,7 +1950,7 @@ end subroutine variational_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 double precision, parameter :: dt = ! TODO integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -2014,7 +2025,7 @@ def MonteCarlo(a,nmax,dt): # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 dt = 1.3 @@ -2112,8 +2123,8 @@ end subroutine variational_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 - double precision, parameter :: dt = 1.0 + double precision, parameter :: a = 1.2d0 + double precision, parameter :: dt = 1.0d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -2442,7 +2453,7 @@ def MonteCarlo(a, nmax, dt, Eref): # TODO # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 dt = 0.01 E_ref = -0.5 @@ -2484,7 +2495,7 @@ end subroutine pdmc program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 double precision, parameter :: dt = 0.1d0 double precision, parameter :: E_ref = -0.5d0 double precision, parameter :: tau = 10.d0 @@ -2576,7 +2587,7 @@ def MonteCarlo(a, nmax, dt, tau, Eref): # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 dt = 0.1 tau = 10. @@ -2694,7 +2705,7 @@ end subroutine pdmc program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 double precision, parameter :: dt = 0.1d0 double precision, parameter :: E_ref = -0.5d0 double precision, parameter :: tau = 10.d0 @@ -2814,7 +2825,7 @@ def MonteCarlo(a,nmax): E += w * e_loc(a,r) return E/N -a = 0.9 +a = 1.2 nmax = 100000 X = [MonteCarlo(a,nmax) for i in range(30)] E, deltaE = ave_error(X) @@ -2860,7 +2871,7 @@ end subroutine gaussian_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -3003,7 +3014,7 @@ def MonteCarlo(a,dt,nmax): return E/N -a = 0.9 +a = 1.2 nmax = 100000 dt = 0.2 X = [MonteCarlo(a,dt,nmax) for i in range(30)] @@ -3046,8 +3057,8 @@ end subroutine variational_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 - double precision, parameter :: dt = 0.2 + double precision, parameter :: a = 1.2d0 + double precision, parameter :: dt = 0.2d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 @@ -3093,4 +3104,6 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc | <2021-02-04 Thu 12:15-12:30> | 2.4 | |------------------------------+---------| | <2021-02-04 Thu 14:00-14:10> | 3.1 | + | <2021-02-04 Thu 14:10-14:30> | 3.2 | + | <2021-02-04 Thu 14:30-15:30> | 3.3 | |------------------------------+---------| diff --git a/energy_hydrogen.py b/energy_hydrogen.py index 14be354..9168a71 100644 --- a/energy_hydrogen.py +++ b/energy_hydrogen.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + import numpy as np from hydrogen import e_loc, psi diff --git a/hydrogen.f90 b/hydrogen.f90 index 338f42b..0eedc53 100644 --- a/hydrogen.f90 +++ b/hydrogen.f90 @@ -43,7 +43,8 @@ double precision function e_loc(a,r) implicit none double precision, intent(in) :: a, r(3) - double precision, external :: kinetic, potential + double precision, external :: kinetic + double precision, external :: potential e_loc = kinetic(a,r) + potential(r) diff --git a/hydrogen.py b/hydrogen.py index b3746c5..ebbd960 100644 --- a/hydrogen.py +++ b/hydrogen.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 import numpy as np def potential(r): diff --git a/plot_hydrogen.py b/plot_hydrogen.py index 25701f4..33a76f0 100644 --- a/plot_hydrogen.py +++ b/plot_hydrogen.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + import numpy as np import matplotlib.pyplot as plt diff --git a/qmc_gaussian.f90 b/qmc_gaussian.f90 index 6adfc3f..9a4b32b 100644 --- a/qmc_gaussian.f90 +++ b/qmc_gaussian.f90 @@ -32,7 +32,7 @@ end subroutine gaussian_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 diff --git a/qmc_gaussian.py b/qmc_gaussian.py index 129c79f..1afdbaa 100644 --- a/qmc_gaussian.py +++ b/qmc_gaussian.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + from hydrogen import * from qmc_stats import * @@ -16,7 +18,7 @@ def MonteCarlo(a,nmax): E += w * e_loc(a,r) return E/N -a = 0.9 +a = 1.2 nmax = 100000 X = [MonteCarlo(a,nmax) for i in range(30)] E, deltaE = ave_error(X) diff --git a/qmc_metropolis.f90 b/qmc_metropolis.f90 index a3db63f..0f1f3b7 100644 --- a/qmc_metropolis.f90 +++ b/qmc_metropolis.f90 @@ -49,8 +49,8 @@ end subroutine metropolis_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9d0 - double precision, parameter :: dt = 1.3d0 + double precision, parameter :: a = 1.2d0 + double precision, parameter :: dt = 1.0d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 diff --git a/qmc_metropolis.py b/qmc_metropolis.py index 8581d30..022c601 100644 --- a/qmc_metropolis.py +++ b/qmc_metropolis.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + from hydrogen import * from qmc_stats import * @@ -25,9 +27,9 @@ def MonteCarlo(a,nmax,dt): return energy/nmax, N_accep/nmax # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 -dt = 1.3 +dt = 1.0 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)] diff --git a/qmc_stats.py b/qmc_stats.py index e21ebb9..2888195 100644 --- a/qmc_stats.py +++ b/qmc_stats.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + from math import sqrt def ave_error(arr): M = len(arr) diff --git a/qmc_uniform.f90 b/qmc_uniform.f90 index ec95768..b2dec69 100644 --- a/qmc_uniform.f90 +++ b/qmc_uniform.f90 @@ -31,7 +31,7 @@ end subroutine uniform_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 + double precision, parameter :: a = 1.2d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 diff --git a/qmc_uniform.py b/qmc_uniform.py index 66a8616..ce938f0 100644 --- a/qmc_uniform.py +++ b/qmc_uniform.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + from hydrogen import * from qmc_stats import * @@ -16,7 +18,7 @@ def MonteCarlo(a, nmax): return energy / normalization -a = 0.9 +a = 1.2 nmax = 100000 X = [MonteCarlo(a,nmax) for i in range(30)] diff --git a/variance_hydrogen.py b/variance_hydrogen.py index abaf3bd..07b7d17 100644 --- a/variance_hydrogen.py +++ b/variance_hydrogen.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + import numpy as np from hydrogen import e_loc, psi diff --git a/vmc.f90 b/vmc.f90 index 6a3dde8..e3a4514 100644 --- a/vmc.f90 +++ b/vmc.f90 @@ -28,8 +28,8 @@ end subroutine variational_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 - double precision, parameter :: dt = 0.2 + double precision, parameter :: a = 1.2d0 + double precision, parameter :: dt = 0.2d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 diff --git a/vmc.py b/vmc.py index b34321e..e4fba20 100644 --- a/vmc.py +++ b/vmc.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + from hydrogen import * from qmc_stats import * @@ -19,7 +21,7 @@ def MonteCarlo(a,dt,nmax): return E/N -a = 0.9 +a = 1.2 nmax = 100000 dt = 0.2 X = [MonteCarlo(a,dt,nmax) for i in range(30)] diff --git a/vmc_metropolis.f90 b/vmc_metropolis.f90 index 499f09b..4bbb498 100644 --- a/vmc_metropolis.f90 +++ b/vmc_metropolis.f90 @@ -73,8 +73,8 @@ end subroutine variational_montecarlo program qmc implicit none - double precision, parameter :: a = 0.9 - double precision, parameter :: dt = 1.0 + double precision, parameter :: a = 1.2d0 + double precision, parameter :: dt = 1.0d0 integer*8 , parameter :: nmax = 100000 integer , parameter :: nruns = 30 diff --git a/vmc_metropolis.py b/vmc_metropolis.py index 4181fe4..93cd5b6 100644 --- a/vmc_metropolis.py +++ b/vmc_metropolis.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python3 + from hydrogen import * from qmc_stats import * @@ -41,7 +43,7 @@ def MonteCarlo(a,nmax,dt): # Run simulation -a = 0.9 +a = 1.2 nmax = 100000 dt = 1.3