Improved async DMC

This commit is contained in:
Anthony Scemama 2022-01-06 17:44:55 +01:00
parent 5b1379fd9c
commit a0323922a8
1 changed files with 77 additions and 29 deletions

View File

@ -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