From bb774c319fc96b6f75a946e488d808d155eaf655 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 3 May 2016 21:10:25 +0200 Subject: [PATCH] Separated PDMC and SRMC --- src/PROPERTIES/properties_general.irp.f | 15 ++++--- src/SAMPLING/pdmc_step.irp.f | 52 +++++++++++++------------ src/SAMPLING/srmc_step.irp.f | 32 +++++++-------- 3 files changed, 53 insertions(+), 46 deletions(-) diff --git a/src/PROPERTIES/properties_general.irp.f b/src/PROPERTIES/properties_general.irp.f index 3ef9e2f..1cf0bdc 100644 --- a/src/PROPERTIES/properties_general.irp.f +++ b/src/PROPERTIES/properties_general.irp.f @@ -47,15 +47,20 @@ BEGIN_PROVIDER [ double precision, wf_extension ] SOFT_TOUCH wf_extension_min wf_extension_max END_PROVIDER -BEGIN_PROVIDER [ double precision, srmc_pop_weight ] +BEGIN_PROVIDER [ double precision, pop_weight ] implicit none BEGIN_DOC ! Weight of the SRMC population END_DOC - srmc_pop_weight = pop_weight_mult - srmc_pop_weight_min = min(srmc_pop_weight,srmc_pop_weight_min) - srmc_pop_weight_max = max(srmc_pop_weight,srmc_pop_weight_max) - SOFT_TOUCH srmc_pop_weight_min srmc_pop_weight_max + include '../types.F' + if (qmc_method == t_SRMC) then + pop_weight = srmc_pop_weight_mult + else if (qmc_method == t_PDMC) then + pop_weight = pdmc_pop_weight_mult + endif + pop_weight_min = min(pop_weight,pop_weight_min) + pop_weight_max = max(pop_weight,pop_weight_max) + SOFT_TOUCH pop_weight_min pop_weight_max END_PROVIDER diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index 00eb9a0..de64aa2 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -157,7 +157,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) * pop_weight_mult + elec_coord(elec_num+1,3) = pdmc_weight(i_walk) * pdmc_pop_weight_mult do l=1,3 do i=1,elec_num+1 elec_coord_full(i,l,i_walk) = elec_coord(i,l) @@ -172,8 +172,8 @@ END_SHELL psi_value_save(i_walk) = psi_value E_loc_save(i_walk) = E_loc - if (dabs(pdmc_weight(i_walk)*pop_weight_mult) > 1.d-6) then - dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pop_weight_mult) + if (dabs(pdmc_weight(i_walk)*pdmc_pop_weight_mult) > 1.d-6) then + dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult) else dmc_zv_weight = 0.d0 endif @@ -184,18 +184,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 * pop_weight_mult * pdmc_weight(i_walk) - ! $X_2_pdmc_block_walk += $X_2 * pop_weight_mult * pdmc_weight(i_walk) + ! $X_pdmc_block_walk += $X * pdmc_pop_weight_mult * pdmc_weight(i_walk) + ! $X_2_pdmc_block_walk += $X_2 * pdmc_pop_weight_mult * pdmc_weight(i_walk) ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm - $X_pdmc_block_walk_kahan($D2 3) = $X * pop_weight_mult * pdmc_weight(i_walk) - $X_pdmc_block_walk_kahan($D2 1) + $X_pdmc_block_walk_kahan($D2 3) = $X * pdmc_pop_weight_mult * 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 * pop_weight_mult * 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_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) @@ -213,10 +213,12 @@ for p in properties: END_SHELL - block_weight += pop_weight_mult * pdmc_weight(i_walk) + block_weight += pdmc_pop_weight_mult * pdmc_weight(i_walk) else - pdmc_weight(i_walk) = 0.d0 + pdmc_weight(i_walk) = 1.d0 + pdmc_pop_weight(:) = 1.d0 + pdmc_pop_weight_mult = 1.d0 endif @@ -231,15 +233,15 @@ END_SHELL ! Eventually, recompute the weight of the population if (pdmc_projection_step == 1) then - pop_weight_mult = 1.d0 + pdmc_pop_weight_mult = 1.d0 do k=1,pdmc_projection - pop_weight_mult *= pop_weight(k) + pdmc_pop_weight_mult *= pdmc_pop_weight(k) enddo endif ! Remove contribution of the old value of the weight at the new ! projection step - pop_weight_mult *= 1.d0/pop_weight(pdmc_projection_step) + pdmc_pop_weight_mult *= 1.d0/pdmc_pop_weight(pdmc_projection_step) ! Compute the new weight of the population double precision :: sum_weight @@ -247,10 +249,10 @@ END_SHELL do k=1,walk_num sum_weight += pdmc_weight(k) enddo - pop_weight(pdmc_projection_step) = sum_weight/dble(walk_num) + pdmc_pop_weight(pdmc_projection_step) = sum_weight/dble(walk_num) ! Update the running population weight - pop_weight_mult *= pop_weight(pdmc_projection_step) + pdmc_pop_weight_mult *= pdmc_pop_weight(pdmc_projection_step) call system_clock(cpu1, count_rate, count_max) if (cpu1 < cpu0) then @@ -264,7 +266,7 @@ END_SHELL cpu2 = cpu1 endif - SOFT_TOUCH elec_coord_full pop_weight_mult + SOFT_TOUCH elec_coord_full pdmc_pop_weight_mult first_loop = .False. @@ -291,14 +293,6 @@ END_SHELL END_PROVIDER -BEGIN_PROVIDER [ double precision, pop_weight_mult ] - implicit none - BEGIN_DOC -! Population weight of PDMC - END_DOC - pop_weight_mult = pop_weight(pdmc_projection) -END_PROVIDER - BEGIN_PROVIDER [ integer, pdmc_projection ] &BEGIN_PROVIDER [ integer, pdmc_projection_step ] implicit none @@ -312,12 +306,20 @@ END_PROVIDER pdmc_projection_step = 0 END_PROVIDER -BEGIN_PROVIDER [ double precision, pop_weight, (0:pdmc_projection+1) ] +BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection+1) ] implicit none BEGIN_DOC ! Population weight of PDMC END_DOC - pop_weight(:) = 1.d0 + pdmc_pop_weight(:) = 1.d0 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult ] + implicit none + BEGIN_DOC +! Population weight of PDMC + END_DOC + pdmc_pop_weight_mult = pdmc_pop_weight(pdmc_projection) END_PROVIDER diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f index 57e7210..f97c055 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -154,7 +154,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) = srmc_weight(i_walk) * pop_weight_mult + elec_coord(elec_num+1,3) = srmc_weight(i_walk) * srmc_pop_weight_mult do l=1,3 do i=1,elec_num+1 elec_coord_full(i,l,i_walk) = elec_coord(i,l) @@ -175,18 +175,18 @@ from properties import * t = """ if (calc_$X) then ! Kahan's summation algorithm to compute these sums reducing the rounding error: - ! $X_srmc_block_walk += $X * pop_weight_mult * srmc_weight(i_walk) - ! $X_2_srmc_block_walk += $X_2 * pop_weight_mult * srmc_weight(i_walk) + ! $X_srmc_block_walk += $X * srmc_pop_weight_mult * srmc_weight(i_walk) + ! $X_2_srmc_block_walk += $X_2 * srmc_pop_weight_mult * srmc_weight(i_walk) ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm - $X_srmc_block_walk_kahan($D2 3) = $X * pop_weight_mult * srmc_weight(i_walk) - $X_srmc_block_walk_kahan($D2 1) + $X_srmc_block_walk_kahan($D2 3) = $X * srmc_pop_weight_mult * srmc_weight(i_walk) - $X_srmc_block_walk_kahan($D2 1) $X_srmc_block_walk_kahan($D2 2) = $X_srmc_block_walk $D1 + $X_srmc_block_walk_kahan($D2 3) $X_srmc_block_walk_kahan($D2 1) = ($X_srmc_block_walk_kahan($D2 2) - $X_srmc_block_walk $D1 ) & - $X_srmc_block_walk_kahan($D2 3) $X_srmc_block_walk $D1 = $X_srmc_block_walk_kahan($D2 2) - $X_2_srmc_block_walk_kahan($D2 3) = $X_2 * pop_weight_mult * srmc_weight(i_walk) - $X_2_srmc_block_walk_kahan($D2 1) + $X_2_srmc_block_walk_kahan($D2 3) = $X_2 * srmc_pop_weight_mult * srmc_weight(i_walk) - $X_2_srmc_block_walk_kahan($D2 1) $X_2_srmc_block_walk_kahan($D2 2) = $X_2_srmc_block_walk $D1 + $X_2_srmc_block_walk_kahan($D2 3) $X_2_srmc_block_walk_kahan($D2 1) = ($X_2_srmc_block_walk_kahan($D2 2) - $X_2_srmc_block_walk $D1 ) & - $X_2_srmc_block_walk_kahan($D2 3) @@ -204,7 +204,7 @@ for p in properties: END_SHELL - block_weight += pop_weight_mult * srmc_weight(i_walk) + block_weight += srmc_pop_weight_mult * srmc_weight(i_walk) else srmc_weight(i_walk) = 0.d0 @@ -222,15 +222,15 @@ END_SHELL ! Eventually, recompute the weight of the population if (srmc_projection_step == 1) then - pop_weight_mult = 1.d0 + srmc_pop_weight_mult = 1.d0 do k=1,srmc_projection - pop_weight_mult *= pop_weight(k) + srmc_pop_weight_mult *= srmc_pop_weight(k) enddo endif ! Remove contribution of the old value of the weight at the new ! projection step - pop_weight_mult *= 1.d0/pop_weight(srmc_projection_step) + srmc_pop_weight_mult *= 1.d0/srmc_pop_weight(srmc_projection_step) ! Compute the new weight of the population double precision :: sum_weight @@ -238,10 +238,10 @@ END_SHELL do k=1,walk_num sum_weight += srmc_weight(k) enddo - pop_weight(srmc_projection_step) = sum_weight/dble(walk_num) + srmc_pop_weight(srmc_projection_step) = sum_weight/dble(walk_num) ! Update the running population weight - pop_weight_mult *= pop_weight(srmc_projection_step) + srmc_pop_weight_mult *= srmc_pop_weight(srmc_projection_step) ! Reconfiguration integer :: ipos(walk_num) @@ -288,7 +288,7 @@ END_SHELL cpu2 = cpu1 endif - SOFT_TOUCH elec_coord_full pop_weight_mult + SOFT_TOUCH elec_coord_full srmc_pop_weight_mult first_loop = .False. @@ -315,12 +315,12 @@ END_SHELL END_PROVIDER -BEGIN_PROVIDER [ double precision, pop_weight_mult ] +BEGIN_PROVIDER [ double precision, srmc_pop_weight_mult ] implicit none BEGIN_DOC ! Population weight of SRMC END_DOC - pop_weight_mult = pop_weight(srmc_projection) + srmc_pop_weight_mult = srmc_pop_weight(srmc_projection) END_PROVIDER BEGIN_PROVIDER [ integer, srmc_projection ] @@ -336,12 +336,12 @@ END_PROVIDER srmc_projection_step = 0 END_PROVIDER -BEGIN_PROVIDER [ double precision, pop_weight, (0:srmc_projection+1) ] +BEGIN_PROVIDER [ double precision, srmc_pop_weight, (0:srmc_projection+1) ] implicit none BEGIN_DOC ! Population weight of SRMC END_DOC - pop_weight = 1.d0 + srmc_pop_weight = 1.d0 END_PROVIDER