LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2025-03-25 10:06:32 +01:00
Code Releases Activity
qmcchem/src/SAMPLING/fkmc_step.irp.f

376 lines
11 KiB
Fortran
Raw Permalink Blame History

! Providers of *_fkmc_block_walk
!==============================
BEGIN_SHELL [ /usr/bin/env python2 ]
from properties import *
t = """
BEGIN_PROVIDER [ $T, $X_fkmc_block_walk $D1 ]
&BEGIN_PROVIDER [ $T, $X_fkmc_block_walk_kahan $D2 ]
&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk $D1 ]
&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk_kahan $D2 ]
implicit none
BEGIN_DOC
! fkMC averages of $X. Computed in E_loc_fkmc_block_walk
END_DOC
$X_fkmc_block_walk = 0.d0
$X_fkmc_block_walk_kahan = 0.d0
$X_2_fkmc_block_walk = 0.d0
$X_2_fkmc_block_walk_kahan = 0.d0
END_PROVIDER
"""
for p in properties:
if p[1] != 'e_loc':
if p[2] == "":
D1 = ""
D2 = ", (3)"
else:
D1 = ", ("+p[2][1:-1]+")"
D2 = ", ("+p[2][1:-1]+",3)"
print t.replace("$X",p[1]).replace("$T",p[0]).replace("$D1",D1).replace("$D2",D2)
END_SHELL
BEGIN_PROVIDER [ double precision, E_loc_fkmc_block_walk ]
&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk ]
&BEGIN_PROVIDER [ double precision, E_loc_fkmc_block_walk_kahan, (3) ]
&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk_kahan, (3) ]
implicit none
include '../types.F'
BEGIN_DOC
! Properties averaged over the block using the FKMC method
END_DOC
integer, parameter :: BIRTH=1, DEATH=2
real, allocatable :: elec_coord_tmp(:,:,:)
integer :: mod_align
double precision :: E_loc_save(walk_num_dmc_max)
double precision :: E_loc_save_tmp(walk_num_dmc_max)
double precision :: psi_value_save(walk_num)
double precision :: psi_value_save_tmp(walk_num)
double precision :: fkmc_weight(walk_num)
double precision :: delta(walk_num)
double precision, allocatable :: psi_grad_psi_inv_save(:,:,:)
double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:)
double precision, allocatable :: fkmc_clock_tmp(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save_tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save_tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_weight
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_clock_tmp
allocate ( psi_grad_psi_inv_save(elec_num_8,3,walk_num), &
psi_grad_psi_inv_save_tmp(elec_num_8,3,walk_num), &
elec_coord_tmp(mod_align(elec_num+1),3,walk_num), &
fkmc_clock_tmp(2,walk_num) )
psi_value_save = 0.d0
psi_value_save_tmp = 0.d0
fkmc_weight = 1.d0
! Initialization
if (vmc_algo /= t_Brownian) then
call abrt(irp_here,'FKMC should run with Brownian algorithm')
endif
integer :: k, i_walk, i_step
BEGIN_SHELL [ /usr/bin/env python2 ]
from properties import *
t = """
if (calc_$X) then
!DIR$ VECTOR ALIGNED
$X_fkmc_block_walk = 0.d0
!DIR$ VECTOR ALIGNED
$X_fkmc_block_walk_kahan = 0.d0
!DIR$ VECTOR ALIGNED
$X_2_fkmc_block_walk = 0.d0
!DIR$ VECTOR ALIGNED
$X_2_fkmc_block_walk_kahan = 0.d0
endif
"""
for p in properties:
print t.replace("$X",p[1])
END_SHELL
logical :: loop
integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max
loop = .True.
call system_clock(cpu0, count_rate, count_max)
cpu2 = cpu0
block_weight = 0.d0
real, external :: accep_rate
double precision :: thr
thr = 2.d0/time_step_sq
logical :: first_loop
first_loop = .True.
do while (loop)
! Every walker makes a step
do i_walk=1,walk_num
if (.not.first_loop) then
integer :: i,j,l
do l=1,3
do i=1,elec_num+1
elec_coord(i,l) = elec_coord_full(i,l,i_walk)
enddo
do i=1,elec_num
psi_grad_psi_inv_x(i) = psi_grad_psi_inv_save(i,1,i_walk)
psi_grad_psi_inv_y(i) = psi_grad_psi_inv_save(i,2,i_walk)
psi_grad_psi_inv_z(i) = psi_grad_psi_inv_save(i,3,i_walk)
enddo
psi_value = psi_value_save(i_walk)
E_loc = E_loc_save(i_walk)
enddo
SOFT_TOUCH elec_coord psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z psi_value E_loc
else
do l=1,3
do i=1,elec_num+1
elec_coord(i,l) = elec_coord_full(i,l,i_walk)
enddo
enddo
TOUCH elec_coord
E_loc_save(i_walk) = E_loc
psi_value_save(i_walk) = psi_value
endif
double precision :: p,q
real :: delta_x
logical :: accepted
call brownian_step(p,q,accepted,delta_x)
if ( psi_value * psi_value_save(i_walk) >= 0.d0 ) then
delta(i_walk) = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p
if ( delta(i_walk) > thr ) then
delta(i_walk) = thr
else if ( delta(i_walk) < -thr ) then
delta(i_walk) = -thr
endif
fkmc_weight(i_walk) = dexp(-dtime_step*delta(i_walk))
elec_coord(elec_num+1,1) += p*time_step
elec_coord(elec_num+1,2) = E_loc
elec_coord(elec_num+1,3) = fkmc_weight(i_walk)
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,i_walk) = elec_coord(i,l)
enddo
enddo
do i=1,elec_num
psi_grad_psi_inv_save(i,1,i_walk) = psi_grad_psi_inv_x(i)
psi_grad_psi_inv_save(i,2,i_walk) = psi_grad_psi_inv_y(i)
psi_grad_psi_inv_save(i,3,i_walk) = psi_grad_psi_inv_z(i)
enddo
psi_value_save(i_walk) = psi_value
E_loc_save(i_walk) = E_loc
BEGIN_SHELL [ /usr/bin/env python2 ]
from properties import *
t = """
if (calc_$X) then
! Kahan's summation algorithm to compute these sums reducing the rounding error:
! $X_fkmc_block_walk += $X * fkmc_weight(i_walk)
! $X_2_fkmc_block_walk += $X_2 * fkmc_weight(i_walk)
! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
$X_fkmc_block_walk_kahan($D2 3) = $X * fkmc_weight(i_walk) - $X_fkmc_block_walk_kahan($D2 1)
$X_fkmc_block_walk_kahan($D2 2) = $X_fkmc_block_walk $D1 + $X_fkmc_block_walk_kahan($D2 3)
$X_fkmc_block_walk_kahan($D2 1) = ($X_fkmc_block_walk_kahan($D2 2) - $X_fkmc_block_walk $D1 ) &
- $X_fkmc_block_walk_kahan($D2 3)
$X_fkmc_block_walk $D1 = $X_fkmc_block_walk_kahan($D2 2)
$X_2_fkmc_block_walk_kahan($D2 3) = $X_2 * fkmc_weight(i_walk) - $X_2_fkmc_block_walk_kahan($D2 1)
$X_2_fkmc_block_walk_kahan($D2 2) = $X_2_fkmc_block_walk $D1 + $X_2_fkmc_block_walk_kahan($D2 3)
$X_2_fkmc_block_walk_kahan($D2 1) = ($X_2_fkmc_block_walk_kahan($D2 2) - $X_2_fkmc_block_walk $D1 ) &
- $X_2_fkmc_block_walk_kahan($D2 3)
$X_2_fkmc_block_walk $D1 = $X_2_fkmc_block_walk_kahan($D2 2)
endif
"""
for p in properties:
if p[2] == "":
D1 = ""
D2 = ""
else:
D1 = "("+":"*(p[2].count(',')+1)+")"
D2 = ":"*(p[2].count(',')+1)+","
print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2)
END_SHELL
block_weight += fkmc_weight(i_walk)
else
fkmc_weight(i_walk) = 0.d0
delta(i_walk) = 1.d5
endif
enddo
! Compute the new weight of the population
double precision :: sum_weight
sum_weight = 0.d0
do k=1,walk_num
sum_weight += fkmc_weight(k)
enddo
do k=1,walk_num
do l=1,3
do i=1,elec_num+1
elec_coord_tmp(i,l,k) = elec_coord_full(i,l,k)
enddo
do i=1,elec_num
psi_grad_psi_inv_save_tmp(i,l,k) = psi_grad_psi_inv_save(i,l,k)
enddo
enddo
psi_value_save_tmp(k) = psi_value_save(k)
E_loc_save_tmp(k) = E_loc_save(k)
if (fkmc_weight(k) == 0.d0) then
fkmc_clock(DEATH,k) = -1.d0
endif
if ( delta(k) <= 0.d0 ) then
fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k) +time_step * delta(k)
fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k)
else
fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k)
fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k) -time_step * delta(k)
endif
enddo
! Reconfiguration
! ===============
! Identify first which walkers will be killed to place branched walkers there
! later
double precision, external :: qmc_ranf
integer :: ipm, m
integer :: killed(walk_num)
m=1
do k=1,walk_num
fkmc_clock(DEATH,k) = fkmc_clock_tmp(DEATH,k)
if (fkmc_clock_tmp(DEATH,k) <= 0.d0) then
killed(m) = k
m += 1
fkmc_clock(DEATH,k) = -dlog(qmc_ranf())
fkmc_clock(BIRTH,k) = -dlog(qmc_ranf())
ipm = k
do while (ipm == k)
ipm = 1 + int (walk_num*qmc_ranf())
enddo
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm)
enddo
do i=1,elec_num
psi_grad_psi_inv_save(i,l,k) = psi_grad_psi_inv_save_tmp(i,l,ipm)
enddo
enddo
psi_value_save(k) = psi_value_save_tmp(ipm)
E_loc_save(k) = E_loc_save_tmp(ipm)
endif
enddo
killed(m) = 0
m=1
do k=1,walk_num
fkmc_clock(BIRTH,k) = fkmc_clock_tmp(BIRTH,k)
if (fkmc_clock_tmp(BIRTH,k) <= 0.d0) then
fkmc_clock(BIRTH,k) = -dlog(qmc_ranf())
if (killed(m) == 0) then
ipm = k
do while (ipm == k)
ipm = 1 + int (walk_num*qmc_ranf())
enddo
else
ipm = killed(m)
m +=1
endif
fkmc_clock(BIRTH,ipm) = -dlog(qmc_ranf())
fkmc_clock(DEATH,ipm) = -dlog(qmc_ranf())
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,ipm) = elec_coord_tmp(i,l,k)
enddo
do i=1,elec_num
psi_grad_psi_inv_save(i,l,ipm) = psi_grad_psi_inv_save_tmp(i,l,k)
enddo
enddo
psi_value_save(ipm) = psi_value_save_tmp(k)
E_loc_save(ipm) = E_loc_save_tmp(k)
endif
enddo
call system_clock(cpu1, count_rate, count_max)
if (cpu1 < cpu0) then
cpu1 = cpu1+cpu0
endif
loop = dble(cpu1-cpu0)/dble(count_rate) < block_time
if (cpu1-cpu2 > count_rate) then
integer :: do_run
call get_running(do_run)
loop = loop.and.(do_run == t_Running)
cpu2 = cpu1
endif
! Update E_ref to take into account the weight of the population
E_ref -= dlog(sum_weight / dble(walk_num) ) / time_step
SOFT_TOUCH elec_coord_full E_ref
first_loop = .False.
enddo
double precision :: factor
factor = 1.d0/block_weight
SOFT_TOUCH block_weight
BEGIN_SHELL [ /usr/bin/env python2 ]
from properties import *
t = """
if (calc_$X) then
$X_fkmc_block_walk *= factor
$X_2_fkmc_block_walk *= factor
endif
"""
for p in properties:
print t.replace("$X",p[1])
END_SHELL
deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp, &
fkmc_clock_tmp )
END_PROVIDER
BEGIN_PROVIDER [ double precision, fkmc_clock, (2,walk_num) ]
implicit none
BEGIN_DOC
! Branching clocks for the FKMC algotithm. (1,:) is the birth clock and
! (2,:) is the death clock.
END_DOC
integer :: i
double precision, external :: qmc_ranf
do i=1, walk_num
fkmc_clock(1,i) = -dlog(qmc_ranf())
fkmc_clock(2,i) = -dlog(qmc_ranf())
enddo
END_PROVIDER
View Git Blame Copy Permalink