From a0323922a8f3457881a3c1da50c02fe62a006bc5 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 6 Jan 2022 17:44:55 +0100
Subject: [PATCH] Improved async DMC

---
 src/MAIN/admc.org | 106 +++++++++++++++++++++++++++++++++-------------
 1 file changed, 77 insertions(+), 29 deletions(-)

diff --git a/src/MAIN/admc.org b/src/MAIN/admc.org
index b4b1912..ba4432c 100644
--- a/src/MAIN/admc.org
+++ b/src/MAIN/admc.org
@@ -20,13 +20,15 @@
   integer          :: iw      ! Number of copies in branching
   integer          :: l
 
-  double precision :: w, E_out, w_sum
+  real, allocatable :: elec_coord_new(:,:,:)
 
-  double precision, allocatable :: elec_coord_new(:,:,:)
+  double precision :: w
+  double precision, allocatable :: E_out(:), w_sum(:)
 
   double precision, external :: qmc_ranf                                                  
 
   allocate(elec_coord_new(elec_num+1,3,walk_num))
+  allocate(E_out(walk_num), w_sum(walk_num))
   #+end_src
 
 ** Main flow
@@ -54,30 +56,33 @@ subroutine run
      call abrt(irp_here,'DMC should run with Brownian algorithm')
   endif
 
-  do iter=1,2
-     call read_coords()
-     
-     k_new = 1
-     do while (k_new < walk_num)
-        call pdmc_trajectory(k_full, w, E_out, w_sum)
-        write(*,*) 'E', E_out, w_sum
-        k_full = k_full+1
-        if (k_full > walk_num) k_full = 1
-        
+  do iter=1,1000
+!    call read_coords()
+     k_new  = 1
+     do k_full=1,walk_num
+        call pdmc_trajectory(k_full, w, E_out(k_full), w_sum(k_full))
+
         <<branching>>
 
+       if (k_new >= walk_num) then
+         w_sum(k_full+1:) = 0.d0
+         exit
+       end if
      end do
   
-     elec_coord_full(1:elec_num+1,1:3,1:walk_num) = &
-          elec_coord_new(1:elec_num+1,1:3,1:walk_num)
+     k_new = k_new-1
+     elec_coord_full(1:elec_num+1,1:3,1:k_new) = &
+          elec_coord_new(1:elec_num+1,1:3,1:k_new)
 
-     call write_coords()
+!     call write_coords(k_new)
+     call write_energy(walk_num, E_out, w_sum)
   end do
 
 end subroutine run
 
 <<read_coords>>
 <<write_coords>>
+<<write_energy>>
 <<pdmc_trajectory>>
   #+end_src
 
@@ -95,16 +100,16 @@ end subroutine run
   ! Duplicate walker
   do l=1,iw
     elec_coord_new(1:elec_num+1,1:3,k_new) = &
-         elec_coord_full(1:elec_num+1,1:3,k_full) 
+         elec_coord(1:elec_num+1,1:3) 
     k_new = k_new+1
-    if (k_new > walk_num) exit
+    if (k_new >= walk_num) exit
   end do
 
     #+end_src
 
-* Read/write coordinates
+* Read/write
   
-** Read
+** Read coordinates
 
    Fetch a new set of coordinates for ~walk_num~ walkers from the pool of coordinates.
   
@@ -126,7 +131,7 @@ subroutine read_coords()
 end subroutine read_coords
    #+end_src
 
-** Write
+** Write coordinates
 
    Send the current set of coordinates for ~walk_num~ walkers to the pool of coordinates.
   
@@ -146,6 +151,35 @@ subroutine write_coords()
 end subroutine write_coords
    #+end_src
 
+** Write energy
+
+   Compute the weighted average over the computed energies.
+  \[
+     E = \frac{\sum_i w_i E_i}{\sum_i w_i}
+  \]
+  
+   #+NAME: write_energy
+   #+begin_src f90
+subroutine write_energy(walk_num_, E_out, w_sum)
+  implicit none
+  integer, intent(in) :: walk_num_
+  double precision, intent(in) :: E_out(walk_num_)
+  double precision, intent(in) :: w_sum(walk_num_)
+
+  integer :: i, k
+  double precision :: E, S
+
+  E = 0.d0
+  S = 0.d0
+  do k=1,walk_num
+     S = S + w_sum(k)
+     E = E + w_sum(k) * E_out(k)
+  end do
+  write(*,*) 'E', E/S, S
+
+end subroutine write_energy
+   #+end_src
+
 * PDMC trajectory
   
   Computes a PDMC trajectory until the weight ~w~ is $1/2 < w < 3/2$.
@@ -155,7 +189,7 @@ end subroutine write_coords
   \]
 
   The function returns:
-  - ~pdmc_weight~: the last of all $w_i$
+  - ~w~: the last of all $w_i$
   - ~E_out~: The average energy $E$ of the trajectory
   - ~w_sum~: The sum of the weights
     
@@ -167,6 +201,10 @@ end subroutine write_coords
   ! If true, continue to make more steps
   logical :: loop
 
+  ! Max number of steps
+  integer :: imax
+  integer, parameter :: nmax=10000
+
   ! Brownian step variables
   double precision               :: p,q
   real                           :: delta_x
@@ -174,6 +212,7 @@ end subroutine write_coords
 
   ! Local energies from the past
   double precision :: E_loc_save(4)
+  double precision :: w
 
   #+end_src
   
@@ -198,32 +237,41 @@ subroutine pdmc_trajectory(k_full, pdmc_weight, E_out, w_sum)
   pdmc_weight = 1.d0
   loop = .True.
 
-  do while (loop)
+  do imax = 1, nmax
 
      call brownian_step(p,q,accepted,delta_x)
 
-     delta = (9.d0*E_loc+19.d0*E_loc_save(1)-5.d0*E_loc_save(2)+E_loc_save(3))/24.d0
+!    delta = (9.d0*E_loc+19.d0*E_loc_save(1)-5.d0*E_loc_save(2)+E_loc_save(3))/24.d0
+     delta = E_loc
      delta = (delta - E_ref)*p
 
      if (delta >= 0.d0) then
-       pdmc_weight = dexp(-dtime_step*delta)
+       w = dexp(-dtime_step*delta)
      else
-       pdmc_weight = 2.d0-dexp(dtime_step*delta)
+       w = 2.d0-dexp(dtime_step*delta)
      endif
      elec_coord(elec_num+1,1) += p*time_step
      elec_coord(elec_num+1,2)  = E_loc
      elec_coord(elec_num+1,3)  = pdmc_weight
+     if (accepted) then
+        E_loc_save(4) = E_loc_save(3) 
+        E_loc_save(3) = E_loc_save(2) 
+        E_loc_save(2) = E_loc_save(1) 
+        E_loc_save(1) = E_loc
+     endif
 
      w_sum = w_sum + pdmc_weight
      E_out = E_out + pdmc_weight * E_loc
 
-     loop = pdmc_weight > 0.5d0 .or. pdmc_weight < 1.5d0
+     pdmc_weight = pdmc_weight * w
+
+     loop = pdmc_weight > 0.5d0 .and. pdmc_weight < 2.0d0
+     if (.not.loop) exit
      
   end do
 
-  elec_coord_full(1:elec_num+1,1:3,k_full) = elec_coord(1:elec_num+1,1:3)
-  SOFT_TOUCH elec_coord_full
-  
+  E_out = E_out / w_sum
+
 end subroutine pdmc_trajectory
    #+end_src