Eorkflow

.github/workflows/gh-pages.yml

@@ -0,0 +1,34 @@
+name: github pages
+
+on:
+  push:
+    branches:
+      - master
+
+jobs:
+  deploy:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: install extra repository
+        run: sudo add-apt-repository ppa:kelleyk/emacs
+
+      - name: refresh apt
+        run: sudo apt-get update
+
+      - name: install dependencies
+        run: sudo apt-get install emacs26
+
+      - name: install htmlize
+        run: git clone https://github.com/hniksic/emacs-htmlize && cp emacs-htmlize/htmlize.el docs/
+
+      - name: make
+        run: cd docs ; ./create.sh ../QMC.org
+
+      - name: Deploy
+        uses: peaceiris/actions-gh-pages@v3
+        with:
+          github_token: ${{ secrets.GITHUB_TOKEN }}
+          publish_dir: ./docs
+

QMC.org

@@ -1253,7 +1253,7 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
     that the acceptance rate is around 0.5 for a good efficiency of
     the simulation.
 
-**** TODO Exercise 
+**** Exercise 
 
      #+begin_exercise
      Modify the previous program to introduce the accept/reject step.
@@ -1303,11 +1303,12 @@ tau = 1.0
 X = [MonteCarlo(a,tau,nmax) for i in range(30)]
 E, deltaE = ave_error([x[0] for x in X])
 A, deltaA = ave_error([x[1] for x in X])
-print(f"E = {E} +/- {deltaE}  {A} +/- {deltaA}")
+print(f"E = {E} +/- {deltaE}\nA = {A} +/- {deltaA}")
       #+END_SRC
 
       #+RESULTS:
-      : E = -0.49387078389332206 +/- 0.0033326460286729792  0.4983000000000001 +/- 0.006825097363627021
+      : E = -0.4949730317138491 +/- 0.00012478601801760644
+      : A = 0.7887163333333334 +/- 0.00026834549840347617
    
  *Fortran*
       #+BEGIN_SRC f90
@@ -1399,17 +1400,168 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
    :header-args:python: :tangle dmc.py
    :header-args:f90: :tangle dmc.f90
    :END:
+   
+** Hydrogen atom
+   
+**** Exercise 
 
-  We will now consider the H_2 molecule in a minimal basis composed of the
-  $1s$ orbitals of the hydrogen atoms:
+     #+begin_exercise
+     Modify the Metropolis VMC program to introduce the PDMC weight.
+     In the limit $\tau \rightarrow 0$, you should recover the exact
+     energy of H for any value of $a$.
+     #+end_exercise
+   
+ *Python*
+      #+BEGIN_SRC python :results output
+from hydrogen  import *
+from qmc_stats import *
 
-  $$
-  \Psi(\mathbf{r}_1, \mathbf{r}_2) =
-  \exp(-(\mathbf{r}_1 - \mathbf{R}_A)) + 
-  $$
-  where $\mathbf{r}_1$ and $\mathbf{r}_2$ denote the electron
-  coordinates and $\mathbf{R}_A$ and $\mathbf{R}_B$ the coordinates of
-  the nuclei.
+def MonteCarlo(a,tau,nmax,Eref):
+    E = 0.
+    N = 0.
+    accep_rate = 0.
+    sq_tau = np.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)
+    w = 1.0
+    for istep in range(nmax):
+        chi = np.random.normal(loc=0., scale=1.0, size=(3))
+        el = e_loc(a,r_old)
+        w *= np.exp(-tau*(el - Eref))
+        N += w
+        E += w * el
+
+        r_new = r_old + tau * d_old + sq_tau * chi
+        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
+        # PDMC weight
+        if np.random.uniform() < q:
+            accep_rate += w
+            r_old = r_new
+            d_old = d_new
+            d2_old = d2_new
+            psi_old = psi_new
+    return E/N, accep_rate/N
+
+
+a = 0.9
+nmax = 10000
+tau = .1
+X = [MonteCarlo(a,tau,nmax,-0.5) for i in range(30)]
+E, deltaE = ave_error([x[0] for x in X])
+A, deltaA = ave_error([x[1] for x in X])
+print(f"E = {E} +/- {deltaE}\nA = {A} +/- {deltaA}")
+      #+END_SRC
+
+      #+RESULTS:
+      : E = -0.49654807434947584 +/- 0.0006868522447409156
+      : A = 0.9876193891840709 +/- 0.00041857361650995804
+   
+ *Fortran*
+      #+BEGIN_SRC f90
+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 :: norm, sq_tau, chi(3), d2_old, prod, u
+  double precision :: psi_old, psi_new, d2_new, argexpo, q
+  double precision :: r_old(3), r_new(3)
+  double precision :: d_old(3), d_new(3)
+  double precision, external :: e_loc, psi
+
+  sq_tau = dsqrt(tau)
+  
+  ! Initialization
+  energy = 0.d0
+  norm   = 0.d0
+  accep_rate = 0.d0
+  call random_gauss(r_old,3)
+  call drift(a,r_old,d_old)
+  d2_old = d_old(1)*d_old(1) + d_old(2)*d_old(2) + d_old(3)*d_old(3)
+  psi_old = psi(a,r_old)
+
+  do istep = 1,nmax
+     call random_gauss(chi,3)
+     r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau
+     call drift(a,r_new,d_new)
+     d2_new = d_new(1)*d_new(1) + d_new(2)*d_new(2) + d_new(3)*d_new(3)
+     psi_new = psi(a,r_new)
+     ! Metropolis
+     prod = (d_new(1) + d_old(1))*(r_new(1) - r_old(1)) + &
+            (d_new(2) + d_old(2))*(r_new(2) - r_old(2)) + &
+            (d_new(3) + d_old(3))*(r_new(3) - r_old(3))
+     argexpo = 0.5d0 * (d2_new - d2_old)*tau + prod
+     q = psi_new / psi_old
+     q = dexp(-argexpo) * q*q
+     call random_number(u)
+     if (u<q) then
+        accep_rate = accep_rate + 1.d0
+        r_old(:) = r_new(:)
+        d_old(:) = d_new(:)
+        d2_old = d2_new
+        psi_old = psi_new
+     end if
+     norm = norm + 1.d0
+     energy = energy + e_loc(a,r_old)
+  end do
+  energy = energy / norm
+  accep_rate = accep_rate / norm
+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
+
+      #+RESULTS:
+      :  E =  -0.49499990423528023      +/-   1.5958250761863871E-004
+      :  A =   0.78861366666666655      +/-   3.5096729498002445E-004
+     
+
+** Dihydrogen
+
+   We will now consider the H_2 molecule in a minimal basis composed of the
+   $1s$ orbitals of the hydrogen atoms:
+
+   $$
+   \Psi(\mathbf{r}_1, \mathbf{r}_2) =
+   \exp(-(\mathbf{r}_1 - \mathbf{R}_A)) + 
+   $$
+   where $\mathbf{r}_1$ and $\mathbf{r}_2$ denote the electron
+   coordinates and $\mathbf{R}_A$ and $\mathbf{R}_B$ the coordinates of
+   the nuclei.
 
 
   

docs/config.el

@@ -0,0 +1,71 @@
+;; Thanks to Tobias's answer on Emacs Stack Exchange:
+;; https://emacs.stackexchange.com/questions/38437/org-mode-batch-export-missing-syntax-highlighting
+
+(package-initialize)
+(require 'htmlize)
+(require 'font-lock)
+(require 'subr-x) ;; for `when-let'
+
+(unless (boundp 'maximal-integer)
+  (defconst maximal-integer (lsh -1 -1)
+    "Maximal integer value representable natively in emacs lisp."))
+
+(defun face-spec-default (spec)
+  "Get list containing at most the default entry of face SPEC.
+Return nil if SPEC has no default entry."
+  (let* ((first (car-safe spec))
+     (display (car-safe first)))
+    (when (eq display 'default)
+      (list (car-safe spec)))))
+
+(defun face-spec-min-color (display-atts)
+  "Get min-color entry of DISPLAY-ATTS pair from face spec."
+  (let* ((display (car-safe display-atts)))
+    (or (car-safe (cdr (assoc 'min-colors display)))
+    maximal-integer)))
+
+(defun face-spec-highest-color (spec)
+  "Search face SPEC for highest color.
+That means the DISPLAY entry of SPEC
+with class 'color and highest min-color value."
+  (let ((color-list (cl-remove-if-not
+             (lambda (display-atts)
+               (when-let ((display (car-safe display-atts))
+                  (class (and (listp display)
+                          (assoc 'class display)))
+                  (background (assoc 'background display)))
+             (and (member 'light (cdr background))
+                  (member 'color (cdr class)))))
+             spec)))
+    (cl-reduce (lambda (display-atts1 display-atts2)
+         (if (> (face-spec-min-color display-atts1)
+            (face-spec-min-color display-atts2))
+             display-atts1
+           display-atts2))
+           (cdr color-list)
+           :initial-value (car color-list))))
+
+(defun face-spec-t (spec)
+  "Search face SPEC for fall back."
+  (cl-find-if (lambda (display-atts)
+        (eq (car-safe display-atts) t))
+          spec))
+
+(defun my-face-attribute (face attribute &optional frame inherit)
+  "Get FACE ATTRIBUTE from `face-user-default-spec' and not from `face-attribute'."
+  (let* ((face-spec (face-user-default-spec face))
+     (display-attr (or (face-spec-highest-color face-spec)
+               (face-spec-t face-spec)))
+     (attr (cdr display-attr))
+     (val (or (plist-get attr attribute) (car-safe (cdr (assoc attribute attr))))))
+    ;; (message "attribute: %S" attribute) ;; for debugging
+    (when (and (null (eq attribute :inherit))
+           (null val))
+      (let ((inherited-face (my-face-attribute face :inherit)))
+    (when (and inherited-face
+           (null (eq inherited-face 'unspecified)))
+      (setq val (my-face-attribute inherited-face attribute)))))
+    ;; (message "face: %S attribute: %S display-attr: %S, val: %S" face attribute display-attr val) ;; for debugging
+    (or val 'unspecified)))
+
+(advice-add 'face-attribute :override #'my-face-attribute)

docs/create.sh

@@ -0,0 +1,14 @@
+#!/bin/bash                                                                                                                   
+INPUT=$1
+
+if [[ -f ../docs/htmlize.el ]]
+then
+    emacs --batch --load ../docs/htmlize.el --load ../docs/config.el $INPUT -f org-html-export-to-html
+else
+    emacs --batch --load ../docs/config.el $INPUT -f org-html-export-to-html
+fi
+
+mv ../QMC.html index.html
+
+
+