diff --git a/ocaml/Default.ml b/ocaml/Default.ml index 65b542d..12f5f86 100644 --- a/ocaml/Default.ml +++ b/ocaml/Default.ml @@ -3,7 +3,7 @@ open Core.Std;; let simulation_nucl_fitcusp_factor = lazy( let default = - 0.5 + 1. in if (Ezfio.has_pseudo_do_pseudo ()) then if (Ezfio.get_pseudo_do_pseudo ()) then @@ -14,8 +14,8 @@ let simulation_nucl_fitcusp_factor = lazy( default ) -let electrons_elec_walk_num = lazy ( 30 ) -let electrons_elec_walk_num_tot = lazy ( 10000 ) +let electrons_elec_walk_num = lazy ( 100 ) +let electrons_elec_walk_num_tot = lazy ( 1000 ) let jastrow_jast_type = lazy ( "None" ) let simulation_block_time = lazy ( 30 ) let simulation_ci_threshold = lazy ( 1.e-8 ) diff --git a/ocaml/Input.ml b/ocaml/Input.ml index ad9c770..01e8f97 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -387,7 +387,7 @@ end module Method : sig - type t = VMC | DMC | SRMC + type t = VMC | DMC | SRMC | FKMC val doc : string val read : unit -> t val write : t -> unit @@ -396,21 +396,23 @@ module Method : sig end = struct - type t = VMC | DMC | SRMC + type t = VMC | DMC | SRMC | FKMC - let doc = "QMC Method : [ VMC | DMC | SRMC ]" + let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC ]" let of_string = function | "VMC" | "vmc" -> VMC | "DMC" | "dmc" -> DMC | "SRMC" | "srmc" -> SRMC - | x -> failwith ("Method should be [ VMC | DMC | SRMC ], not "^x^".") + | "FKMC" | "fkmc" -> FKMC + | x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC ], not "^x^".") let to_string = function | VMC -> "VMC" | DMC -> "DMC" | SRMC -> "SRMC" + | FKMC -> "FKMC" let read () = @@ -851,10 +853,12 @@ let validate () = let () = match (sampling, meth, Pseudo.to_bool do_pseudo) with | (Sampling.Brownian, Method.DMC, true) + | (Sampling.Brownian, Method.FKMC, true) | (Sampling.Brownian, Method.SRMC, true) -> if ( (Time_step.to_float ts) >= 0.5 ) then warn ( "Time step seems large for "^(Method.to_string meth) ) | (Sampling.Brownian, Method.SRMC, false) + | (Sampling.Brownian, Method.FKMC, false) | (Sampling.Brownian, Method.DMC, false) -> if ( (Time_step.to_float ts) >= 0.01 ) then warn ( "Time step seems large for "^(Method.to_string meth) ) @@ -865,8 +869,9 @@ let validate () = if ( (Time_step.to_float ts) <= 0.01 ) then warn "Time step seems small for Langevin sampling." | (Sampling.Langevin, Method.SRMC, _) + | (Sampling.Langevin, Method.FKMC, _) | (Sampling.Langevin, Method.DMC, _) -> - failwith "Lanvegin sampling is incompatible with DMC/SRMC" + failwith "Lanvegin sampling is incompatible with DMC" in @@ -874,6 +879,7 @@ let validate () = let () = match (meth, Ref_energy.(read () |> to_float) ) with | (Method.SRMC,0.) + | (Method.FKMC,0.) | (Method.DMC,0.) -> failwith ("E_ref should not be zero in "^(Method.to_string meth) ) | _ -> () in @@ -888,6 +894,7 @@ let validate () = let () = match (meth, Property.(calc E_loc)) with | (Method.SRMC, false) + | (Method.FKMC, false) | (Method.DMC, false) -> failwith ( "E_loc should be sampled in "^(Method.to_string meth) ) | (Method.VMC, false) -> warn "Sampling of E_loc is not activated in input" | _ -> () diff --git a/src/SAMPLING/block.irp.f b/src/SAMPLING/block.irp.f index 092eb17..a732941 100644 --- a/src/SAMPLING/block.irp.f +++ b/src/SAMPLING/block.irp.f @@ -30,6 +30,12 @@ t = """ $X_block_walk = $X_srmc_block_walk $X_2_block_walk = $X_2_srmc_block_walk endif + else if (qmc_method == t_FKMC) then + PROVIDE E_loc_fkmc_block_walk + if (calc_$X) then + $X_block_walk = $X_fkmc_block_walk + $X_2_block_walk = $X_2_fkmc_block_walk + endif endif END_PROVIDER diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f index e11662e..2764609 100644 --- a/src/SAMPLING/dmc_step.irp.f +++ b/src/SAMPLING/dmc_step.irp.f @@ -195,7 +195,7 @@ END_SHELL do k=1,walk_num_dmc sum_weight += dmc_weight(k) enddo - E0 = E_ref - log(real(walk_num_dmc)/real(walk_num)) * 0.1d0/dtime_step + E0 = E_ref - log(sum_weight/real(walk_num)) * 0.1d0 /dtime_step ! Branching integer :: ipos(walk_num_dmc_max), walk_num_dmc_new @@ -214,20 +214,22 @@ END_SHELL ipos(k) = k enddo + walk_num_dmc_new = walk_num_dmc do k=1,walk_num_dmc r = qmc_ranf() if (dmc_weight(k) > 1.d0) then if ( 1.d0+r < dmc_weight(k) ) then - walk_num_dmc = walk_num_dmc+1 - ipos(walk_num_dmc) = k + walk_num_dmc_new = walk_num_dmc_new+1 + ipos(walk_num_dmc_new) = k endif else if ( r > dmc_weight(k) ) then - ipos(k) = ipos(walk_num_dmc) - walk_num_dmc = walk_num_dmc-1 + ipos(k) = ipos(walk_num_dmc_new) + walk_num_dmc_new = walk_num_dmc_new-1 endif endif enddo + walk_num_dmc = walk_num_dmc_new integer :: ipm do k=1,walk_num_dmc @@ -328,7 +330,7 @@ BEGIN_PROVIDER [ integer, walk_num_dmc_max ] BEGIN_DOC ! Max number of walkers in DMC END_DOC - walk_num_dmc_max = 3 * walk_num + walk_num_dmc_max = max(3 * walk_num, 30) END_PROVIDER diff --git a/src/SAMPLING/fkmc_step.irp.f b/src/SAMPLING/fkmc_step.irp.f new file mode 100644 index 0000000..878af69 --- /dev/null +++ b/src/SAMPLING/fkmc_step.irp.f @@ -0,0 +1,375 @@ +! Providers of *_fkmc_block_walk +!============================== +BEGIN_SHELL [ /usr/bin/python ] +from properties import * + +t = """ + BEGIN_PROVIDER [ $T, $X_fkmc_block_walk $D1 ] +&BEGIN_PROVIDER [ $T, $X_fkmc_block_walk_kahan $D2 ] +&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk $D1 ] +&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk_kahan $D2 ] + implicit none + BEGIN_DOC +! fkMC averages of $X. Computed in E_loc_fkmc_block_walk + END_DOC + $X_fkmc_block_walk = 0.d0 + $X_fkmc_block_walk_kahan = 0.d0 + $X_2_fkmc_block_walk = 0.d0 + $X_2_fkmc_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_fkmc_block_walk ] +&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk ] +&BEGIN_PROVIDER [ double precision, E_loc_fkmc_block_walk_kahan, (3) ] +&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk_kahan, (3) ] + implicit none + include '../types.F' + BEGIN_DOC +! Properties averaged over the block using the FKMC method + END_DOC + + integer, parameter :: BIRTH=1, DEATH=2 + 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 :: psi_value_save(walk_num) + double precision :: psi_value_save_tmp(walk_num) + double precision :: fkmc_weight(walk_num) + double precision :: delta(walk_num) + double precision, allocatable :: psi_grad_psi_inv_save(:,:,:) + double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:) + double precision, allocatable :: fkmc_clock_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 :: E_loc_save_tmp + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_weight + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_clock_tmp + 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), & + fkmc_clock_tmp(2,walk_num) ) + psi_value_save = 0.d0 + psi_value_save_tmp = 0.d0 + fkmc_weight = 1.d0 + +! Initialization + if (vmc_algo /= t_Brownian) then + call abrt(irp_here,'FKMC 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_fkmc_block_walk = 0.d0 + !DIR$ VECTOR ALIGNED + $X_fkmc_block_walk_kahan = 0.d0 + !DIR$ VECTOR ALIGNED + $X_2_fkmc_block_walk = 0.d0 + !DIR$ VECTOR ALIGNED + $X_2_fkmc_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 :: thr + + thr = 2.d0/time_step_sq + + logical :: first_loop + first_loop = .True. + + do while (loop) + + ! Every walker makes a step + do i_walk=1,walk_num + + 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(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 + E_loc_save(i_walk) = E_loc + psi_value_save(i_walk) = psi_value + 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 + delta(i_walk) = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p + if ( delta(i_walk) > thr ) then + delta(i_walk) = thr + else if ( delta(i_walk) < -thr ) then + delta(i_walk) = -thr + endif + fkmc_weight(i_walk) = dexp(-dtime_step*delta(i_walk)) + elec_coord(elec_num+1,1) += p*time_step + elec_coord(elec_num+1,2) = E_loc + elec_coord(elec_num+1,3) = fkmc_weight(i_walk) + 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(i_walk) = E_loc + +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_fkmc_block_walk += $X * fkmc_weight(i_walk) + ! $X_2_fkmc_block_walk += $X_2 * fkmc_weight(i_walk) + ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm + + $X_fkmc_block_walk_kahan($D2 3) = $X * fkmc_weight(i_walk) - $X_fkmc_block_walk_kahan($D2 1) + $X_fkmc_block_walk_kahan($D2 2) = $X_fkmc_block_walk $D1 + $X_fkmc_block_walk_kahan($D2 3) + $X_fkmc_block_walk_kahan($D2 1) = ($X_fkmc_block_walk_kahan($D2 2) - $X_fkmc_block_walk $D1 ) & + - $X_fkmc_block_walk_kahan($D2 3) + $X_fkmc_block_walk $D1 = $X_fkmc_block_walk_kahan($D2 2) + + + $X_2_fkmc_block_walk_kahan($D2 3) = $X_2 * fkmc_weight(i_walk) - $X_2_fkmc_block_walk_kahan($D2 1) + $X_2_fkmc_block_walk_kahan($D2 2) = $X_2_fkmc_block_walk $D1 + $X_2_fkmc_block_walk_kahan($D2 3) + $X_2_fkmc_block_walk_kahan($D2 1) = ($X_2_fkmc_block_walk_kahan($D2 2) - $X_2_fkmc_block_walk $D1 ) & + - $X_2_fkmc_block_walk_kahan($D2 3) + $X_2_fkmc_block_walk $D1 = $X_2_fkmc_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 += fkmc_weight(i_walk) + + else + fkmc_weight(i_walk) = 0.d0 + delta(i_walk) = 1.d5 + endif + + enddo + + ! Compute the new weight of the population + double precision :: sum_weight + sum_weight = 0.d0 + do k=1,walk_num + sum_weight += fkmc_weight(k) + enddo + + do k=1,walk_num + do l=1,3 + do i=1,elec_num+1 + elec_coord_tmp(i,l,k) = elec_coord_full(i,l,k) + enddo + do i=1,elec_num + psi_grad_psi_inv_save_tmp(i,l,k) = psi_grad_psi_inv_save(i,l,k) + enddo + enddo + psi_value_save_tmp(k) = psi_value_save(k) + E_loc_save_tmp(k) = E_loc_save(k) + if (fkmc_weight(k) == 0.d0) then + fkmc_clock(DEATH,k) = -1.d0 + endif + if ( delta(k) <= 0.d0 ) then + fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k) +time_step * delta(k) + fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k) + else + fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k) + fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k) -time_step * delta(k) + endif + enddo + +! Reconfiguration +! =============== + + ! Identify first which walkers will be killed to place branched walkers there + ! later + + double precision, external :: qmc_ranf + integer :: ipm, m + integer :: killed(walk_num) + + m=1 + do k=1,walk_num + fkmc_clock(DEATH,k) = fkmc_clock_tmp(DEATH,k) + if (fkmc_clock_tmp(DEATH,k) <= 0.d0) then + killed(m) = k + m += 1 + fkmc_clock(DEATH,k) = -dlog(qmc_ranf()) + fkmc_clock(BIRTH,k) = -dlog(qmc_ranf()) + ipm = k + do while (ipm == k) + ipm = 1 + int (walk_num*qmc_ranf()) + enddo + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm) + enddo + do i=1,elec_num + psi_grad_psi_inv_save(i,l,k) = psi_grad_psi_inv_save_tmp(i,l,ipm) + enddo + enddo + psi_value_save(k) = psi_value_save_tmp(ipm) + E_loc_save(k) = E_loc_save_tmp(ipm) + endif + enddo + killed(m) = 0 + + m=1 + do k=1,walk_num + fkmc_clock(BIRTH,k) = fkmc_clock_tmp(BIRTH,k) + if (fkmc_clock_tmp(BIRTH,k) <= 0.d0) then + fkmc_clock(BIRTH,k) = -dlog(qmc_ranf()) + if (killed(m) == 0) then + ipm = k + do while (ipm == k) + ipm = 1 + int (walk_num*qmc_ranf()) + enddo + else + ipm = killed(m) + m +=1 + endif + fkmc_clock(BIRTH,ipm) = -dlog(qmc_ranf()) + fkmc_clock(DEATH,ipm) = -dlog(qmc_ranf()) + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,ipm) = elec_coord_tmp(i,l,k) + enddo + do i=1,elec_num + psi_grad_psi_inv_save(i,l,ipm) = psi_grad_psi_inv_save_tmp(i,l,k) + enddo + enddo + psi_value_save(ipm) = psi_value_save_tmp(k) + E_loc_save(ipm) = E_loc_save_tmp(k) + + endif + + enddo + + + + call system_clock(cpu1, count_rate, count_max) + if (cpu1 < cpu0) then + cpu1 = cpu1+cpu0 + endif + loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate) + if (cpu1-cpu2 > count_rate) then + integer :: do_run + call get_running(do_run) + loop = do_run == t_Running + cpu2 = cpu1 + endif + +! Update E_ref to take into account the weight of the population + E_ref -= dlog(sum_weight / dble(walk_num) ) / time_step + SOFT_TOUCH elec_coord_full E_ref + + 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_fkmc_block_walk *= factor + $X_2_fkmc_block_walk *= factor + endif +""" +for p in properties: + print t.replace("$X",p[1]) +END_SHELL + + deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp, & + fkmc_clock_tmp ) + +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, fkmc_clock, (2,walk_num) ] + implicit none + BEGIN_DOC + ! Branching clocks for the FKMC algotithm. (1,:) is the birth clock and + ! (2,:) is the death clock. + END_DOC + integer :: i + double precision, external :: qmc_ranf + do i=1, walk_num + fkmc_clock(1,i) = -dlog(qmc_ranf()) + fkmc_clock(2,i) = -dlog(qmc_ranf()) + enddo + +END_PROVIDER + diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f index 3afaab9..fc98ab5 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -133,7 +133,9 @@ END_SHELL elec_coord(i,l) = elec_coord_full(i,l,i_walk) enddo enddo - SOFT_TOUCH elec_coord + TOUCH elec_coord + psi_value_save(i_walk) = psi_value + E_loc_save(i_walk) = E_loc endif double precision :: p,q diff --git a/src/mo.irp.f b/src/mo.irp.f index 2dfa9c6..5e43659 100644 --- a/src/mo.irp.f +++ b/src/mo.irp.f @@ -47,7 +47,6 @@ END_PROVIDER END_DOC mo_scale = 1.d0/(0.4d0*log(float(elec_num+1))) mo_norm = mo_scale*mo_scale - END_PROVIDER @@ -536,7 +535,7 @@ BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num) TOUCH point nucl_fitcusp_param primitives_reduced mo_coef r = point(3) f = dfour_pi*r*r*dr - do i=1,mo_tot_num + do i=1,mo_num t = 0.d0 do j=1,ao_num if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then diff --git a/src/simulation.irp.f b/src/simulation.irp.f index 3608e08..0995e2d 100644 --- a/src/simulation.irp.f +++ b/src/simulation.irp.f @@ -148,7 +148,7 @@ BEGIN_PROVIDER [ integer, qmc_method ] implicit none include 'types.F' BEGIN_DOC - ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC + ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC, t_FKMC END_DOC character*(32) :: method method = types(t_VMC) @@ -160,8 +160,10 @@ BEGIN_PROVIDER [ integer, qmc_method ] qmc_method = t_DMC else if (method == types(t_SRMC)) then qmc_method = t_SRMC + else if (method == types(t_FKMC)) then + qmc_method = t_FKMC else - call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC )') + call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC )') endif call cinfo(irp_here,'qmc_method',trim(method)) @@ -288,6 +290,9 @@ BEGIN_PROVIDER [ integer, vmc_algo ] if (qmc_method == t_SRMC) then stop 'Langevin incompatible with SRMC' endif + if (qmc_method == t_FKMC) then + stop 'Langevin incompatible with FKMC' + endif else if (Sampling == types(t_MTM)) then vmc_algo = t_MTM else diff --git a/src/types.F b/src/types.F index 26ce25d..f2684ec 100644 --- a/src/types.F +++ b/src/types.F @@ -6,6 +6,7 @@ integer, parameter :: t_VMC = 7 integer, parameter :: t_DMC = 8 integer, parameter :: t_SRMC = 9 + integer, parameter :: t_FKMC = 10 integer, parameter :: t_Simple = 11 integer, parameter :: t_None = 12 @@ -26,7 +27,7 @@ 'VMC ', & 'DMC ', & 'SRMC ', & - ' ', & + 'FKMC ', & 'Simple ', & 'None ', & ' ', &