From 8cc2c6a24bb7647bf99e3736260249fa746c945b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 2 May 2016 10:52:29 +0200 Subject: [PATCH 01/20] Implemented ZV DMC --- src/PROPERTIES/properties_energy.irp.f | 8 ++++++++ src/SAMPLING/dmc_step.irp.f | 5 +++++ 2 files changed, 13 insertions(+) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index d21f973..e6d40e9 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -248,5 +248,13 @@ BEGIN_PROVIDER [ double precision, E_loc ] END_PROVIDER +BEGIN_PROVIDER [ double precision, E_loc_zv ] + implicit none + BEGIN_DOC + ! Zero-variance parameter on E_loc + END_DOC + E_loc_zv = E_ref - E_loc +END_PROVIDER + diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f index 2764609..0694a68 100644 --- a/src/SAMPLING/dmc_step.irp.f +++ b/src/SAMPLING/dmc_step.irp.f @@ -150,6 +150,11 @@ END_SHELL psi_value_save(i_walk) = psi_value E_loc_save(i_walk) = E_loc + if (calc_E_loc_zv) then + E_loc_zv = (E_ref - E_loc)/dmc_weight(i_walk) + TOUCH E_loc_zv + endif + BEGIN_SHELL [ /usr/bin/python ] from properties import * t = """ From 571df84d9ddec079f559fed2d14de5707dd6722e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 2 May 2016 21:19:36 +0200 Subject: [PATCH 02/20] ZV energy --- ocaml/Input.ml | 32 ++- ocaml/Qmcchem_result.ml | 9 + ocaml/Random_variable.ml | 33 ++- src/PROPERTIES/properties_energy.irp.f | 2 +- src/SAMPLING/block.irp.f | 6 + src/SAMPLING/dmc_step.irp.f | 5 - src/SAMPLING/pdmc_step.irp.f | 325 +++++++++++++++++++++++++ src/SAMPLING/reconfigure.irp.f | 2 +- src/SAMPLING/srmc_step.irp.f | 1 + src/simulation.irp.f | 4 +- src/types.F | 49 ++-- 11 files changed, 416 insertions(+), 52 deletions(-) create mode 100644 src/SAMPLING/pdmc_step.irp.f diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 01e8f97..894754b 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -387,7 +387,7 @@ end module Method : sig - type t = VMC | DMC | SRMC | FKMC + type t = VMC | DMC | SRMC | FKMC | PDMC val doc : string val read : unit -> t val write : t -> unit @@ -396,22 +396,24 @@ module Method : sig end = struct - type t = VMC | DMC | SRMC | FKMC + type t = VMC | DMC | SRMC | FKMC | PDMC - let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC ]" + let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC | PDMC ]" let of_string = function | "VMC" | "vmc" -> VMC | "DMC" | "dmc" -> DMC | "SRMC" | "srmc" -> SRMC + | "PDMC" | "pdmc" -> PDMC | "FKMC" | "fkmc" -> FKMC - | x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC ], not "^x^".") + | x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC | PDMC ], not "^x^".") let to_string = function | VMC -> "VMC" | DMC -> "DMC" | SRMC -> "SRMC" + | PDMC -> "PDMC" | FKMC -> "FKMC" @@ -852,25 +854,19 @@ let validate () = (* Check sampling and time steps *) 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) ) | (Sampling.Brownian, Method.VMC, _) -> if ( (Time_step.to_float ts) >= 10. ) then warn "Time step seems large for VMC." | (Sampling.Langevin, Method.VMC, _) -> 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, _) -> + | (Sampling.Brownian, _, true) -> + if ( (Time_step.to_float ts) >= 0.5 ) then + warn ( "Time step seems large for "^(Method.to_string meth) ) + | (Sampling.Brownian, _, false) -> + if ( (Time_step.to_float ts) >= 0.01 ) then + warn ( "Time step seems large for "^(Method.to_string meth) ) + | (Sampling.Langevin, _, _) -> failwith "Lanvegin sampling is incompatible with DMC" in @@ -879,6 +875,7 @@ let validate () = let () = match (meth, Ref_energy.(read () |> to_float) ) with | (Method.SRMC,0.) + | (Method.PDMC,0.) | (Method.FKMC,0.) | (Method.DMC,0.) -> failwith ("E_ref should not be zero in "^(Method.to_string meth) ) | _ -> () @@ -894,6 +891,7 @@ let validate () = let () = match (meth, Property.(calc E_loc)) with | (Method.SRMC, false) + | (Method.PDMC, 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/ocaml/Qmcchem_result.ml b/ocaml/Qmcchem_result.ml index 23e0660..c8cbc04 100644 --- a/ocaml/Qmcchem_result.ml +++ b/ocaml/Qmcchem_result.ml @@ -145,6 +145,15 @@ let display_summary ~range = in List.iter properties ~f:print_property ; +(** TODO *) + let open Random_variable in + let p = (of_raw_data ~range Property.E_loc_zv) + +! (of_raw_data ~range Property.E_loc) + in + Printf.printf "%20s : %s\n" + ("E_loc_zv(+)") + (Random_variable.to_string p); +(** TODO *) let cpu = Random_variable.of_raw_data ~range Property.Cpu diff --git a/ocaml/Random_variable.ml b/ocaml/Random_variable.ml index 01fc04a..94c46b2 100644 --- a/ocaml/Random_variable.ml +++ b/ocaml/Random_variable.ml @@ -162,6 +162,8 @@ let average { property ; data } = + + (** Compute sum (for CPU/Wall time) *) let sum { property ; data } = List.fold data ~init:0. ~f:(fun accu x -> @@ -342,7 +344,9 @@ let max_block = (** Create a hash table for merging *) -let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t = +let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) + ?(update_value=(fun wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw) ) + ?(update_weight=(fun wc wb -> wc +. wb) ) t = let table = Hashtbl.create ~hashable:hashable () in List.iter t.data ~f:(fun block -> @@ -356,7 +360,7 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t = Weight.to_float block.weight in let sw = - wc +. wb + update_weight wc wb in if (Property.is_scalar current.property) then let vc, vb = @@ -365,7 +369,7 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t = in Some { property = current.property ; weight = Weight.of_float sw ; - value = Sample.of_float ((wc *. vc +. wb *. vb) /. sw); + value = Sample.of_float (update_value wc vc wb vb sw); block_id = update_block_id block.block_id; pid = block.pid ; compute_node = block.compute_node; @@ -380,7 +384,7 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t = { property = current.property ; weight = Weight.of_float sw ; value = - Array.init dim ~f:(fun i -> ((wc *. vc.(i) +. wb *. vb.(i)) /. sw)) + Array.init dim ~f:(fun i -> update_value wc vc.(i) wb vb.(i) sw) |> Sample.of_float_array ~dim ; block_id = update_block_id block.block_id; pid = block.pid ; @@ -401,9 +405,8 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t = (** Genergic merge function *) -let merge ~hashable ~create_key ?update_block_id t = - let table = create_hash ~hashable:hashable ~create_key:create_key - ?update_block_id:update_block_id t +let merge ~hashable ~create_key ?update_block_id ?update_value ?update_weight t = + let table = create_hash ~hashable ~create_key ?update_block_id ?update_value ?update_weight t in { property = t.property ; data = Hashtbl.to_alist table @@ -454,6 +457,20 @@ let merge_per_compute_node_and_block_id = (Block_id.to_int block.Block.block_id) ) +(** Create a new random variable which is a sum of 2 *) +let (+!) p1 p2 = + merge + ~hashable:String.hashable + ~create_key:(fun block -> + Printf.sprintf "%s %10.10d %10.10d" + (Compute_node.to_string block.Block.compute_node) + (Block_id.to_int block.Block.block_id) + (Pid.to_int block.Block.pid) ) + ~update_value: (fun wc vc wb vb sw -> (vc +. vb) ) + ~update_weight:(fun wc wb -> wc ) + { property = p1.property ; + data = List.concat [ p1.data ; p2.data ] } + (** Merge two consecutive blocks *) @@ -467,7 +484,7 @@ let compress = ((Block_id.to_int block_id)+1)/2 |> Block_id.of_int ) - + (** Last value on each compute node (for wall_time) *) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index e6d40e9..727feec 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ] BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC - E_loc_zv = E_ref - E_loc + E_loc_zv = 0.d0 END_PROVIDER diff --git a/src/SAMPLING/block.irp.f b/src/SAMPLING/block.irp.f index a732941..2492e12 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_PDMC) then + PROVIDE E_loc_pdmc_block_walk + if (calc_$X) then + $X_block_walk = $X_pdmc_block_walk + $X_2_block_walk = $X_2_pdmc_block_walk + endif else if (qmc_method == t_FKMC) then PROVIDE E_loc_fkmc_block_walk if (calc_$X) then diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f index 0694a68..2764609 100644 --- a/src/SAMPLING/dmc_step.irp.f +++ b/src/SAMPLING/dmc_step.irp.f @@ -150,11 +150,6 @@ END_SHELL psi_value_save(i_walk) = psi_value E_loc_save(i_walk) = E_loc - if (calc_E_loc_zv) then - E_loc_zv = (E_ref - E_loc)/dmc_weight(i_walk) - TOUCH E_loc_zv - endif - BEGIN_SHELL [ /usr/bin/python ] from properties import * t = """ diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f new file mode 100644 index 0000000..8dd011b --- /dev/null +++ b/src/SAMPLING/pdmc_step.irp.f @@ -0,0 +1,325 @@ +! 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(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 :: 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 :: 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 :: 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 + + 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 + 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 + delta = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p + if ( delta > thr ) then + pdmc_weight(i_walk) = dexp(-dtime_step*thr) + else if ( delta < -thr ) then + pdmc_weight(i_walk) = dexp(dtime_step*thr) + else + pdmc_weight(i_walk) = 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) * 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 + 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 + + if (calc_E_loc_zv) then + if (dabs(pdmc_weight(i_walk)*pop_weight_mult) > 1.d-6) then + E_loc_zv = (E_ref-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult) + else + E_loc_zv = 0.d0 + endif + TOUCH E_loc_zv + 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_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) + ! 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 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 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 += pop_weight_mult * pdmc_weight(i_walk) + + else + pdmc_weight(i_walk) = 0.d0 + endif + + + enddo + + ! Move to the next projection step + if (pdmc_projection > 0) then + pdmc_projection_step = mod(pdmc_projection_step,pdmc_projection)+1 + else + pdmc_projection_step = 1 + endif + + ! Eventually, recompute the weight of the population + if (pdmc_projection_step == 1) then + pop_weight_mult = 1.d0 + do k=1,pdmc_projection + pop_weight_mult *= 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) + + ! Compute the new weight of the population + double precision :: sum_weight + sum_weight = 0.d0 + do k=1,walk_num + sum_weight += pdmc_weight(k) + enddo + pop_weight(pdmc_projection_step) = sum_weight/dble(walk_num) + + ! Update the running population weight + pop_weight_mult *= pop_weight(pdmc_projection_step) + + 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 + + SOFT_TOUCH elec_coord_full 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 + + deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp ) + +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 + 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 = int( pdmc_projection_time/time_step) + pdmc_projection_step = 0 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pop_weight, (0:pdmc_projection+1) ] + implicit none + BEGIN_DOC +! Population weight of PDMC + END_DOC + pop_weight(:) = 1.d0 +END_PROVIDER + + 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..57e7210 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -169,6 +169,7 @@ END_SHELL psi_value_save(i_walk) = psi_value E_loc_save(i_walk) = E_loc + BEGIN_SHELL [ /usr/bin/python ] from properties import * t = """ diff --git a/src/simulation.irp.f b/src/simulation.irp.f index 0995e2d..ceb66d9 100644 --- a/src/simulation.irp.f +++ b/src/simulation.irp.f @@ -162,8 +162,10 @@ BEGIN_PROVIDER [ integer, qmc_method ] qmc_method = t_SRMC else if (method == types(t_FKMC)) then qmc_method = t_FKMC + else if (method == types(t_PDMC)) then + qmc_method = t_PDMC else - call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC )') + call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC | PDMC )') endif call cinfo(irp_here,'qmc_method',trim(method)) diff --git a/src/types.F b/src/types.F index f2684ec..b3e1cac 100644 --- a/src/types.F +++ b/src/types.F @@ -7,30 +7,41 @@ integer, parameter :: t_DMC = 8 integer, parameter :: t_SRMC = 9 integer, parameter :: t_FKMC = 10 + integer, parameter :: t_PDMC = 11 - integer, parameter :: t_Simple = 11 - integer, parameter :: t_None = 12 - integer, parameter :: t_Core = 14 + integer, parameter :: t_Simple = 21 + integer, parameter :: t_None = 22 + integer, parameter :: t_Core = 24 integer, parameter :: t_Stopped = 0 integer, parameter :: t_Queued = 1 integer, parameter :: t_Running = 2 integer, parameter :: t_Stopping = 3 - character*(32) :: types(15) = & - (/ ' ', & - ' ', & - 'Brownian ', & - 'Langevin ', & - ' ', & - ' ', & - 'VMC ', & - 'DMC ', & - 'SRMC ', & - 'FKMC ', & - 'Simple ', & - 'None ', & - ' ', & - 'Core ', & - ' '/) + character*(32) :: types(25) = & + (/ ' ', & + ' ', & + 'Brownian ', & + 'Langevin ', & + ' ', & + ' ', & + 'VMC ', & + 'DMC ', & + 'SRMC ', & + 'FKMC ', & + 'PDMC ', & + ' ', & + ' ', & + ' ', & + ' ', & + ' ', & + ' ', & + ' ', & + ' ', & + ' ', & + 'Simple ', & + 'None ', & + ' ', & + 'Core ', & + ' '/) From 54f2bae5f6405ebdd41144848629ee5d7228a1ae Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 2 May 2016 21:51:09 +0200 Subject: [PATCH 03/20] Added operators on random variables --- ocaml/Qmcchem_result.ml | 8 ++--- ocaml/Random_variable.ml | 45 ++++++++++++++++++++++++-- src/PROPERTIES/properties_energy.irp.f | 2 +- src/SAMPLING/pdmc_step.irp.f | 2 +- 4 files changed, 48 insertions(+), 9 deletions(-) diff --git a/ocaml/Qmcchem_result.ml b/ocaml/Qmcchem_result.ml index c8cbc04..c1737d6 100644 --- a/ocaml/Qmcchem_result.ml +++ b/ocaml/Qmcchem_result.ml @@ -145,15 +145,15 @@ let display_summary ~range = in List.iter properties ~f:print_property ; -(** TODO *) +(* let open Random_variable in - let p = (of_raw_data ~range Property.E_loc_zv) - +! (of_raw_data ~range Property.E_loc) + let p = (of_raw_data ~range Property.E_loc) + +! (of_raw_data ~range Property.E_loc_zv) in Printf.printf "%20s : %s\n" ("E_loc_zv(+)") (Random_variable.to_string p); -(** TODO *) + *) let cpu = Random_variable.of_raw_data ~range Property.Cpu diff --git a/ocaml/Random_variable.ml b/ocaml/Random_variable.ml index 94c46b2..c1ae8ad 100644 --- a/ocaml/Random_variable.ml +++ b/ocaml/Random_variable.ml @@ -457,20 +457,59 @@ let merge_per_compute_node_and_block_id = (Block_id.to_int block.Block.block_id) ) -(** Create a new random variable which is a sum of 2 *) -let (+!) p1 p2 = +(** Create float, variable operators *) +let one_variable_operator ~update_value p f = merge + ~update_value + ~hashable:String.hashable + ~create_key:(fun block -> + Printf.sprintf "%s %10.10d %10.10d" + (Compute_node.to_string block.Block.compute_node) + (Block_id.to_int block.Block.block_id) + (Pid.to_int block.Block.pid) ) + ~update_weight:(fun wc wb -> wc ) + p + +let ( +@ ) p f = one_variable_operator p f + ~update_value: (fun wc vc wb vb sw -> f +. (wc *. vc +. wb *. vb) /. sw) + +let ( *@ ) p f = one_variable_operator p f + ~update_value: (fun wc vc wb vb sw -> f *. (wc *. vc +. wb *. vb) /. sw) + +let ( -@ ) p f = one_variable_operator p f + ~update_value: (fun wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw -. f) + +let ( /@ ) p f = one_variable_operator p f + ~update_value: (fun wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw /. f) + + +(** Create two variable operators *) +let two_variable_operator ~update_value p1 p2 = + merge + ~update_value ~hashable:String.hashable ~create_key:(fun block -> Printf.sprintf "%s %10.10d %10.10d" (Compute_node.to_string block.Block.compute_node) (Block_id.to_int block.Block.block_id) (Pid.to_int block.Block.pid) ) - ~update_value: (fun wc vc wb vb sw -> (vc +. vb) ) ~update_weight:(fun wc wb -> wc ) { property = p1.property ; data = List.concat [ p1.data ; p2.data ] } +let ( +! ) = two_variable_operator + ~update_value: (fun wc vc wb vb sw -> (vc +. vb) ) + +let ( *! ) = two_variable_operator + ~update_value: (fun wc vc wb vb sw -> (vc *. vb) ) + +let ( -! ) = two_variable_operator + ~update_value: (fun wc vc wb vb sw -> (vc -. vb) ) + +let ( /! ) = two_variable_operator + ~update_value: (fun wc vc wb vb sw -> (vc /. vb) ) + + (** Merge two consecutive blocks *) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 727feec..1924d49 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ] BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC - E_loc_zv = 0.d0 + E_loc_zv = E_loc END_PROVIDER diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index 8dd011b..9b3758b 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -174,7 +174,7 @@ END_SHELL if (calc_E_loc_zv) then if (dabs(pdmc_weight(i_walk)*pop_weight_mult) > 1.d-6) then - E_loc_zv = (E_ref-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult) + E_loc_zv = E_loc + (E_ref-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult) else E_loc_zv = 0.d0 endif From b6b9a85cb2c7ab2a8f5951ee8b00a0125cec08bd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 3 May 2016 09:01:05 +0200 Subject: [PATCH 04/20] Added E_trial --- ezfio_config/qmc.config | 1 + ocaml/Input.ml | 57 ++++++++++++++++++++++++++ ocaml/Qmcchem_edit.ml | 38 ++++++++++------- src/PROPERTIES/properties_energy.irp.f | 2 +- src/SAMPLING/pdmc_step.irp.f | 2 +- src/ezfio_interface.irp.f | 1 + src/simulation.irp.f | 8 ++++ 7 files changed, 92 insertions(+), 17 deletions(-) diff --git a/ezfio_config/qmc.config b/ezfio_config/qmc.config index 29aaac9..3789a8b 100644 --- a/ezfio_config/qmc.config +++ b/ezfio_config/qmc.config @@ -64,6 +64,7 @@ simulation ci_threshold double precision md5_key character*(32) E_ref double precision + E_trial double precision srmc_projection_time real jastrow diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 894754b..6a1ff43 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -490,6 +490,63 @@ end +module Trial_wf_energy : sig + + type t = float + val doc : string + val read : unit -> t + val write : t -> unit + val to_float : t -> float + val of_float : float -> t + val to_string : t -> string + val of_string : string -> t + +end = struct + + type t = float + let doc = "Energy of the trial wave function (au)" + + let of_float x = + if (x > 0.) then + failwith "Reference energy should not be positive."; + if (x <= -1_000_000.) then + failwith "Reference energy is too low."; + x + + + let to_float x = x + + let read () = + let _ = + Lazy.force Qputils.ezfio_filename + in + if (not (Ezfio.has_simulation_e_trial ())) then + to_float 0. + |> Ezfio.set_simulation_e_trial; + Ezfio.get_simulation_e_trial () + |> of_float + + + let write t = + let _ = + Lazy.force Qputils.ezfio_filename + in + to_float t + |> Ezfio.set_simulation_e_trial + + + let of_string x = + Float.of_string x + |> of_float + + + let to_string x = + to_float x + |> Float.to_string + + +end + module Ref_energy : sig type t = float diff --git a/ocaml/Qmcchem_edit.ml b/ocaml/Qmcchem_edit.ml index 2f6a6e7..bcad9c7 100644 --- a/ocaml/Qmcchem_edit.ml +++ b/ocaml/Qmcchem_edit.ml @@ -24,6 +24,7 @@ type field = | Method | Sampling | Ref_energy + | Trial_wf_energy | CI_threshold | Time_step | SRMC_projection_time @@ -62,6 +63,8 @@ let get field = option_to_string Sampling.read Sampling.to_string Sampling.doc | Ref_energy -> option_to_string Ref_energy.read Ref_energy.to_string Ref_energy.doc + | Trial_wf_energy -> + option_to_string Trial_wf_energy.read Trial_wf_energy.to_string Trial_wf_energy.doc | CI_threshold -> option_to_string CI_threshold.read CI_threshold.to_string CI_threshold.doc | Time_step -> @@ -105,7 +108,7 @@ let write_input_in_ezfio ezfio_filename fields = (** Run the edit command *) -let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = +let run ~c ?f ?t ?l ?m ?e ?et ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = let interactive = ref ( if c then @@ -130,18 +133,19 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = in (); in - handle_option Input.Ref_energy.(of_float , write) e; - handle_option Input.Jastrow_type.(of_string, write) j; - handle_option Input.Block_time.(of_int , write) l; - handle_option Input.Method.(of_string, write) m; - handle_option Input.Stop_time.(of_int , write) t; - handle_option Input.Sampling.(of_string, write) s; - handle_option Input.Fitcusp_factor.(of_float , write) f; - handle_option Input.Time_step.(of_float , write) ts; - handle_option Input.Walk_num.(of_int , write) w; - handle_option Input.Walk_num_tot.(of_int , write) wt; - handle_option Input.CI_threshold.(of_float , write) n; - handle_option Input.SRMC_projection_time.(of_float , write) p; + handle_option Input.Ref_energy.(of_float , write) e; + handle_option Input.Trial_wf_energy.(of_float , write) et; + handle_option Input.Jastrow_type.(of_string, write) j; + handle_option Input.Block_time.(of_int , write) l; + handle_option Input.Method.(of_string, write) m; + handle_option Input.Stop_time.(of_int , write) t; + handle_option Input.Sampling.(of_string, write) s; + handle_option Input.Fitcusp_factor.(of_float , write) f; + handle_option Input.Time_step.(of_float , write) ts; + handle_option Input.Walk_num.(of_int , write) w; + handle_option Input.Walk_num_tot.(of_int , write) wt; + handle_option Input.CI_threshold.(of_float , write) n; + handle_option Input.SRMC_projection_time.(of_float , write) p; let fields = @@ -153,6 +157,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = Time_step ; SRMC_projection_time ; Ref_energy ; + Trial_wf_energy ; Walk_num ; Walk_num_tot ; Fitcusp_factor ; @@ -225,6 +230,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = | Walk_num_tot -> Walk_num_tot.(of_string s |> write) | CI_threshold -> CI_threshold.(of_string s |> write) | Jastrow_type -> Jastrow_type.(of_string s |> write) + | Trial_wf_energy -> Trial_wf_energy.(of_string s |> write) | Properties -> Properties.(of_string s |> write) end with @@ -281,6 +287,8 @@ let spec = ~doc:("method "^Input.Method.doc) +> flag "e" (optional float) ~doc:("energy "^Input.Ref_energy.doc) + +> flag "et" (optional float) + ~doc:("energy "^Input.Trial_wf_energy.doc) +> flag "s" (optional string) ~doc:("sampling "^Input.Sampling.doc) +> flag "ts" (optional float) @@ -307,8 +315,8 @@ let command = Edit input data ") spec - (fun c f t l m e s ts w wt n j p ezfio_file input () -> - run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_file ) + (fun c f t l m e et s ts w wt n j p ezfio_file input () -> + run ~c ?f ?t ?l ?m ?e ?et ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_file ) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 1924d49..4dd387f 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ] BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC - E_loc_zv = E_loc + E_loc_zv = E_trial END_PROVIDER diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index 9b3758b..592c1a9 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -174,7 +174,7 @@ END_SHELL if (calc_E_loc_zv) then if (dabs(pdmc_weight(i_walk)*pop_weight_mult) > 1.d-6) then - E_loc_zv = E_loc + (E_ref-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult) + E_loc_zv = E_loc + (E_trial-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult) else E_loc_zv = 0.d0 endif diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index 08c9cfe..0bb76be 100644 --- a/src/ezfio_interface.irp.f +++ b/src/ezfio_interface.irp.f @@ -45,6 +45,7 @@ data = [ \ ("simulation_http_server" , "character*(128)", "" ), ("simulation_md5_key" , "character*(32)" , "" ), ("simulation_e_ref" , "double precision" , "" ), +("simulation_e_trial" , "double precision" , "" ), ("simulation_do_run" , "logical " , "" ), ("pseudo_do_pseudo" , "logical " , "" ), diff --git a/src/simulation.irp.f b/src/simulation.irp.f index ceb66d9..d3ba571 100644 --- a/src/simulation.irp.f +++ b/src/simulation.irp.f @@ -365,3 +365,11 @@ BEGIN_PROVIDER [ character*(32), md5_key ] endif END_PROVIDER +BEGIN_PROVIDER [ double precision, E_trial ] + implicit none + BEGIN_DOC + ! Energy of the trial wave function + END_DOC + call get_simulation_e_trial(E_trial) +END_PROVIDER + From 2d269aa1fc4924da4d0a2c6af0fbc278d68c3898 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 3 May 2016 09:15:43 +0200 Subject: [PATCH 05/20] Created pdmc_weight --- src/PROPERTIES/properties_energy.irp.f | 14 +++++++++++++- src/SAMPLING/pdmc_step.irp.f | 14 ++++++-------- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 4dd387f..f4a47d1 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -156,6 +156,14 @@ BEGIN_PROVIDER [ double precision, E_kin_elec, (elec_num) ] END_PROVIDER +BEGIN_PROVIDER [ double precision, dmc_zv_weight ] + implicit none + BEGIN_DOC + ! Weight for Zero-variance in DMC + END_DOC + dmc_zv_weight = 1.d0 +END_PROVIDER + !==========================================================================! ! PROPERTIES ! @@ -253,7 +261,11 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ] BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC - E_loc_zv = E_trial + E_loc_zv = E_loc + (E_trial-E_loc) * dmc_zv_weight + + E_loc_zv_min = min(E_loc_zv,E_loc_zv_min) + E_loc_zv_max = max(E_loc_zv,E_loc_zv_max) + SOFT_TOUCH E_loc_zv_min E_loc_zv_max END_PROVIDER diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index 592c1a9..00eb9a0 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -172,14 +172,12 @@ END_SHELL psi_value_save(i_walk) = psi_value E_loc_save(i_walk) = E_loc - if (calc_E_loc_zv) then - if (dabs(pdmc_weight(i_walk)*pop_weight_mult) > 1.d-6) then - E_loc_zv = E_loc + (E_trial-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult) - else - E_loc_zv = 0.d0 - endif - TOUCH E_loc_zv - endif + 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) + else + dmc_zv_weight = 0.d0 + endif + TOUCH dmc_zv_weight BEGIN_SHELL [ /usr/bin/python ] from properties import * From bb774c319fc96b6f75a946e488d808d155eaf655 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 3 May 2016 21:10:25 +0200 Subject: [PATCH 06/20] 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 From 3b94f15bcef4a29a3c743dc841c9cd7d5dbdcf76 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 6 May 2016 21:43:42 +0200 Subject: [PATCH 07/20] E_diag --- src/PROPERTIES/properties_energy.irp.f | 27 ++++- src/PROPERTIES/properties_general.irp.f | 2 +- src/SAMPLING/pdmc_step.irp.f | 143 +++++++++++++++--------- 3 files changed, 115 insertions(+), 57 deletions(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index f4a47d1..d76d474 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -164,6 +164,14 @@ BEGIN_PROVIDER [ double precision, dmc_zv_weight ] dmc_zv_weight = 1.d0 END_PROVIDER +BEGIN_PROVIDER [ double precision, dmc_zv_weight_half ] + implicit none + BEGIN_DOC + ! Weight for Zero-variance in DMC + END_DOC + dmc_zv_weight_half = 1.d0 +END_PROVIDER + !==========================================================================! ! PROPERTIES ! @@ -256,17 +264,26 @@ BEGIN_PROVIDER [ double precision, E_loc ] END_PROVIDER -BEGIN_PROVIDER [ double precision, E_loc_zv ] +BEGIN_PROVIDER [ double precision, E_loc_zv, (3) ] implicit none BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC - E_loc_zv = E_loc + (E_trial-E_loc) * dmc_zv_weight + E_loc_zv(1) = E_loc + (E_trial-E_loc) * dmc_zv_weight_half + E_loc_zv(2) = E_loc + (E_trial-E_loc) * dmc_zv_weight + E_loc_zv(3) = dmc_zv_weight_half + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, E_loc_zv_diag ] + implicit none + BEGIN_DOC + ! Zero-variance parameter on E_loc + END_DOC + E_loc_zv_diag = E_trial - E_loc_zv_min = min(E_loc_zv,E_loc_zv_min) - E_loc_zv_max = max(E_loc_zv,E_loc_zv_max) - SOFT_TOUCH E_loc_zv_min E_loc_zv_max END_PROVIDER + diff --git a/src/PROPERTIES/properties_general.irp.f b/src/PROPERTIES/properties_general.irp.f index 1cf0bdc..f545b71 100644 --- a/src/PROPERTIES/properties_general.irp.f +++ b/src/PROPERTIES/properties_general.irp.f @@ -56,7 +56,7 @@ BEGIN_PROVIDER [ double precision, pop_weight ] 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 + pop_weight = pdmc_pop_weight_mult(1) endif pop_weight_min = min(pop_weight,pop_weight_min) pop_weight_max = max(pop_weight,pop_weight_max) diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index de64aa2..eedd919 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -110,10 +110,14 @@ END_SHELL call abrt(irp_here,'walk_num > 1') 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) + H = 0.d0 + S = 0.d0 + do while (loop) - ! Every walker makes a step - do i_walk=1,walk_num + i_walk = 1 if (.not.first_loop) then integer :: i,j,l @@ -157,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 + elec_coord(elec_num+1,3) = pdmc_weight(i_walk) * pdmc_pop_weight_mult(1) do l=1,3 do i=1,elec_num+1 elec_coord_full(i,l,i_walk) = elec_coord(i,l) @@ -172,30 +176,32 @@ 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.d-6) then - dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult) + 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)) else dmc_zv_weight = 0.d0 + dmc_zv_weight_half = 0.d0 endif - TOUCH dmc_zv_weight + TOUCH dmc_zv_weight dmc_zv_weight_half 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_weight(i_walk) - ! $X_2_pdmc_block_walk += $X_2 * pdmc_pop_weight_mult * pdmc_weight(i_walk) + ! $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) ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm - $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 3) = $X * pdmc_pop_weight_mult(1) * 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_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(1) * 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,46 +219,58 @@ for p in properties: END_SHELL - block_weight += pdmc_pop_weight_mult * pdmc_weight(i_walk) + block_weight += pdmc_pop_weight_mult(1) * pdmc_weight(i_walk) - else +! 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(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) + + 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 - pdmc_pop_weight_mult = 1.d0 - endif + 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 > 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 - ! Move to the next projection step - if (pdmc_projection > 0) then - pdmc_projection_step = mod(pdmc_projection_step,pdmc_projection)+1 - else - pdmc_projection_step = 1 - endif - - ! Eventually, recompute the weight of the population - if (pdmc_projection_step == 1) then - pdmc_pop_weight_mult = 1.d0 - do k=1,pdmc_projection - pdmc_pop_weight_mult *= pdmc_pop_weight(k) - enddo - endif - - ! Remove contribution of the old value of the weight at the new - ! 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 - sum_weight = 0.d0 - do k=1,walk_num - sum_weight += pdmc_weight(k) - enddo - pdmc_pop_weight(pdmc_projection_step) = sum_weight/dble(walk_num) - - ! Update the running population weight - pdmc_pop_weight_mult *= pdmc_pop_weight(pdmc_projection_step) call system_clock(cpu1, count_rate, count_max) if (cpu1 < cpu0) then @@ -266,12 +284,13 @@ END_SHELL cpu2 = cpu1 endif - SOFT_TOUCH elec_coord_full pdmc_pop_weight_mult + 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 @@ -288,13 +307,25 @@ for p in properties: print t.replace("$X",p[1]) 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 +! 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) & +! + 5.d0*E_trial*E_trial & +! - 8.d0*E_loc_zv_pdmc_block_walk(1)*E_trial & +! + 4.d0*E_loc_zv_pdmc_block_walk(1)*E_loc_zv_pdmc_block_walk(1) & +! ) - E_loc_zv_pdmc_block_walk(2) - E_trial) + deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp ) END_PROVIDER BEGIN_PROVIDER [ integer, pdmc_projection ] -&BEGIN_PROVIDER [ integer, pdmc_projection_step ] +&BEGIN_PROVIDER [ integer, pdmc_projection_step, (pdmc_n_diag) ] implicit none BEGIN_DOC ! Number of projection steps for PDMC @@ -303,23 +334,33 @@ END_PROVIDER pdmc_projection_time = 1. call get_simulation_srmc_projection_time(pdmc_projection_time) pdmc_projection = int( pdmc_projection_time/time_step) - pdmc_projection_step = 0 + pdmc_projection_step(:) = 0 END_PROVIDER -BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection+1) ] +BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection+1,pdmc_n_diag) ] implicit none BEGIN_DOC ! Population weight of PDMC END_DOC - pdmc_pop_weight(:) = 1.d0 + pdmc_pop_weight(:,:) = 1.d0 END_PROVIDER -BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult ] +BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (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(:) = pdmc_pop_weight(pdmc_projection,:) END_PROVIDER +BEGIN_PROVIDER [ integer, pdmc_n_diag ] + implicit none + BEGIN_DOC +! Size of the matrix to diagonalize + END_DOC + pdmc_n_diag = 2 +END_PROVIDER + + + From ca21d44409a6d8cc931e06393451686a194eed78 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 7 May 2016 00:58:45 +0200 Subject: [PATCH 08/20] 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 From 09bf6140c86ff9e0d65d157520026a3528f4638b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 9 May 2016 09:21:10 +0200 Subject: [PATCH 09/20] General diag --- src/PROPERTIES/properties_energy.irp.f | 9 ++--- src/PROPERTIES/properties_general.irp.f | 2 +- src/SAMPLING/pdmc_step.irp.f | 48 ++++++++++++++----------- 3 files changed, 33 insertions(+), 26 deletions(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index d76d474..07680e1 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -264,14 +264,15 @@ BEGIN_PROVIDER [ double precision, E_loc ] END_PROVIDER -BEGIN_PROVIDER [ double precision, E_loc_zv, (3) ] +BEGIN_PROVIDER [ double precision, E_loc_zv, ((pdmc_n_diag+1)*2) ] implicit none BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC - E_loc_zv(1) = E_loc + (E_trial-E_loc) * dmc_zv_weight_half - E_loc_zv(2) = E_loc + (E_trial-E_loc) * dmc_zv_weight - E_loc_zv(3) = dmc_zv_weight_half +! E_loc_zv(1) = E_loc + (E_trial-E_loc) * dmc_zv_weight_half +! E_loc_zv(2) = E_loc + (E_trial-E_loc) * dmc_zv_weight +! E_loc_zv(3) = dmc_zv_weight_half + E_loc_zv(:) = 0.d0 END_PROVIDER diff --git a/src/PROPERTIES/properties_general.irp.f b/src/PROPERTIES/properties_general.irp.f index f545b71..b15258c 100644 --- a/src/PROPERTIES/properties_general.irp.f +++ b/src/PROPERTIES/properties_general.irp.f @@ -56,7 +56,7 @@ BEGIN_PROVIDER [ double precision, pop_weight ] 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(1) + pop_weight = pdmc_pop_weight_mult(pdmc_n_diag) endif pop_weight_min = min(pop_weight,pop_weight_min) pop_weight_max = max(pop_weight,pop_weight_max) diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index dc18daa..86e710f 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(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) + 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 @@ -150,7 +150,7 @@ END_SHELL logical :: accepted call brownian_step(p,q,accepted,delta_x) - if ( psi_value * psi_value_save(i_walk) >= 0.d0 ) then +! 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 pdmc_weight(i_walk) = dexp(-dtime_step*thr) @@ -185,6 +185,11 @@ END_SHELL 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 = """ @@ -221,9 +226,6 @@ END_SHELL block_weight += pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - H(0,0) += E_loc - S(0,0) += 1.d0 - 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 @@ -231,14 +233,13 @@ END_SHELL S(k,l) += pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk) enddo enddo + H = H + (E_trial - E_loc) - 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 +! 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 @@ -303,17 +304,22 @@ for p in properties: print t.replace("$X",p[1]) END_SHELL + H = 0.d0 + H(0,0) = E_loc_zv_pdmc_block_walk(1) + H(1,0) = E_loc_zv_pdmc_block_walk(3) + H(0,1) = E_loc_zv_pdmc_block_walk(3) + H(1,1) = E_loc_zv_pdmc_block_walk(5) + + S = 0.d0 + S(0,0) = E_loc_zv_pdmc_block_walk(6) + S(1,0) = E_loc_zv_pdmc_block_walk(8) + S(0,1) = E_loc_zv_pdmc_block_walk(8) + S(1,1) = E_loc_zv_pdmc_block_walk(10) + S(2,2) = 1.d0 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) & -! + 5.d0*E_trial*E_trial & -! - 8.d0*E_loc_zv_pdmc_block_walk(1)*E_trial & -! + 4.d0*E_loc_zv_pdmc_block_walk(1)*E_loc_zv_pdmc_block_walk(1) & -! ) - E_loc_zv_pdmc_block_walk(2) - E_trial) + deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp ) @@ -359,7 +365,7 @@ BEGIN_PROVIDER [ integer, pdmc_n_diag ] BEGIN_DOC ! Size of the matrix to diagonalize END_DOC - pdmc_n_diag = 2 + pdmc_n_diag = 8 END_PROVIDER From 605ee8018c0cda9f953606c86edf40ba175d725c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 26 May 2016 19:58:43 +0200 Subject: [PATCH 10/20] ZV estimator --- src/PROPERTIES/properties_energy.irp.f | 8 ++--- src/SAMPLING/pdmc_step.irp.f | 41 +++++++++++++------------- 2 files changed, 25 insertions(+), 24 deletions(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 07680e1..357b21d 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -264,15 +264,15 @@ BEGIN_PROVIDER [ double precision, E_loc ] END_PROVIDER -BEGIN_PROVIDER [ double precision, E_loc_zv, ((pdmc_n_diag+1)*2) ] +!BEGIN_PROVIDER [ double precision, E_loc_zv, ((pdmc_n_diag+1)*2) ] +BEGIN_PROVIDER [ double precision, E_loc_zv ] implicit none BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC -! E_loc_zv(1) = E_loc + (E_trial-E_loc) * dmc_zv_weight_half -! E_loc_zv(2) = E_loc + (E_trial-E_loc) * dmc_zv_weight + E_loc_zv = E_loc + (E_trial-E_loc) * dmc_zv_weight ! E_loc_zv(3) = dmc_zv_weight_half - E_loc_zv(:) = 0.d0 +! E_loc_zv(:) = 0.d0 END_PROVIDER diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index 86e710f..c415ee2 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -185,10 +185,10 @@ END_SHELL 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 +! 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 * @@ -304,22 +304,23 @@ for p in properties: print t.replace("$X",p[1]) END_SHELL - H = 0.d0 - H(0,0) = E_loc_zv_pdmc_block_walk(1) - H(1,0) = E_loc_zv_pdmc_block_walk(3) - H(0,1) = E_loc_zv_pdmc_block_walk(3) - H(1,1) = E_loc_zv_pdmc_block_walk(5) - - S = 0.d0 - S(0,0) = E_loc_zv_pdmc_block_walk(6) - S(1,0) = E_loc_zv_pdmc_block_walk(8) - S(0,1) = E_loc_zv_pdmc_block_walk(8) - S(1,1) = E_loc_zv_pdmc_block_walk(10) - S(2,2) = 1.d0 - - 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) - +! 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 ) From ec1b79489303a795426c31b3c0eeb1fc1cdb1211 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 24 Jun 2016 09:14:18 +0200 Subject: [PATCH 11/20] Better time step error --- src/SAMPLING/pdmc_step.irp.f | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index c415ee2..0efb370 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -44,8 +44,7 @@ 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 :: psi_value_save(walk_num) double precision :: psi_value_save_tmp(walk_num) double precision :: pdmc_weight(walk_num) @@ -54,7 +53,6 @@ END_SHELL !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 :: pdmc_weight @@ -131,7 +129,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 @@ -142,7 +140,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 @@ -151,13 +149,19 @@ 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 - pdmc_weight(i_walk) = dexp(-dtime_step*thr) - else if ( delta < -thr ) then - pdmc_weight(i_walk) = dexp(dtime_step*thr) - else + +!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 @@ -174,7 +178,10 @@ END_SHELL enddo psi_value_save(i_walk) = psi_value - E_loc_save(i_walk) = E_loc + 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)) From 619db89adc40ae2aa90bd0620c5f202d7f2f3021 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 1 Jul 2016 21:56:15 +0200 Subject: [PATCH 12/20] Conflicts: src/PROPERTIES/properties_general.irp.f --- README.md | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index 24fb732..ba76183 100644 --- a/README.md +++ b/README.md @@ -90,12 +90,12 @@ $ ninja Example of a QMC=Chem calculation --------------------------------- -Calculation with the [quantum package](http://github.com/LCPQ/quantum_package) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +1.Calculation with the [quantum package](http://github.com/LCPQ/quantum_package) -1) Create the ``xyz`` file containing the nuclear coordinates of the system -``` +Create the `xyz` file containing the nuclear coordinates of the system + +```bash $ cat > h2o.xyz << EOF 3 Water molecule @@ -105,31 +105,32 @@ H -0.239987 0.926627 0. EOF ``` -2) Choose a suitable basis set and create the [EZFIO database](https://github.com/LCPQ/ezfio) +Choose a suitable basis set and create the [EZFIO database](https://github.com/LCPQ/ezfio) ```bash $ qp_create_ezfio_from_xyz -b cc-pvdz h2o.xyz -o h2o ``` -3) Run the SCF calculation +Run the SCF calculation ```bash $ qp_run SCF h2o ``` -4) Run the CIPSI calculation +Run the CIPSI calculation ```bash $ qp_run full_ci h2o ``` -5) Transform the input for use in QMC=Chem +Transform the input for use in QMC=Chem ```bash $ qp_run save_for_qmcchem h2o ``` -FN-DMC calculation with QMC=Chem -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +2.FN-DMC calculation with QMC=Chem + Before using QMC=Chem, you need to load the environment variables: From d71c80afcae88bedb9cf3a6733cb472dc8d54572 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Jul 2016 16:04:49 +0200 Subject: [PATCH 13/20] Added D and F spherical harmonics in pseudo --- src/pseudo.irp.f | 57 +++++++++++++++++++++++++++++++++--------------- 1 file changed, 40 insertions(+), 17 deletions(-) diff --git a/src/pseudo.irp.f b/src/pseudo.irp.f index dd01e01..156c25a 100644 --- a/src/pseudo.irp.f +++ b/src/pseudo.irp.f @@ -333,31 +333,54 @@ real function ylm(l,m,x,y,z,r_inv) continue case(1) + ylm = ylm * sq3 * r_inv select case (m) case (-1) - ylm = ylm * sq3 * y * r_inv + ylm = ylm * y case (0) - ylm = ylm * sq3 * z * r_inv + ylm = ylm * z case (1) - ylm = ylm * sq3 * x * r_inv + ylm = ylm * x + end select + + case(2) + ylm = ylm * r_inv * r_inv * sqrt(5.) + select case (m) + case(-2) + ylm = ylm * sqrt(3.) * x * y + case(-1) + ylm = ylm * sqrt(3.) * y * z + case(0) + ylm = ylm * 0.5 * (2.*z*z - x*x - y*y) + case(1) + ylm = ylm * sqrt(3.) * z * x + case(2) + ylm = ylm * 0.5 * sqrt(3.) * (x*x - y*y) + end select + + case(3) + ylm = ylm * r_inv * r_inv * r_inv * sqrt(7.) + select case (m) + case(-3) + ylm = ylm * 0.25 * sqrt(10.) * (3.*x*x - y*y) + case(-2) + ylm = ylm * sqrt(15.) * x*y*z + case(-1) + ylm = ylm * 0.25*sqrt(6.) * y * (4.*z*z - x*x -y*y) + case(0) + ylm = ylm * 0.5 * z * (2.*z*z - 3.*(x*x + y*y)) + case(1) + ylm = ylm * 0.25*sqrt(6.) * x * (4.*z*z - x*x -y*y) + case(2) + ylm = ylm * 0.5*sqrt(15.) * z * (x*x -y*y) + case(3) + ylm = ylm * 0.25*sqrt(10.) * x * (x*x - 3. * y*y) end select -! case(2) -! select case (m) -! case(-2) -! ylm = ylm * sqrt(15.) * x * y * r_inv * r_inv -! case(-1) -! ylm = ylm * sqrt(15.) * y * z * r_inv * r_inv -! case(0) -! ylm = 0.5 * ylm * sqrt(15.) * (2.*z*z - x*x - y*y) * r_inv * r_inv -! case(1) -! ylm = ylm * sqrt(15.) * z * x * r_inv * r_inv -! case(2) -! ylm = 0.5 * ylm * sqrt(15.) * (x*x - y*y) * r_inv * r_inv -! end select case default - stop 'problem in Ylm of pseudo' + print *, 'l=', l + stop 'problem in Ylm of pseudo : Ylm not implemented (pseudo.irp.f)' end select From 8733aa742ec6afb7142c9d5c562023def38d7a53 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jul 2016 13:51:47 +0200 Subject: [PATCH 14/20] Removed check of Jastrow and pseudo --- ocaml/Input.ml | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 01e8f97..ba7a001 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -758,13 +758,6 @@ end = struct let _ = Lazy.force Qputils.ezfio_filename in - let () = - match (Pseudo.read () |> Pseudo.to_bool, t) with - | (false, _) - | (true , None) -> () - | _ -> failwith "Jastrow and Pseudopotentials are incompatible for now" - in - to_string t |> Ezfio.set_jastrow_jast_type @@ -843,8 +836,6 @@ let validate () = Sampling.read () and ts = Time_step.read () - and jast_type = - Jastrow_type.read () and do_pseudo = Pseudo.read () in @@ -900,14 +891,6 @@ let validate () = | _ -> () in - (* Pseudo and Jastrow are incompatible *) - let () = - match (Pseudo.to_bool do_pseudo, jast_type) with - | (true, Jastrow_type.Core ) - | (true, Jastrow_type.Simple) -> failwith "Jastrow and Pseudopotentials are incompatible" - | _ -> () - in - (* Fitcusp is incompatible with pseudo *) let () = let f = From b577984de6a0b5a962100269591d9df636486b5b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 3 Oct 2016 00:34:42 +0200 Subject: [PATCH 15/20] Updated ZMQ --- install/build.ninja | 7 +++---- install/scripts/install_zmq.sh | 10 +++++----- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/install/build.ninja b/install/build.ninja index 351db7d..34ab4f7 100644 --- a/install/build.ninja +++ b/install/build.ninja @@ -7,15 +7,14 @@ URL_OPAM ="https://raw.github.com/ocaml/opam/master/shell/opam_installer.sh" URL_IRPF90="https://github.com/scemama/irpf90/archive/v1.6.7.tar.gz" URL_EZFIO ="https://github.com/scemama/EZFIO/archive/v1.3.1.tar.gz" -URL_ZMQ ="http://download.zeromq.org/zeromq-4.0.7.tar.gz" -#URL_ZMQ ="http://download.zeromq.org/zeromq-4.1.3.tar.gz" -URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/v4.1.3.tar.gz" +URL_ZMQ ="http://download.zeromq.org/zeromq-4.1.4.tar.gz" +URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/4.1.4.tar.gz" # Rules ####### rule download - command = [[ -e ${out} ]] || (wget --no-check-certificate ${url} -O ${out}.tmp -o /dev/null && mv ${out}.tmp ${out}) + command = [ -e ${out} ] || (wget --no-check-certificate ${url} -O ${out}.tmp -o /dev/null && mv ${out}.tmp ${out}) description = Downloading ${descr} rule install diff --git a/install/scripts/install_zmq.sh b/install/scripts/install_zmq.sh index 80da3a3..6c7609c 100755 --- a/install/scripts/install_zmq.sh +++ b/install/scripts/install_zmq.sh @@ -3,7 +3,7 @@ TARGET=zmq function _install() { - LIBVERSION=4 + LIBVERSION=5 cd .. ; QMCCHEM_PATH="$PWD" ; cd - set +u export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./ @@ -14,10 +14,10 @@ function _install() make -j 8 cd - rm -f -- "${QMCCHEM_PATH}"/lib/libzmq.{a,so,so.$LIBVERSION} -# cp "${BUILD}"/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ -# cp "${BUILD}"/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION - cp "${BUILD}"/src/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ - cp "${BUILD}"/src/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION + cp "${BUILD}"/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ + cp "${BUILD}"/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION +# cp "${BUILD}"/src/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ +# cp "${BUILD}"/src/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION cp "${BUILD}"/include/{zmq,zmq_utils}.h "${QMCCHEM_PATH}"/lib/ cd "${QMCCHEM_PATH}"/lib ln libzmq.so.$LIBVERSION libzmq.so || cp libzmq.so.$LIBVERSION libzmq.so From dcf19db6d625743dfc7316b537b4c1ddc880a634 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 28 Nov 2016 15:08:55 +0100 Subject: [PATCH 16/20] Commented prefetch instructions --- src/det.irp.f | 7 +++++++ src/mo.irp.f | 28 ++++++++++++++-------------- 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/src/det.irp.f b/src/det.irp.f index 0982d07..85494f2 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -1114,6 +1114,9 @@ end endif !DIR$ FORCEINLINE call bitstring_to_list ( psi_det_alpha(1,det_i), mo_list_alpha_curr, l, N_int ) + if (l /= elec_alpha_num) then + stop 'error in number of alpha electrons' + endif END_PROVIDER @@ -1132,8 +1135,12 @@ END_PROVIDER else mo_list_beta_prev = 0 endif + !DIR$ FORCEINLINE call bitstring_to_list ( psi_det_beta(1,det_j), mo_list_beta_curr, l, N_int ) + if (l /= elec_beta_num) then + stop 'error in number of beta electrons' + endif END_PROVIDER BEGIN_PROVIDER [ double precision, det_alpha_value_curr ] diff --git a/src/mo.irp.f b/src/mo.irp.f index eb300af..4b9bba7 100644 --- a/src/mo.irp.f +++ b/src/mo.irp.f @@ -701,13 +701,13 @@ subroutine sparse_full_mv(A,LDA, & ! LDC and LDA have to be factors of simd_sp - IRP_IF NO_PREFETCH - IRP_ELSE - call MM_PREFETCH (A(j,indices(1)),3) - call MM_PREFETCH (A(j,indices(2)),3) - call MM_PREFETCH (A(j,indices(3)),3) - call MM_PREFETCH (A(j,indices(4)),3) - IRP_ENDIF +! IRP_IF NO_PREFETCH +! IRP_ELSE +! call MM_PREFETCH (A(1,indices(1)),3) +! call MM_PREFETCH (A(1,indices(2)),3) +! call MM_PREFETCH (A(1,indices(3)),3) +! call MM_PREFETCH (A(1,indices(4)),3) +! IRP_ENDIF !DIR$ SIMD do j=1,LDC @@ -757,13 +757,13 @@ subroutine sparse_full_mv(A,LDA, & !DIR$ VECTOR ALIGNED !DIR$ SIMD FIRSTPRIVATE(d11,d21,d31,d41) do j=1,$IRP_ALIGN/4 - IRP_IF NO_PREFETCH - IRP_ELSE - call MM_PREFETCH (A(j+k,indices(kao+4)),3) - call MM_PREFETCH (A(j+k,indices(kao+5)),3) - call MM_PREFETCH (A(j+k,indices(kao+6)),3) - call MM_PREFETCH (A(j+k,indices(kao+7)),3) - IRP_ENDIF +! IRP_IF NO_PREFETCH +! IRP_ELSE +! call MM_PREFETCH (A(j+k,indices(kao+4)),3) +! call MM_PREFETCH (A(j+k,indices(kao+5)),3) +! call MM_PREFETCH (A(j+k,indices(kao+6)),3) +! call MM_PREFETCH (A(j+k,indices(kao+7)),3) +! IRP_ENDIF C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11 + A(j+k,k_vec(2))*d21& + A(j+k,k_vec(3))*d31 + A(j+k,k_vec(4))*d41 enddo From 4c42654401513d709721f90ea37c144f52007e9b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 28 Nov 2016 19:45:09 +0100 Subject: [PATCH 17/20] Cleaned PDMC --- src/PROPERTIES/properties_energy.irp.f | 4 +++- src/SAMPLING/pdmc_step.irp.f | 20 ++++++++++---------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 357b21d..7f42053 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -270,7 +270,9 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ] BEGIN_DOC ! Zero-variance parameter on E_loc END_DOC - E_loc_zv = E_loc + (E_trial-E_loc) * dmc_zv_weight + E_loc_zv = E_loc + E_loc_zv += (E_trial-E_loc) * dmc_zv_weight +! E_loc_zv += - time_step*(E_trial**2 + 1.44341217940434 - E_loc**2)*dmc_zv_weight ! E_loc_zv(3) = dmc_zv_weight_half ! E_loc_zv(:) = 0.d0 diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index 0efb370..9c9b4a8 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -109,9 +109,9 @@ END_SHELL 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 +! 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) @@ -234,13 +234,13 @@ 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) +! 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 From bb4b75e86b2fd1130d8dc069b02a8a84d24a7983 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 28 Dec 2016 16:57:53 +0100 Subject: [PATCH 18/20] 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 From a1192beb8b8672bd84c17ed4d376484eb0194cda Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 29 Dec 2016 01:43:34 +0100 Subject: [PATCH 19/20] Promela model --- promela/qmcchem.pml | 271 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 271 insertions(+) create mode 100644 promela/qmcchem.pml diff --git a/promela/qmcchem.pml b/promela/qmcchem.pml new file mode 100644 index 0000000..a0c001b --- /dev/null +++ b/promela/qmcchem.pml @@ -0,0 +1,271 @@ +#define NPROC 2 +#define BUFSIZE 4 +#define not_here false + +mtype = { NONE, TERMINATE, OK, TEST, ERROR, PROPERTY, WALKERS, EZFIO, GETWALKERS, REGISTER, + EZFIO_REPLY, UNREGISTER, STOPPING, STOPPED, QUEUED, RUNNING }; + +typedef message_req { + mtype m = NONE; + byte value = 0; + chan reply = [BUFSIZE] of { mtype }; +} + +typedef message_pull { + mtype m = NONE; + byte value = 0; +} + +chan dataserver_pull = [NPROC] of { message_pull }; +chan dataserver_req = [NPROC] of { message_req }; + +byte dataserver_status_pub; +bit http_address = 0; +bit killall_qmc = 0; +bit killall_dataserver = 0; +byte dataserver_status = QUEUED; +byte dataserver_status_n_connected = 0; + + + + +/* qmcchem process */ +active proctype qmcchem() { + byte reply = NONE; + byte dataserver_pid; + byte i,j; + message_req msg; + + dataserver_pid = run dataserver(); + + /* Wait until ZMQ socket is open */ + (http_address == 1); + do + :: (reply == OK) -> break + :: (reply == NONE) -> + msg.m = TEST; + dataserver_req ! msg; + msg.reply ? reply ; + assert (reply == OK || reply == NONE) + od; + printf("Dataserver is ready.\n"); + + /* Start the QMC processes */ + + printf("qmcchem: Starting qmc processes.\n"); + atomic { + i=0; + do + :: (i < NPROC) -> + run qmc(); i++ + :: else -> break + od; + } + printf("qmcchem: qmc processes started.\n"); + +} + + + + + + + +/* dataserver process */ +proctype dataserver() { + byte reply = 0; + byte request = 0; + byte cont = 0; + byte reply_pid = 0; + message_req msg; + + /* Simulate initialization */ + http_address = 1; + dataserver_req ? msg; + msg.reply ! NONE ; + + /* Status thread */ + run dataserver_status_thread(); + run dataserver_main_thread(); +} + +#define delay 5 +#define stop_time 100 + + +proctype dataserver_status_thread() { + byte count=0; + byte n_connected = 0; + byte time=0; + + dataserver_status_pub = dataserver_status; + do + :: (dataserver_status == STOPPED) -> break + :: else -> + time = (time < stop_time -> time+1 : time); + count++; + if + :: (count != delay) -> skip + :: else -> + count = 0; + if + :: (dataserver_status == RUNNING && + n_connected == dataserver_status_n_connected && + time >= stop_time) -> + dataserver_status = STOPPING; + printf("Stop time reached : STOPPING\n") + + :: (dataserver_status == STOPPING && + n_connected != dataserver_status_n_connected && + dataserver_status_n_connected == 0) -> + dataserver_status = STOPPED; + printf("No more connected clients : STOPPED\n") + + :: (n_connected != dataserver_status_n_connected && + dataserver_status_n_connected > 0) -> + n_connected = dataserver_status_n_connected; + + :: else -> skip + fi + fi + dataserver_status_pub = dataserver_status; + od + printf ("End of dataserver_status_thread\n"); +} + + +proctype dataserver_main_thread() { + byte time = 0; + mtype reply; + dataserver_status = QUEUED; + message_req msg; + message_pull pull; + + /* Inform main process that the qmc processes can start (blocking recv) */ + dataserver_req ? msg; + assert (msg.m == TEST); + msg.reply ! OK; + + do + :: (dataserver_status == STOPPED && (!dataserver_pull ?[pull]) && (!dataserver_req ?[msg])) -> break + :: else -> + do + :: (dataserver_pull ?[pull]) -> + dataserver_pull ? pull + printf("pull: "); printm(pull.m); printf("\n"); + if + :: (pull.m == ERROR) -> skip; + :: (pull.m == WALKERS) -> skip + :: (pull.m == PROPERTY) -> skip; + fi + :: else -> break + od + if + :: (dataserver_req ?[msg]) -> + dataserver_req ? msg; + printf("req : "); printm(msg.m); printf("\n"); + if + :: (msg.m == TEST) -> reply = OK + :: (msg.m == EZFIO) -> reply = EZFIO_REPLY + :: (msg.m == GETWALKERS) -> reply = WALKERS + :: (msg.m == REGISTER && dataserver_status == QUEUED ) -> + dataserver_status_n_connected++; + dataserver_status = RUNNING; + reply = OK; + printf("Status changed to RUNNING\n") + :: (msg.m == REGISTER && dataserver_status == RUNNING ) -> + dataserver_status_n_connected++; + reply = OK + :: (msg.m == REGISTER && + (dataserver_status == STOPPED || dataserver_status == STOPPING) ) -> + dataserver_status_n_connected++; reply = ERROR; + printf("dataserver_req: register failed \n") + :: (msg.m == UNREGISTER) -> + dataserver_status_n_connected--; + reply = OK; + if + :: (dataserver_status_n_connected == 0) -> + dataserver_status = STOPPED + printf("Status changed to STOPPED\n") + :: else -> skip + fi + :: else -> skip + fi; + msg.reply ! reply + :: else -> skip + fi + od +} + +/* qmc processes */ +proctype qmc() { + byte status; + mtype reply; + message_req msg; + message_pull pull; + + /* Init */ + status = dataserver_status_pub; + + msg.m = REGISTER; + dataserver_req ! msg; +end: msg.reply ? reply; + if + :: (reply == ERROR) -> goto exit; + :: else -> assert (reply == OK); + fi; + + msg.m = EZFIO; + dataserver_req ! msg; + msg.reply ? reply; + if + :: (reply == ERROR) -> goto exit; + :: else -> assert (reply == EZFIO_REPLY); + fi; + + msg.m = GETWALKERS; + dataserver_req ! msg; + msg.reply ? reply; + if + :: (reply == ERROR) -> goto exit; + :: else -> assert (reply == WALKERS); + fi; + + + + /* Equilibration */ + (dataserver_status_pub == RUNNING); + + msg.m = EZFIO; + dataserver_req ! msg; + msg.reply ? reply; + if + :: (reply == ERROR) -> goto exit; + :: else -> assert (reply == EZFIO_REPLY); + fi; + + status = dataserver_status_pub; + + /* Cycles */ + do + :: (status != RUNNING) -> break + :: else -> + pull.m = PROPERTY; pull.value = 0; + dataserver_pull ! pull; + pull.m = PROPERTY; pull.value =1 ; + dataserver_pull ! pull; + pull.m = WALKERS; + dataserver_pull ! pull; + status = dataserver_status_pub; + od; + + /* Termination */ + msg.m = UNREGISTER; + dataserver_req ! msg; + msg.reply ? reply; + assert (reply == OK); + +exit: skip +} + + From c59606bb36863dfa59309081f6b71bac943f95b6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 29 Dec 2016 01:48:59 +0100 Subject: [PATCH 20/20] Removed useless property --- src/PROPERTIES/properties_energy.irp.f | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 7f42053..12ebb73 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -278,15 +278,6 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ] END_PROVIDER -BEGIN_PROVIDER [ double precision, E_loc_zv_diag ] - implicit none - BEGIN_DOC - ! Zero-variance parameter on E_loc - END_DOC - E_loc_zv_diag = E_trial - -END_PROVIDER -