From ca21d44409a6d8cc931e06393451686a194eed78 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 7 May 2016 00:58:45 +0200 Subject: [PATCH] Generalized n_diag --- src/SAMPLING/pdmc_step.irp.f | 68 ++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index eedd919..dc18daa 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -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