mirror of https://gitlab.com/scemama/qmcchem.git
284 lines
7.7 KiB
Fortran
284 lines
7.7 KiB
Fortran
! Providers of *_dmc_block_walk
|
|
!==============================
|
|
BEGIN_SHELL [ /usr/bin/python ]
|
|
from properties import *
|
|
|
|
t = """
|
|
BEGIN_PROVIDER [ $T, $X_dmc_block_walk $D1 ]
|
|
&BEGIN_PROVIDER [ $T, $X_dmc_block_walk_kahan $D2 ]
|
|
&BEGIN_PROVIDER [ $T, $X_2_dmc_block_walk $D1 ]
|
|
&BEGIN_PROVIDER [ $T, $X_2_dmc_block_walk_kahan $D2 ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! VMC averages of $X
|
|
END_DOC
|
|
$X_dmc_block_walk = 0.d0
|
|
$X_dmc_block_walk_kahan = 0.d0
|
|
$X_2_dmc_block_walk = 0.d0
|
|
$X_2_dmc_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_dmc_block_walk ]
|
|
&BEGIN_PROVIDER [ double precision, E_loc_2_dmc_block_walk ]
|
|
&BEGIN_PROVIDER [ double precision, E_loc_dmc_block_walk_kahan, (3) ]
|
|
&BEGIN_PROVIDER [ double precision, E_loc_2_dmc_block_walk_kahan, (3)
|
|
implicit none
|
|
include '../types.F'
|
|
BEGIN_DOC
|
|
! Properties averaged over the block using the DMC method
|
|
END_DOC
|
|
|
|
real, allocatable :: elec_coord_tmp(:,:,:)
|
|
integer :: mod_align
|
|
double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:)
|
|
double precision :: psi_value_save_tmp(walk_num)
|
|
integer :: trapped_walk_tmp(walk_num)
|
|
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save_tmp
|
|
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp
|
|
allocate ( elec_coord_tmp(mod_align(elec_num+1),3,walk_num) )
|
|
allocate ( psi_grad_psi_inv_save_tmp(elec_num_8,3,walk_num) )
|
|
|
|
! Initialization
|
|
if (vmc_algo /= t_Brownian) then
|
|
call abrt(irp_here,'DMC should run with Brownian algorithm')
|
|
endif
|
|
PROVIDE E_loc_vmc_block_walk
|
|
|
|
integer :: k, i_walk, i_step
|
|
|
|
BEGIN_SHELL [ /usr/bin/python ]
|
|
from properties import *
|
|
t = """
|
|
if (calc_$X) then
|
|
$X_dmc_block_walk = 0.d0
|
|
$X_dmc_block_walk_kahan = 0.d0
|
|
$X_2_dmc_block_walk = 0.d0
|
|
$X_2_dmc_block_walk_kahan = 0.d0
|
|
$X_min = huge(1.)
|
|
$X_max =-huge(1.)
|
|
endif
|
|
"""
|
|
for p in properties:
|
|
print t.replace("$X",p[1])
|
|
END_SHELL
|
|
|
|
double precision :: icount
|
|
|
|
icount = 0.d0
|
|
|
|
logical :: loop
|
|
integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max
|
|
|
|
loop = .True.
|
|
call system_clock(cpu0, count_rate, count_max)
|
|
cpu2 = cpu0
|
|
do while (loop)
|
|
! Move to the next projection step
|
|
dmc_projection_step = mod(dmc_projection_step,dmc_projection)+1
|
|
|
|
! Remove contribution of the old value of the weight at the new
|
|
! projection step
|
|
pop_weight_mult *= 1.d0/pop_weight(dmc_projection_step)
|
|
|
|
! Compute the new weight of the population
|
|
pop_weight(dmc_projection_step) = 0.d0
|
|
do k=1,walk_num
|
|
pop_weight(dmc_projection_step) += dmc_weight(k)
|
|
enddo
|
|
|
|
! Normalize the weight of the walkers by the weight of the population
|
|
do k=1,walk_num
|
|
dmc_weight(k) = dmc_weight(k)/pop_weight(dmc_projection_step)
|
|
enddo
|
|
|
|
! Normalize the weight of the population at the current projection step by
|
|
! the number of walkers
|
|
pop_weight(dmc_projection_step) = pop_weight(dmc_projection_step)/dble(walk_num)
|
|
|
|
! Update the running population weight
|
|
pop_weight_mult *= pop_weight(dmc_projection_step)
|
|
SOFT_TOUCH pop_weight_mult
|
|
|
|
BEGIN_SHELL [ /usr/bin/python ]
|
|
from properties import *
|
|
t = """
|
|
if (calc_$X) then
|
|
! Kahan's summation algorithm to compute these sums reducing the rounding error:
|
|
! $X_dmc_block_walk($D) += $X * pop_weight_mult
|
|
! $X_2_dmc_block_walk($D) += $X_2 * pop_weight_mult
|
|
! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
|
|
|
|
$X_dmc_block_walk_kahan($D2 3) = $X * pop_weight_mult - $X_dmc_block_walk_kahan($D2 1)
|
|
$X_dmc_block_walk_kahan($D2 2) = $X_dmc_block_walk $D1 + $X_dmc_block_walk_kahan($D2 3)
|
|
$X_dmc_block_walk_kahan($D2 1) = ($X_dmc_block_walk_kahan($D2 2) - $X_dmc_block_walk $D1 ) &
|
|
- $X_dmc_block_walk_kahan($D2 3)
|
|
$X_dmc_block_walk $D1 = $X_dmc_block_walk_kahan($D2 2)
|
|
|
|
|
|
$X_2_dmc_block_walk_kahan($D2 3) = $X_2 * pop_weight_mult - $X_2_dmc_block_walk_kahan($D2 1)
|
|
$X_2_dmc_block_walk_kahan($D2 2) = $X_2_dmc_block_walk $D1 + $X_2_dmc_block_walk_kahan($D2 3)
|
|
$X_2_dmc_block_walk_kahan($D2 1) = ($X_2_dmc_block_walk_kahan($D2 2) - $X_2_dmc_block_walk $D1 ) &
|
|
- $X_2_dmc_block_walk_kahan($D2 3)
|
|
$X_2_dmc_block_walk $D1 = $X_2_dmc_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
|
|
icount += pop_weight_mult
|
|
|
|
! Reconfiguration
|
|
|
|
integer :: ipos(walk_num)
|
|
call reconfigure(ipos,dmc_weight)
|
|
|
|
do k=1,walk_num
|
|
integer :: i, l
|
|
do l=1,3
|
|
do i=1,elec_num+1
|
|
elec_coord_tmp(i,l,k) = elec_coord_full(i,l,k)
|
|
enddo
|
|
!DIR$ VECTOR ALIGNED
|
|
!DIR$ LOOP COUNT(200)
|
|
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)
|
|
trapped_walk_tmp(k) = trapped_walk(k)
|
|
enddo
|
|
|
|
integer :: ipm
|
|
do k=1,walk_num
|
|
ipm = ipos(k)
|
|
do l=1,3
|
|
do i=1,elec_num+1
|
|
elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm)
|
|
enddo
|
|
!DIR$ VECTOR ALIGNED
|
|
!DIR$ LOOP COUNT(200)
|
|
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)
|
|
trapped_walk(k) = trapped_walk_tmp(ipm)
|
|
enddo
|
|
|
|
! Set 1st walker
|
|
!DIR$ VECTOR ALIGNED
|
|
!DIR$ LOOP COUNT(200)
|
|
do i=1,elec_num
|
|
psi_grad_psi_inv_x(i) = psi_grad_psi_inv_save(i,1,1)
|
|
psi_grad_psi_inv_y(i) = psi_grad_psi_inv_save(i,2,1)
|
|
psi_grad_psi_inv_z(i) = psi_grad_psi_inv_save(i,3,1)
|
|
enddo
|
|
|
|
!DIR$ VECTOR UNALIGNED
|
|
!DIR$ LOOP COUNT(200)
|
|
do i=1,elec_num
|
|
elec_coord(i,1) = elec_coord_full(i,1,1)
|
|
elec_coord(i,2) = elec_coord_full(i,2,1)
|
|
elec_coord(i,3) = elec_coord_full(i,3,1)
|
|
enddo
|
|
psi_value = psi_value_save(1)
|
|
|
|
TOUCH elec_coord_full psi_value_save psi_grad_psi_inv_save psi_value psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z elec_coord
|
|
|
|
call system_clock(cpu1, count_rate, count_max)
|
|
if (cpu1 < cpu0) then
|
|
cpu1 = cpu1+cpu0
|
|
endif
|
|
loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate)
|
|
if (cpu1-cpu2 > count_rate) then
|
|
integer :: do_run
|
|
call get_running(do_run)
|
|
loop = do_run == t_Running
|
|
cpu2 = cpu1
|
|
endif
|
|
|
|
enddo
|
|
|
|
|
|
|
|
double precision :: factor
|
|
factor = 1.d0/icount
|
|
block_weight *= icount
|
|
SOFT_TOUCH block_weight
|
|
BEGIN_SHELL [ /usr/bin/python ]
|
|
from properties import *
|
|
t = """
|
|
if (calc_$X) then
|
|
$X_dmc_block_walk *= factor
|
|
$X_2_dmc_block_walk *= factor
|
|
endif
|
|
"""
|
|
for p in properties:
|
|
print t.replace("$X",p[1])
|
|
END_SHELL
|
|
|
|
deallocate ( elec_coord_tmp )
|
|
deallocate ( psi_grad_psi_inv_save_tmp )
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, E_ref ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Weight of the DMC population
|
|
END_DOC
|
|
E_ref = 0.d0
|
|
call get_simulation_E_ref(E_ref)
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, pop_weight_mult ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Population weight of DMC
|
|
END_DOC
|
|
pop_weight_mult = pop_weight(dmc_projection)
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ integer, dmc_projection ]
|
|
&BEGIN_PROVIDER [ integer, dmc_projection_step ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Number of projection steps for SRMC
|
|
END_DOC
|
|
real :: dmc_projection_time
|
|
dmc_projection_time = 1.
|
|
call get_simulation_dmc_projection_time(dmc_projection_time)
|
|
dmc_projection = int( dmc_projection_time/time_step)
|
|
dmc_projection_step = 0
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, pop_weight, (dmc_projection) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Population weight of DMC
|
|
END_DOC
|
|
pop_weight = 1.d0
|
|
pop_weight(dmc_projection) = 1.d0/dble(dmc_projection)
|
|
END_PROVIDER
|
|
|