LCPQ
/
qmcchem
mirror of https://gitlab.com/scemama/qmcchem.git
10
1
Fork 1
Code Releases Activity

Generalized n_diag

This commit is contained in:
Anthony Scemama 2016-05-07 00:58:45 +02:00
parent 3b94f15bce
commit ca21d44409
1 changed files with 34 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
src/SAMPLING/pdmc_step.irp.f
View File

@ -111,7 +111,7 @@ END_SHELL
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)
double precision :: H(0:pdmc_n_diag/2,0:pdmc_n_diag/2), S(0:pdmc_n_diag/2,0:pdmc_n_diag/2), w(0:pdmc_n_diag/2), work(3*pdmc_n_diag+1)
H = 0.d0
S = 0.d0
@ -161,7 +161,7 @@ END_SHELL
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)
elec_coord(elec_num+1,3) = pdmc_weight(i_walk) * pdmc_pop_weight_mult(pdmc_n_diag)
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,i_walk) = elec_coord(i,l)
@ -176,9 +176,9 @@ END_SHELL
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))
if (dabs(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag)) > 1.d-15) then
dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag))
dmc_zv_weight_half = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag/2))
else
dmc_zv_weight = 0.d0
dmc_zv_weight_half = 0.d0
@ -190,18 +190,18 @@ 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)
! $X_pdmc_block_walk += $X * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk)
! $X_2_pdmc_block_walk += $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * 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 3) = $X * pdmc_pop_weight_mult(pdmc_n_diag) * 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 3) = $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * 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)
@ -219,25 +219,21 @@ for p in properties:
END_SHELL
block_weight += pdmc_pop_weight_mult(1) * pdmc_weight(i_walk)
block_weight += pdmc_pop_weight_mult(pdmc_n_diag) * 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(0,0) += E_loc
S(0,0) += 1.d0
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)
pdmc_pop_weight_mult(0) = 1.d0/pdmc_weight(i_walk)
do k=0,pdmc_n_diag/2
do l=0,pdmc_n_diag/2
H(k,l) += E_loc*pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
S(k,l) += pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
enddo
enddo
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
@ -247,8 +243,8 @@ END_SHELL
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
if (pdmc_projection(pdmc_n_diag) > 0) then
pdmc_projection_step(k) = mod(pdmc_projection_step(k),pdmc_projection(k))+1
else
pdmc_projection_step(k) = 1
endif
@ -256,7 +252,7 @@ END_SHELL
! 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
do l=1,pdmc_projection(k)
pdmc_pop_weight_mult(k) *= pdmc_pop_weight(l,k)
enddo
endif
@ -308,9 +304,9 @@ for p in properties:
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
call dsygv(1, 'N', 'U', pdmc_n_diag/2+1, H, pdmc_n_diag/2+1, S, pdmc_n_diag/2+1, w, work, 3*(pdmc_n_diag+1), info)
E_loc_zv_diag_pdmc_block_walk = w(0)
print *, info, 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) &
@ -324,7 +320,7 @@ END_SHELL
END_PROVIDER
BEGIN_PROVIDER [ integer, pdmc_projection ]
BEGIN_PROVIDER [ integer, pdmc_projection, (pdmc_n_diag) ]
&BEGIN_PROVIDER [ integer, pdmc_projection_step, (pdmc_n_diag) ]
implicit none
BEGIN_DOC
@ -333,11 +329,15 @@ END_PROVIDER
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(pdmc_n_diag) = int( pdmc_projection_time/time_step)
integer :: k
do k=1,pdmc_n_diag-1
pdmc_projection(k) = k*pdmc_projection(pdmc_n_diag)/pdmc_n_diag
enddo
pdmc_projection_step(:) = 0
END_PROVIDER
BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection+1,pdmc_n_diag) ]
BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection(pdmc_n_diag)+1,pdmc_n_diag) ]
implicit none
BEGIN_DOC
! Population weight of PDMC
@ -345,12 +345,12 @@ BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection+1,pdmc_n_
pdmc_pop_weight(:,:) = 1.d0
END_PROVIDER
BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (pdmc_n_diag) ]
BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (0:pdmc_n_diag) ]
implicit none
BEGIN_DOC
! Population weight of PDMC
END_DOC
pdmc_pop_weight_mult(:) = pdmc_pop_weight(pdmc_projection,:)
pdmc_pop_weight_mult(:) = 1.d0
END_PROVIDER