From bb4b75e86b2fd1130d8dc069b02a8a84d24a7983 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 28 Dec 2016 16:57:53 +0100 Subject: [PATCH] Fixed Bug when block is too long --- src/SAMPLING/dmc_step.irp.f | 5 +- src/SAMPLING/fkmc_step.irp.f | 4 +- src/SAMPLING/pdmc_step.irp.f | 380 +++++++++++++++++++++++++++++++++++ src/SAMPLING/srmc_step.irp.f | 4 +- src/SAMPLING/vmc_step.irp.f | 4 +- 5 files changed, 389 insertions(+), 8 deletions(-) create mode 100644 src/SAMPLING/pdmc_step.irp.f diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f index 2764609..7044080 100644 --- a/src/SAMPLING/dmc_step.irp.f +++ b/src/SAMPLING/dmc_step.irp.f @@ -248,14 +248,15 @@ END_SHELL if (cpu1 < cpu0) then cpu1 = cpu1+cpu0 endif - loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate) + loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time if (cpu1-cpu2 > count_rate) then integer :: do_run call get_running(do_run) - loop = do_run == t_Running + loop = loop.and.(do_run == t_Running) cpu2 = cpu1 endif + SOFT_TOUCH elec_coord_full_dmc psi_value psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z elec_coord enddo diff --git a/src/SAMPLING/fkmc_step.irp.f b/src/SAMPLING/fkmc_step.irp.f index 878af69..3f57d20 100644 --- a/src/SAMPLING/fkmc_step.irp.f +++ b/src/SAMPLING/fkmc_step.irp.f @@ -319,11 +319,11 @@ END_SHELL if (cpu1 < cpu0) then cpu1 = cpu1+cpu0 endif - loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate) + loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time if (cpu1-cpu2 > count_rate) then integer :: do_run call get_running(do_run) - loop = do_run == t_Running + loop = loop.and.(do_run == t_Running) cpu2 = cpu1 endif diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f new file mode 100644 index 0000000..e4c3e70 --- /dev/null +++ b/src/SAMPLING/pdmc_step.irp.f @@ -0,0 +1,380 @@ +! 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(4,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 :: 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(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 + + 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(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 + 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 + +!2 delta = (E_loc+E_loc_save(1,i_walk))*0.5d0 +!3 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 + pdmc_weight(i_walk) = dexp(-dtime_step*delta) + else + pdmc_weight(i_walk) = 2.d0-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(pdmc_n_diag) + 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(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 + + 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 + endif + TOUCH dmc_zv_weight dmc_zv_weight_half + +! do i=1,pdmc_n_diag+1 +! E_loc_zv(i) = E_loc * pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight + (E_trial-E_loc) * dmc_zv_weight +! E_loc_zv(i+pdmc_n_diag+1) = pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight +! enddo + +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(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(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(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) + $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(pdmc_n_diag) * pdmc_weight(i_walk) + + 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) + +! 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(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 + + ! 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(walk_num)/dble(count_rate) < block_time + if (cpu1-cpu2 > count_rate) then + integer :: do_run + call get_running(do_run) + loop = loop.and.(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 + +! H(0,0) = H(3,3) +! H(1,0) = H(4,3) +! H(0,1) = H(3,4) +! H(1,1) = H(4,4) +! S(0,0) = S(3,3) +! S(1,0) = S(4,3) +! S(0,1) = S(3,4) +! S(1,1) = S(4,4) +! +! print *, H(0,0)/S(0,0) +! print *, H(1,1)/S(1,1) +! print *, '' +! +! 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) +! call dsygv(1, 'N', 'U', 2, 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 *, w + + deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp ) + +END_PROVIDER + + + BEGIN_PROVIDER [ integer, pdmc_projection, (pdmc_n_diag) ] +&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(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(pdmc_n_diag)+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, (0:pdmc_n_diag) ] + implicit none + BEGIN_DOC +! Population weight of PDMC + END_DOC + pdmc_pop_weight_mult(:) = 1.d0 +END_PROVIDER + + +BEGIN_PROVIDER [ integer, pdmc_n_diag ] + implicit none + BEGIN_DOC +! Size of the matrix to diagonalize + END_DOC + pdmc_n_diag = 8 +END_PROVIDER + + + diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f index 836ddae..1300fe2 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -319,11 +319,11 @@ END_SHELL if (cpu1 < cpu0) then cpu1 = cpu1+cpu0 endif - loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate) + loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time if (cpu1-cpu2 > count_rate) then integer :: do_run call get_running(do_run) - loop = do_run == t_Running + loop = loop.and.(do_run == t_Running) cpu2 = cpu1 endif diff --git a/src/SAMPLING/vmc_step.irp.f b/src/SAMPLING/vmc_step.irp.f index 7aedf42..e09ee80 100644 --- a/src/SAMPLING/vmc_step.irp.f +++ b/src/SAMPLING/vmc_step.irp.f @@ -132,11 +132,11 @@ END_SHELL if (cpu1 < cpu0) then cpu1 = cpu1+cpu0 endif - loop = dble(cpu1-cpu0)*dble(walk_num) < dble(block_time)*dble(count_rate) + loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time if (cpu1-cpu2 > count_rate) then integer :: do_run call get_running(do_run) - loop = do_run == t_Running + loop = loop.and.(do_run == t_Running) cpu2 = cpu1 endif