diff --git a/ezfio_config/qmc.config b/ezfio_config/qmc.config index 73812d1..d41db6d 100644 --- a/ezfio_config/qmc.config +++ b/ezfio_config/qmc.config @@ -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 diff --git a/install/build.ninja b/install/build.ninja index 758982d..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.4.tar.gz" -URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/v4.1.3.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 c3b99ea..6c7609c 100755 --- a/install/scripts/install_zmq.sh +++ b/install/scripts/install_zmq.sh @@ -3,7 +3,6 @@ TARGET=zmq function _install() { -# LIBVERSION=4 LIBVERSION=5 cd .. ; QMCCHEM_PATH="$PWD" ; cd - set +u diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 01e8f97..70b321b 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" @@ -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 @@ -758,13 +817,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 +895,6 @@ let validate () = Sampling.read () and ts = Time_step.read () - and jast_type = - Jastrow_type.read () and do_pseudo = Pseudo.read () in @@ -852,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 @@ -879,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) ) | _ -> () @@ -894,20 +939,13 @@ 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" | _ -> () 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 = 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/ocaml/Qmcchem_result.ml b/ocaml/Qmcchem_result.ml index 23e0660..c1737d6 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 ; +(* + 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 diff --git a/ocaml/Random_variable.ml b/ocaml/Random_variable.ml index 01fc04a..c1ae8ad 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,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) *) 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 +} + + diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index d21f973..12ebb73 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -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,20 @@ 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_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 + +END_PROVIDER + + diff --git a/src/PROPERTIES/properties_general.irp.f b/src/PROPERTIES/properties_general.irp.f index b6c0ae9..c5342e6 100644 --- a/src/PROPERTIES/properties_general.irp.f +++ b/src/PROPERTIES/properties_general.irp.f @@ -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 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 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..b34692f 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -209,6 +209,7 @@ END_SHELL E_loc_save(1,i_walk) = E_loc endif + BEGIN_SHELL [ /usr/bin/python ] from properties import * t = """ @@ -319,11 +320,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 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/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index bf537a8..515df24 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/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 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 diff --git a/src/simulation.irp.f b/src/simulation.irp.f index 0995e2d..d3ba571 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)) @@ -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 + 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 ', & + ' '/)