mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-11-07 06:33:38 +01:00
Merge branch 'feature/zveloc' into develop
This commit is contained in:
commit
dc187bea73
@ -61,6 +61,7 @@ simulation
|
||||
ci_threshold double precision
|
||||
md5_key character*(32)
|
||||
E_ref double precision
|
||||
E_trial double precision
|
||||
srmc_projection_time real
|
||||
|
||||
jastrow
|
||||
|
||||
@ -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"
|
||||
|
||||
|
||||
@ -488,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
|
||||
@ -843,25 +902,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
|
||||
|
||||
@ -870,6 +923,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) )
|
||||
| _ -> ()
|
||||
@ -885,6 +939,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"
|
||||
|
||||
@ -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 )
|
||||
|
||||
|
||||
|
||||
|
||||
@ -145,6 +145,15 @@ let display_summary ~range =
|
||||
in
|
||||
List.iter properties ~f:print_property ;
|
||||
|
||||
(*
|
||||
let open Random_variable in
|
||||
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);
|
||||
*)
|
||||
|
||||
let cpu =
|
||||
Random_variable.of_raw_data ~range Property.Cpu
|
||||
|
||||
@ -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,59 @@ let merge_per_compute_node_and_block_id =
|
||||
(Block_id.to_int block.Block.block_id) )
|
||||
|
||||
|
||||
(** 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_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 *)
|
||||
@ -467,7 +523,7 @@ let compress =
|
||||
((Block_id.to_int block_id)+1)/2
|
||||
|> Block_id.of_int )
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
(** Last value on each compute node (for wall_time) *)
|
||||
|
||||
@ -156,6 +156,22 @@ 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
|
||||
|
||||
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 !
|
||||
@ -248,5 +264,27 @@ 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 ]
|
||||
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(3) = dmc_zv_weight_half
|
||||
! E_loc_zv(:) = 0.d0
|
||||
|
||||
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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -47,6 +47,22 @@ BEGIN_PROVIDER [ double precision, wf_extension ]
|
||||
SOFT_TOUCH wf_extension_min wf_extension_max
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, pop_weight ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Weight of the SRMC population
|
||||
END_DOC
|
||||
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(pdmc_n_diag)
|
||||
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
|
||||
|
||||
BEGIN_PROVIDER [ double precision, drift_mod, (size_drift_mod) ]
|
||||
implicit none
|
||||
|
||||
@ -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
|
||||
|
||||
380
src/SAMPLING/pdmc_step.irp.f
Normal file
380
src/SAMPLING/pdmc_step.irp.f
Normal file
@ -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(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 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
|
||||
|
||||
|
||||
|
||||
@ -209,6 +209,7 @@ END_SHELL
|
||||
E_loc_save(1,i_walk) = E_loc
|
||||
endif
|
||||
|
||||
|
||||
BEGIN_SHELL [ /usr/bin/python ]
|
||||
from properties import *
|
||||
t = """
|
||||
|
||||
@ -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 " , "" ),
|
||||
|
||||
|
||||
@ -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))
|
||||
@ -363,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
|
||||
|
||||
|
||||
49
src/types.F
49
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 ', &
|
||||
' '/)
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user