mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-06-02 11:25:18 +02:00
367 lines
11 KiB
Fortran
367 lines
11 KiB
Fortran
! Providers of *_pdmc_block_walk
|
|
!==============================
|
|
BEGIN_SHELL [ /usr/bin/python ]
|
|
from properties import *
|
|
|
|
t = """
|
|
BEGIN_PROVIDER [ $T, $X_pdmc_block_walk $D1 ]
|
|
&BEGIN_PROVIDER [ $T, $X_pdmc_block_walk_kahan $D2 ]
|
|
&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk $D1 ]
|
|
&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk_kahan $D2 ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! pdMC averages of $X. Computed in E_loc_pdmc_block_walk
|
|
END_DOC
|
|
$X_pdmc_block_walk = 0.d0
|
|
$X_pdmc_block_walk_kahan = 0.d0
|
|
$X_2_pdmc_block_walk = 0.d0
|
|
$X_2_pdmc_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_pdmc_block_walk ]
|
|
&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk ]
|
|
&BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk_kahan, (3) ]
|
|
&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk_kahan, (3) ]
|
|
implicit none
|
|
include '../types.F'
|
|
BEGIN_DOC
|
|
! Properties averaged over the block using the PDMC method
|
|
END_DOC
|
|
|
|
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 :: pdmc_weight(walk_num)
|
|
double precision, allocatable :: psi_grad_psi_inv_save(:,:,:)
|
|
double precision, allocatable :: psi_grad_psi_inv_save_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 :: pdmc_weight
|
|
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) )
|
|
psi_value_save = 0.d0
|
|
psi_value_save_tmp = 0.d0
|
|
pdmc_weight = 1.d0
|
|
|
|
! Initialization
|
|
if (vmc_algo /= t_Brownian) then
|
|
call abrt(irp_here,'PDMC should run with Brownian algorithm')
|
|
endif
|
|
|
|
integer :: k, i_walk, i_step
|
|
|
|
BEGIN_SHELL [ /usr/bin/python ]
|
|
from properties import *
|
|
t = """
|
|
if (calc_$X) then
|
|
!DIR$ VECTOR ALIGNED
|
|
$X_pdmc_block_walk = 0.d0
|
|
!DIR$ VECTOR ALIGNED
|
|
$X_pdmc_block_walk_kahan = 0.d0
|
|
!DIR$ VECTOR ALIGNED
|
|
$X_2_pdmc_block_walk = 0.d0
|
|
!DIR$ VECTOR ALIGNED
|
|
$X_2_pdmc_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 :: delta, thr
|
|
|
|
thr = 2.d0/time_step_sq
|
|
|
|
logical :: first_loop
|
|
first_loop = .True.
|
|
if (walk_num > 1) then
|
|
call abrt(irp_here,'walk_num > 1')
|
|
endif
|
|
|
|
integer :: info
|
|
double precision :: H(pdmc_n_diag,pdmc_n_diag), S(pdmc_n_diag,pdmc_n_diag), w(pdmc_n_diag), work(3*pdmc_n_diag)
|
|
H = 0.d0
|
|
S = 0.d0
|
|
|
|
do while (loop)
|
|
|
|
i_walk = 1
|
|
|
|
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
|
|
psi_value_save(i_walk) = psi_value
|
|
E_loc_save(i_walk) = E_loc
|
|
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 = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p
|
|
if ( delta > thr ) then
|
|
pdmc_weight(i_walk) = dexp(-dtime_step*thr)
|
|
else if ( delta < -thr ) then
|
|
pdmc_weight(i_walk) = dexp(dtime_step*thr)
|
|
else
|
|
pdmc_weight(i_walk) = 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(i_walk) * pdmc_pop_weight_mult(1)
|
|
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
|
|
|
|
if (dabs(pdmc_weight(i_walk)*pdmc_pop_weight_mult(1)) > 1.d-15) then
|
|
dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(1))
|
|
dmc_zv_weight_half = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(2))
|
|
else
|
|
dmc_zv_weight = 0.d0
|
|
dmc_zv_weight_half = 0.d0
|
|
endif
|
|
TOUCH dmc_zv_weight dmc_zv_weight_half
|
|
|
|
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_pdmc_block_walk += $X * pdmc_pop_weight_mult(1) * pdmc_weight(i_walk)
|
|
! $X_2_pdmc_block_walk += $X_2 * pdmc_pop_weight_mult(1) * pdmc_weight(i_walk)
|
|
! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
|
|
|
|
$X_pdmc_block_walk_kahan($D2 3) = $X * pdmc_pop_weight_mult(1) * pdmc_weight(i_walk) - $X_pdmc_block_walk_kahan($D2 1)
|
|
$X_pdmc_block_walk_kahan($D2 2) = $X_pdmc_block_walk $D1 + $X_pdmc_block_walk_kahan($D2 3)
|
|
$X_pdmc_block_walk_kahan($D2 1) = ($X_pdmc_block_walk_kahan($D2 2) - $X_pdmc_block_walk $D1 ) &
|
|
- $X_pdmc_block_walk_kahan($D2 3)
|
|
$X_pdmc_block_walk $D1 = $X_pdmc_block_walk_kahan($D2 2)
|
|
|
|
|
|
$X_2_pdmc_block_walk_kahan($D2 3) = $X_2 * pdmc_pop_weight_mult(1) * pdmc_weight(i_walk) - $X_2_pdmc_block_walk_kahan($D2 1)
|
|
$X_2_pdmc_block_walk_kahan($D2 2) = $X_2_pdmc_block_walk $D1 + $X_2_pdmc_block_walk_kahan($D2 3)
|
|
$X_2_pdmc_block_walk_kahan($D2 1) = ($X_2_pdmc_block_walk_kahan($D2 2) - $X_2_pdmc_block_walk $D1 ) &
|
|
- $X_2_pdmc_block_walk_kahan($D2 3)
|
|
$X_2_pdmc_block_walk $D1 = $X_2_pdmc_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 += pdmc_pop_weight_mult(1) * pdmc_weight(i_walk)
|
|
|
|
! H(1,1) += E_loc + (E_trial-E_loc)
|
|
! H(2,2) += E_loc* pdmc_weight(i_walk) + (E_trial-E_loc)
|
|
! H(3,3) += E_loc* pdmc_pop_weight_mult(1) * pdmc_weight(i_walk) + (E_trial-E_loc)
|
|
! H(1,2) += E_loc * pdmc_pop_weight_mult(2) * pdmc_weight(i_walk) + (E_trial-E_loc)
|
|
|
|
H(1,1) += E_loc
|
|
H(2,2) += E_loc*pdmc_pop_weight_mult(1) * pdmc_weight(i_walk)
|
|
H(1,2) += E_loc*pdmc_pop_weight_mult(2) * pdmc_weight(i_walk)
|
|
H(2,1) = H(1,2)
|
|
|
|
H = H+E_trial-E_loc
|
|
|
|
S(1,1) += 1.d0
|
|
S(2,2) += pdmc_pop_weight_mult(1) * pdmc_weight(i_walk)
|
|
S(1,2) += pdmc_pop_weight_mult(2) * pdmc_weight(i_walk)
|
|
S(2,1) = S(1,2)
|
|
|
|
else
|
|
pdmc_weight(i_walk) = 1.d0
|
|
pdmc_pop_weight(:,:) = 1.d0
|
|
pdmc_pop_weight_mult(:) = 1.d0
|
|
endif
|
|
|
|
|
|
do k=1,pdmc_n_diag
|
|
! Move to the next projection step
|
|
if (pdmc_projection > 0) then
|
|
pdmc_projection_step(k) = mod(pdmc_projection_step(k),pdmc_projection/k)+1
|
|
else
|
|
pdmc_projection_step(k) = 1
|
|
endif
|
|
|
|
! Eventually, recompute the weight of the population
|
|
if (pdmc_projection_step(k) == k) then
|
|
pdmc_pop_weight_mult(k) = 1.d0
|
|
do l=1,pdmc_projection/k
|
|
pdmc_pop_weight_mult(k) *= pdmc_pop_weight(l,k)
|
|
enddo
|
|
endif
|
|
! Remove contribution of the old value of the weight at the new
|
|
! projection step
|
|
|
|
pdmc_pop_weight_mult(k) *= 1.d0/pdmc_pop_weight(pdmc_projection_step(k),k)
|
|
pdmc_pop_weight(pdmc_projection_step(k),k) = pdmc_weight(i_walk)/dble(walk_num)
|
|
|
|
! Update the running population weight
|
|
pdmc_pop_weight_mult(k) *= pdmc_pop_weight(pdmc_projection_step(k),k)
|
|
|
|
enddo
|
|
|
|
|
|
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
|
|
|
|
SOFT_TOUCH elec_coord_full pdmc_pop_weight_mult
|
|
|
|
first_loop = .False.
|
|
|
|
enddo
|
|
|
|
|
|
double precision :: factor
|
|
factor = 1.d0/block_weight
|
|
SOFT_TOUCH block_weight
|
|
|
|
BEGIN_SHELL [ /usr/bin/python ]
|
|
from properties import *
|
|
t = """
|
|
if (calc_$X) then
|
|
$X_pdmc_block_walk *= factor
|
|
$X_2_pdmc_block_walk *= factor
|
|
endif
|
|
"""
|
|
for p in properties:
|
|
print t.replace("$X",p[1])
|
|
END_SHELL
|
|
|
|
|
|
call dsygv(1, 'N', 'U', pdmc_n_diag, H, pdmc_n_diag, S, pdmc_n_diag, w, work, 3*pdmc_n_diag, info)
|
|
E_loc_zv_diag_pdmc_block_walk = w(1)
|
|
print *, w
|
|
! E_loc_zv_diag_pdmc_block_walk = 0.5d0 * (sqrt( &
|
|
! E_loc_zv_pdmc_block_walk(2)*E_loc_zv_pdmc_block_walk(2) &
|
|
! - 2.d0*E_trial*E_loc_zv_pdmc_block_walk(2) &
|
|
! + 5.d0*E_trial*E_trial &
|
|
! - 8.d0*E_loc_zv_pdmc_block_walk(1)*E_trial &
|
|
! + 4.d0*E_loc_zv_pdmc_block_walk(1)*E_loc_zv_pdmc_block_walk(1) &
|
|
! ) - E_loc_zv_pdmc_block_walk(2) - E_trial)
|
|
|
|
deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp )
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ integer, pdmc_projection ]
|
|
&BEGIN_PROVIDER [ integer, pdmc_projection_step, (pdmc_n_diag) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Number of projection steps for PDMC
|
|
END_DOC
|
|
real :: pdmc_projection_time
|
|
pdmc_projection_time = 1.
|
|
call get_simulation_srmc_projection_time(pdmc_projection_time)
|
|
pdmc_projection = int( pdmc_projection_time/time_step)
|
|
pdmc_projection_step(:) = 0
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection+1,pdmc_n_diag) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Population weight of PDMC
|
|
END_DOC
|
|
pdmc_pop_weight(:,:) = 1.d0
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (pdmc_n_diag) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Population weight of PDMC
|
|
END_DOC
|
|
pdmc_pop_weight_mult(:) = pdmc_pop_weight(pdmc_projection,:)
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ integer, pdmc_n_diag ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Size of the matrix to diagonalize
|
|
END_DOC
|
|
pdmc_n_diag = 2
|
|
END_PROVIDER
|
|
|
|
|
|
|