2022-01-03 13:18:45 +01:00
|
|
|
#+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
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
real, allocatable :: elec_coord_new(:,:,:)
|
2022-01-03 13:18:45 +01:00
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
double precision :: w
|
|
|
|
double precision, allocatable :: E_out(:), w_sum(:)
|
2022-01-03 13:18:45 +01:00
|
|
|
|
|
|
|
double precision, external :: qmc_ranf
|
|
|
|
|
|
|
|
allocate(elec_coord_new(elec_num+1,3,walk_num))
|
2022-01-06 17:44:55 +01:00
|
|
|
allocate(E_out(walk_num), w_sum(walk_num))
|
2022-01-03 13:18:45 +01:00
|
|
|
#+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
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
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))
|
|
|
|
|
2022-01-03 13:18:45 +01:00
|
|
|
<<branching>>
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
if (k_new >= walk_num) then
|
|
|
|
w_sum(k_full+1:) = 0.d0
|
|
|
|
exit
|
|
|
|
end if
|
2022-01-03 13:18:45 +01:00
|
|
|
end do
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
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)
|
2022-01-03 13:18:45 +01:00
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
! call write_coords(k_new)
|
|
|
|
call write_energy(walk_num, E_out, w_sum)
|
2022-01-03 13:18:45 +01:00
|
|
|
end do
|
|
|
|
|
|
|
|
end subroutine run
|
|
|
|
|
|
|
|
<<read_coords>>
|
|
|
|
<<write_coords>>
|
2022-01-06 17:44:55 +01:00
|
|
|
<<write_energy>>
|
2022-01-03 13:18:45 +01:00
|
|
|
<<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) = &
|
2022-01-06 17:44:55 +01:00
|
|
|
elec_coord(1:elec_num+1,1:3)
|
2022-01-03 13:18:45 +01:00
|
|
|
k_new = k_new+1
|
2022-01-06 17:44:55 +01:00
|
|
|
if (k_new >= walk_num) exit
|
2022-01-03 13:18:45 +01:00
|
|
|
end do
|
|
|
|
|
|
|
|
#+end_src
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
* Read/write
|
2022-01-03 13:18:45 +01:00
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
** Read coordinates
|
2022-01-03 13:18:45 +01:00
|
|
|
|
|
|
|
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
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
** Write coordinates
|
2022-01-03 13:18:45 +01:00
|
|
|
|
|
|
|
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
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
** 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
|
|
|
|
|
2022-01-03 13:18:45 +01:00
|
|
|
* 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:
|
2022-01-06 17:44:55 +01:00
|
|
|
- ~w~: the last of all $w_i$
|
2022-01-03 13:18:45 +01:00
|
|
|
- ~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
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
! Max number of steps
|
|
|
|
integer :: imax
|
|
|
|
integer, parameter :: nmax=10000
|
|
|
|
|
2022-01-03 13:18:45 +01:00
|
|
|
! Brownian step variables
|
|
|
|
double precision :: p,q
|
|
|
|
real :: delta_x
|
|
|
|
logical :: accepted
|
|
|
|
|
|
|
|
! Local energies from the past
|
|
|
|
double precision :: E_loc_save(4)
|
2022-01-06 17:44:55 +01:00
|
|
|
double precision :: w
|
2022-01-03 13:18:45 +01:00
|
|
|
|
|
|
|
#+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.
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
do imax = 1, nmax
|
2022-01-03 13:18:45 +01:00
|
|
|
|
|
|
|
call brownian_step(p,q,accepted,delta_x)
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
! 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
|
2022-01-03 13:18:45 +01:00
|
|
|
delta = (delta - E_ref)*p
|
|
|
|
|
|
|
|
if (delta >= 0.d0) then
|
2022-01-06 17:44:55 +01:00
|
|
|
w = dexp(-dtime_step*delta)
|
2022-01-03 13:18:45 +01:00
|
|
|
else
|
2022-01-06 17:44:55 +01:00
|
|
|
w = 2.d0-dexp(dtime_step*delta)
|
2022-01-03 13:18:45 +01:00
|
|
|
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
|
2022-01-06 17:44:55 +01:00
|
|
|
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
|
2022-01-03 13:18:45 +01:00
|
|
|
|
|
|
|
w_sum = w_sum + pdmc_weight
|
|
|
|
E_out = E_out + pdmc_weight * E_loc
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
pdmc_weight = pdmc_weight * w
|
|
|
|
|
|
|
|
loop = pdmc_weight > 0.5d0 .and. pdmc_weight < 2.0d0
|
|
|
|
if (.not.loop) exit
|
2022-01-03 13:18:45 +01:00
|
|
|
|
|
|
|
end do
|
|
|
|
|
2022-01-06 17:44:55 +01:00
|
|
|
E_out = E_out / w_sum
|
|
|
|
|
2022-01-03 13:18:45 +01:00
|
|
|
end subroutine pdmc_trajectory
|
|
|
|
#+end_src
|
|
|
|
|