From 1a001e620b9421d538e9682d5df344e338528e26 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 24 Jun 2016 09:11:37 +0200 Subject: [PATCH] Better time step error in SRMC --- src/SAMPLING/reconfigure.irp.f | 2 +- src/SAMPLING/srmc_step.irp.f | 94 ++++++++++++++++++++++++---------- 2 files changed, 68 insertions(+), 28 deletions(-) diff --git a/src/SAMPLING/reconfigure.irp.f b/src/SAMPLING/reconfigure.irp.f index c9434fa..6c92502 100644 --- a/src/SAMPLING/reconfigure.irp.f +++ b/src/SAMPLING/reconfigure.irp.f @@ -15,7 +15,7 @@ subroutine reconfigure(ipos,w) tmp = tmp + w(k) enddo dwalk_num = dble(walk_num)/tmp - + integer :: kp, km kp=0 km=0 diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f index fc98ab5..5e06331 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -44,8 +44,8 @@ END_SHELL 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 :: E_loc_save(4,walk_num_dmc_max) + double precision :: E_loc_save_tmp(4,walk_num_dmc_max) double precision :: psi_value_save(walk_num) double precision :: psi_value_save_tmp(walk_num) double precision :: srmc_weight(walk_num) @@ -124,7 +124,7 @@ END_SHELL 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) + E_loc = E_loc_save(1,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 @@ -135,7 +135,7 @@ END_SHELL enddo TOUCH elec_coord psi_value_save(i_walk) = psi_value - E_loc_save(i_walk) = E_loc + E_loc_save(:,i_walk) = E_loc endif double precision :: p,q @@ -144,19 +144,54 @@ END_SHELL 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 - srmc_weight(i_walk) = dexp(-dtime_step*thr) - else if ( delta < -thr ) then - srmc_weight(i_walk) = dexp(dtime_step*thr) - else +! delta = (E_loc+E_loc_save(1,i_walk))*0.5d0 +! delta = (5.d0 * E_loc + 8.d0 * E_loc_save(1,i_walk) - E_loc_save(2,i_walk))/12.d0 + + delta = (9.d0*E_loc+19.d0*E_loc_save(1,i_walk)- & + 5.d0*E_loc_save(2,i_walk)+E_loc_save(3,i_walk))/24.d0 + +! delta = -((-251.d0*E_loc)-646.d0*E_loc_save(1,i_walk)+264.d0*E_loc_save(2,i_walk)-& +! 106.d0*E_loc_save(3,i_walk)+19.d0*E_loc_save(4,i_walk))/720.d0 + + delta = (delta - E_ref)*p + + if (delta >= 0.d0) then srmc_weight(i_walk) = dexp(-dtime_step*delta) + else + srmc_weight(i_walk) = 2.d0-dexp(dtime_step*delta) endif + + +! if (accepted) then +! ! Compute correction to past weights +! double precision :: delta_old, delta_new +! delta_old = (9.d0*E_loc_save(1,i_walk)+19.d0*E_loc_save(2,i_walk)-& +! 5.d0*E_loc_save(3,i_walk)+E_loc_save(4,i_walk))/24.d0 - E_ref +! +! +! if (delta_old >= 0.d0) then +! srmc_weight(i_walk) = srmc_weight(i_walk) * dexp(dtime_step*delta_old) +! else +! srmc_weight(i_walk) = srmc_weight(i_walk) * (2.d0-dexp(-dtime_step*delta_old)) +! endif +! +! delta_new = (-(E_loc_save_tmp(3,i_walk)-13.d0*E_loc_save_tmp(2,i_walk)& +! -13.d0*E_loc_save_tmp(1,i_walk)+E_loc))/24.d0 - E_ref +! +! if (delta_new >= 0.d0) then +! srmc_weight(i_walk) = srmc_weight(i_walk) * dexp(-dtime_step*delta_new) +! else +! srmc_weight(i_walk) = srmc_weight(i_walk) * (2.d0-dexp(dtime_step*delta_new) ) +! endif +! +! 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 do l=1,3 do i=1,elec_num+1 + elec_coord_full(i,l,i_walk) = elec_coord(i,l) enddo enddo @@ -167,25 +202,30 @@ END_SHELL enddo psi_value_save(i_walk) = psi_value - E_loc_save(i_walk) = E_loc + if (accepted) then + E_loc_save(4,i_walk) = E_loc_save(3,i_walk) + E_loc_save(3,i_walk) = E_loc_save(2,i_walk) + E_loc_save(2,i_walk) = E_loc_save(1,i_walk) + E_loc_save(1,i_walk) = E_loc + endif 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_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) @@ -203,7 +243,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 @@ -221,15 +261,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 @@ -237,10 +277,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) @@ -257,7 +297,7 @@ END_SHELL enddo enddo psi_value_save_tmp(k) = psi_value_save(k) - E_loc_save_tmp(k) = E_loc_save(k) + E_loc_save_tmp(:,k) = E_loc_save(:,k) enddo integer :: ipm @@ -272,7 +312,7 @@ END_SHELL enddo enddo psi_value_save(k) = psi_value_save_tmp(ipm) - E_loc_save(k) = E_loc_save_tmp(ipm) + E_loc_save(:,k) = E_loc_save_tmp(:,ipm) enddo call system_clock(cpu1, count_rate, count_max) @@ -287,7 +327,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. @@ -314,7 +354,7 @@ 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 @@ -335,12 +375,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