#+TITLE: Asynchronous DMC
#+AUTHOR: Anthony Scemama
#+EMAIL: scemama@irsamc.ups-tlse.fr


#+PROPERTY: header-args  :tangle no :noweb yes 


* Main program

** Declarations

  #+NAME: declarations
  #+begin_src f90 
  include '../types.F'

  integer          :: iter
  integer          :: k_full  ! Index of walkers in elec_coord_full
  integer          :: k_new   ! Index of walkers in elec_coord_new
  integer          :: iw      ! Number of copies in branching
  integer          :: l

  double precision :: w, E_out, w_sum

  double precision, allocatable :: elec_coord_new(:,:,:)

  double precision, external :: qmc_ranf                                                  

  allocate(elec_coord_new(elec_num+1,3,walk_num))
  #+end_src

** Main flow

   - Fetch ~walk_num~ electron coordinates in ~elec_coord_full~
   - For each set of coordinates,
     - Make a PDMC trajectory, and output the weight ~w~
     - Perform branching depending on the value of the weight
     - Store the new sets of coordinates in ~elec_coord_new~
   - When ~elec_coord_new~ is full, send it to the server
     
  #+begin_src f90 :tangle "admc.irp.f" 
program admc
  call run
  call ezfio_finish
end program admc

subroutine run
  implicit none

  <<declarations>>
  
  ! Initialization
  if (vmc_algo /= t_Brownian) then
     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
        
        <<branching>>

     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)

     call write_coords()
  end do

end subroutine run

<<read_coords>>
<<write_coords>>
<<pdmc_trajectory>>
  #+end_src

** Branching

   #+NAME: branching
    #+begin_src f90
  ! Find number of copies
  iw = int(w)
  w = w - int(w)
  if (qmc_ranf() < w) then
     iw = iw+1
  end if

  ! 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) 
    k_new = k_new+1
    if (k_new > walk_num) exit
  end do

    #+end_src

* Read/write coordinates
  
** Read

   Fetch a new set of coordinates for ~walk_num~ walkers from the pool of coordinates.
  
   #+NAME: read_coords
   #+begin_src f90
subroutine read_coords()
  implicit none

  integer :: i, k

  do k=1,walk_num
     do i=1,elec_num
        read(*,*) elec_coord_full(i,1:3,k)
     end do
  end do

  SOFT_TOUCH elec_coord_full
     
end subroutine read_coords
   #+end_src

** Write

   Send the current set of coordinates for ~walk_num~ walkers to the pool of coordinates.
  
   #+NAME: write_coords
   #+begin_src f90
subroutine write_coords()
  implicit none

  integer :: i, k

  do k=1,walk_num
     do i=1,elec_num
        write(*,*) 'C', elec_coord_full(i,1:3,k)
     end do
  end do

end subroutine write_coords
   #+end_src

* PDMC trajectory
  
  Computes a PDMC trajectory until the weight ~w~ is $1/2 < w < 3/2$.
  The energy of the trajectory is computed as
  \[
     E = \frac{\sum_i w_i E(R_i)}{\sum_i w_i}
  \]

  The function returns:
  - ~pdmc_weight~: the last of all $w_i$
  - ~E_out~: The average energy $E$ of the trajectory
  - ~w_sum~: The sum of the weights
    
  #+NAME: declarations_pdmc
  #+begin_src f90
  integer          :: i,j,l
  double precision :: delta

  ! If true, continue to make more steps
  logical :: loop

  ! Brownian step variables
  double precision               :: p,q
  real                           :: delta_x
  logical                        :: accepted

  ! Local energies from the past
  double precision :: E_loc_save(4)

  #+end_src
  
   #+NAME: pdmc_trajectory
   #+begin_src f90
subroutine pdmc_trajectory(k_full, pdmc_weight, E_out, w_sum)
  implicit none

  integer, intent(in) :: k_full
  double precision, intent(out) :: pdmc_weight, E_out, w_sum

  <<declarations_pdmc>>
  
  elec_coord(1:elec_num+1,1:3) = elec_coord_full(1:elec_num+1,1:3,k_full)
  TOUCH elec_coord

  E_out = 0.d0
  w_sum = 0.d0

  E_loc_save(1:4) = E_loc

  pdmc_weight = 1.d0
  loop = .True.

  do while (loop)

     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 = (delta - E_ref)*p

     if (delta >= 0.d0) then
       pdmc_weight = dexp(-dtime_step*delta)
     else
       pdmc_weight = 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

     w_sum = w_sum + pdmc_weight
     E_out = E_out + pdmc_weight * E_loc

     loop = pdmc_weight > 0.5d0 .or. pdmc_weight < 1.5d0
     
  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
  
end subroutine pdmc_trajectory
   #+end_src