From 591f0306eaf2e69c9dba376608bd77c851354aec Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 5 Feb 2016 00:41:10 +0100 Subject: [PATCH 1/8] Deterministic ZMQ ports --- ocaml/Qmcchem_dataserver.ml | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/ocaml/Qmcchem_dataserver.ml b/ocaml/Qmcchem_dataserver.ml index e424676..8efeb33 100644 --- a/ocaml/Qmcchem_dataserver.ml +++ b/ocaml/Qmcchem_dataserver.ml @@ -71,8 +71,6 @@ let run ?(daemon=true) ezfio_filename = - (* - (** Checks if the port is already open (not working properly yet) *) let check_port n = let adress_prefix = "tcp://*:" @@ -87,7 +85,9 @@ let run ?(daemon=true) ezfio_filename = in let result = try - (ZMQ.Socket.bind socket address; accu ); + ZMQ.Socket.bind socket address; + ZMQ.Socket.unbind socket address; + accu; with | _ -> false; in @@ -100,22 +100,18 @@ let run ?(daemon=true) ezfio_filename = else `Unavailable in - *) (** Random port number between 49152 and 65535 *) let port = let newport = - (* ref (49152 + (Random.int 16383)) *) - ref ( 1024 + (Random.int (49151-1024))) + ref 10000 in - (* while ((check_port !newport) = `Unavailable) do - newport := 49152 + (Random.int 16383) + newport := !newport + 100 done; - *) !newport in From 453a29d60724f77b1b679deb6f469caa2853bd9c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Mar 2016 13:39:06 +0100 Subject: [PATCH 2/8] Try except to clean tmpdir --- ocaml/Qmcchem_forwarder.ml | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/ocaml/Qmcchem_forwarder.ml b/ocaml/Qmcchem_forwarder.ml index deaf50a..4a421f3 100644 --- a/ocaml/Qmcchem_forwarder.ml +++ b/ocaml/Qmcchem_forwarder.ml @@ -465,6 +465,13 @@ let run ezfio_filename dataserver = end; (* Wait for the qmc process to complete *) - ignore (Watchdog.join ()); - terminate () + try + ignore (Watchdog.join ()); + terminate () + with + | error -> + begin + terminate (); + raise error + end From b66cddbe4d955d084e6774e76839e6b8dda041d4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Mar 2016 13:57:33 +0100 Subject: [PATCH 3/8] Cleaning after Ctrl-C --- ocaml/Qmcchem_dataserver.ml | 3 ++- ocaml/Qmcchem_forwarder.ml | 6 +++--- ocaml/Qmcchem_run.ml | 8 ++++---- ocaml/Watchdog.ml | 2 +- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/ocaml/Qmcchem_dataserver.ml b/ocaml/Qmcchem_dataserver.ml index afd4a8a..f63c21b 100644 --- a/ocaml/Qmcchem_dataserver.ml +++ b/ocaml/Qmcchem_dataserver.ml @@ -833,12 +833,13 @@ let run ?(daemon=true) ezfio_filename = (* Handle signals *) let handler s = + Printf.printf "Dataserver received the %s signal... killing\n%!" (Signal.to_string s); Watchdog.kill (); in List.iter [ + Signal.int ; Signal.term ; Signal.quit ; - Signal.int ] ~f:(fun x -> Signal.Expert.handle x handler) ; diff --git a/ocaml/Qmcchem_forwarder.ml b/ocaml/Qmcchem_forwarder.ml index 4a421f3..041b38f 100644 --- a/ocaml/Qmcchem_forwarder.ml +++ b/ocaml/Qmcchem_forwarder.ml @@ -71,6 +71,7 @@ let run ezfio_filename dataserver = let terminate () = (* Clean up the temp directory *) + Watchdog.kill (); Unix.chdir Qmcchem_config.dev_shm; let command = Printf.sprintf "rm -rf -- \"%s\" " tmpdir @@ -96,14 +97,13 @@ let run ezfio_filename dataserver = (* Signal handler to Kill properly all the processes *) let handler s = - Printf.printf "Forwarder received the %s signal... killing\n" (Signal.to_string s); + Printf.printf "Forwarder received the %s signal... killing\n%!" (Signal.to_string s); terminate (); - Watchdog.kill (); in List.iter [ + Signal.int ; Signal.term ; Signal.quit ; - Signal.int ] ~f:(fun x -> Signal.Expert.handle x handler) ; diff --git a/ocaml/Qmcchem_run.ml b/ocaml/Qmcchem_run.ml index f878987..03d9bc1 100644 --- a/ocaml/Qmcchem_run.ml +++ b/ocaml/Qmcchem_run.ml @@ -7,8 +7,8 @@ let full_run ?(start_dataserver=true) ezfio_filename = and scheduler = Scheduler.find () in - Printf.printf "Scheduler : %s\n" (Scheduler.to_string scheduler); - Printf.printf "Launcher : %s\n" (Launcher.to_string launcher ); + Printf.printf "Scheduler : %s\n%!" (Scheduler.to_string scheduler); + Printf.printf "Launcher : %s\n%!" (Launcher.to_string launcher ); (* Create the node file *) @@ -147,13 +147,13 @@ let run a d ?q ?s ezfio_filename = (* Signal handler to Kill properly all the processes *) let handler s = - Printf.printf "Received the %s signal... killing\n" (Signal.to_string s); + Printf.printf "QMC=Chem received the %s signal... killing\n%!" (Signal.to_string s); Watchdog.kill (); in List.iter [ + Signal.int ; Signal.term ; Signal.quit ; - Signal.int ] ~f:(fun x -> Signal.Expert.handle x handler) ; diff --git a/ocaml/Watchdog.ml b/ocaml/Watchdog.ml index a9d2e3d..3a6a0aa 100644 --- a/ocaml/Watchdog.ml +++ b/ocaml/Watchdog.ml @@ -8,7 +8,7 @@ let _threads = ref [] ;; let kill () = let kill pid = Signal.send_i Signal.int (`Pid pid); - Printf.printf "Killed %d\n" (Pid.to_int pid) + Printf.printf "Killed %d\n%!" (Pid.to_int pid) in List.iter ~f:kill (!_list); exit 1 From 58b58acf459c2dcaff5e5ea5d0e5930ab726c40f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 5 Mar 2016 00:25:39 +0100 Subject: [PATCH 4/8] Storing PID in tmpdir --- ocaml/Qmcchem_forwarder.ml | 57 ++++++++++++++++++++++++++++---------- 1 file changed, 43 insertions(+), 14 deletions(-) diff --git a/ocaml/Qmcchem_forwarder.ml b/ocaml/Qmcchem_forwarder.ml index 041b38f..09caa49 100644 --- a/ocaml/Qmcchem_forwarder.ml +++ b/ocaml/Qmcchem_forwarder.ml @@ -1,14 +1,20 @@ open Core.Std;; let bind_socket ~socket_type ~socket ~address = - try - ZMQ.Socket.bind socket address - with - | Unix.Unix_error (_, message, f) -> - failwith @@ Printf.sprintf - "\n%s\nUnable to bind the forwarder's %s socket :\n %s\n%s" - f socket_type address message - | other_exception -> raise other_exception + let rec loop = function + | 0 -> failwith @@ Printf.sprintf + "Unable to bind the forwarder's %s socket : %s\n" + socket_type address + | -1 -> () + | i -> + try + ZMQ.Socket.bind socket address; + loop (-1) + with + | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) ) + | other_exception -> raise other_exception + in loop 10 + let run ezfio_filename dataserver = @@ -45,18 +51,41 @@ let run ezfio_filename dataserver = *) let () = try - Unix.mkdir tmpdir + Unix.mkdir tmpdir; + Unix.chdir tmpdir with | Unix.Unix_error _ -> - (* TODO : wait until the forwarder has started *) begin Unix.chdir tmpdir; - ignore @@ Unix.exec ~prog ~args () + Time.pause @@ Time.Span.of_float 0.1; + match (Sys.file_exists "PID") with + | `No + | `Unknown -> () + | `Yes -> + let pid = + In_channel.with_file "PID" ~f:(fun ic -> + match (In_channel.input_line ic) with + | Some x -> x + | None -> "-1" ) + |> Int.of_string + in + match pid with + | -1 -> () + | pid -> + begin + match Signal.send (Signal.of_system_int 0) (`Pid (Pid.of_int pid)) with + | `No_such_process -> () + | _ -> ignore @@ Unix.exec ~prog ~args () + end end in - Unix.chdir tmpdir; (* Now, only one forwarder will execute the following code *) + Out_channel.with_file "PID" ~f:(fun oc -> + Unix.getpid () + |> Pid.to_int + |> Printf.sprintf "%d\n" + |> Out_channel.output_string oc); (* Fork a qmc *) ignore @@ @@ -71,7 +100,6 @@ let run ezfio_filename dataserver = let terminate () = (* Clean up the temp directory *) - Watchdog.kill (); Unix.chdir Qmcchem_config.dev_shm; let command = Printf.sprintf "rm -rf -- \"%s\" " tmpdir @@ -91,7 +119,8 @@ let run ezfio_filename dataserver = with | _ -> () ; - done + done; + Watchdog.kill () in From 1df8e21ee80f62289d684f6fe99ba220b782f85c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 17 Mar 2016 15:27:31 +0100 Subject: [PATCH 5/8] Random port --- ocaml/Qmcchem_dataserver.ml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ocaml/Qmcchem_dataserver.ml b/ocaml/Qmcchem_dataserver.ml index f63c21b..1f16c2b 100644 --- a/ocaml/Qmcchem_dataserver.ml +++ b/ocaml/Qmcchem_dataserver.ml @@ -106,11 +106,11 @@ let run ?(daemon=true) ezfio_filename = (** Random port number between 49152 and 65535 *) let port = let newport = - ref 10000 + ref ( 1024 + (Random.int (49151-1024))) in while ((check_port !newport) = `Unavailable) do - newport := !newport + 100 + newport := 1024 + (Random.int (49151-1024)) done; !newport in From caf22663b540722444063e135602429c1658f52f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 17 Mar 2016 15:32:17 +0100 Subject: [PATCH 6/8] Changed sexplib.syntax to pa_sexp_conv --- install/scripts/install_ocaml.sh | 2 +- ocaml/build.ninja | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh index 7feeb24..fc31495 100755 --- a/install/scripts/install_ocaml.sh +++ b/install/scripts/install_ocaml.sh @@ -4,7 +4,7 @@ set -u set -e cd .. ; QMCCHEM_PATH="$PWD" ; cd - -PACKAGES="core cryptokit ocamlfind sexplib" +PACKAGES="core cryptokit ocamlfind sexplib pa_sexp_conv" declare -i i i=$(gcc -dumpversion | cut -d '.' -f 2) diff --git a/ocaml/build.ninja b/ocaml/build.ninja index db12102..23814ba 100644 --- a/ocaml/build.ninja +++ b/ocaml/build.ninja @@ -1,7 +1,7 @@ MAIN=qmcchem # Main program to build -PACKAGES=-package core,sexplib.syntax,cryptokit,str,ZMQ +PACKAGES=-package core,pa_sexp_conv,cryptokit,str,ZMQ # Required opam packages, for example: # PACKAGES=-package core,sexplib.syntax From 194b1f750c0203db6a3b50cca5d36a626ccc8c44 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 17 Mar 2016 15:41:47 +0100 Subject: [PATCH 7/8] Changed norm of MOs when the fitcusp is used --- ezfio_config/qmc.config | 3 +- ocaml/Default.ml | 22 +++++--- ocaml/Input.ml | 96 +++++++++++++++++------------------ ocaml/Md5.ml | 3 +- ocaml/Qmcchem_edit.ml | 40 +++++++-------- src/AO/ao_full.irp.f | 1 + src/ezfio_interface.irp.f | 3 +- src/mo.irp.f | 104 ++++++++++++++++++++++++++++++++++++++ src/nuclei.irp.f | 11 ++-- src/simulation.irp.f | 21 +++++--- 10 files changed, 208 insertions(+), 96 deletions(-) diff --git a/ezfio_config/qmc.config b/ezfio_config/qmc.config index f1392ce..29aaac9 100644 --- a/ezfio_config/qmc.config +++ b/ezfio_config/qmc.config @@ -33,7 +33,6 @@ nuclei nucl_label character*(32) (nuclei_nucl_num) nucl_charge real (nuclei_nucl_num) nucl_coord real (nuclei_nucl_num,3) - nucl_fitcusp_radius real (nuclei_nucl_num) spindeterminants n_det_alpha integer @@ -55,7 +54,7 @@ simulation equilibration logical http_server character*(128) do_jast logical - do_nucl_fitcusp logical + nucl_fitcusp_factor real method character*(32) block_time integer sampling character*(32) diff --git a/ocaml/Default.ml b/ocaml/Default.ml index 37904a3..65b542d 100644 --- a/ocaml/Default.ml +++ b/ocaml/Default.ml @@ -1,23 +1,29 @@ open Core.Std;; -let simulation_do_nucl_fitcusp = lazy( - if (not (Ezfio.has_pseudo_do_pseudo ())) then - not (Ezfio.get_pseudo_do_pseudo ()) +let simulation_nucl_fitcusp_factor = lazy( + let default = + 0.5 + in + if (Ezfio.has_pseudo_do_pseudo ()) then + if (Ezfio.get_pseudo_do_pseudo ()) then + 0. + else + default else - true + default ) let electrons_elec_walk_num = lazy ( 30 ) let electrons_elec_walk_num_tot = lazy ( 10000 ) let jastrow_jast_type = lazy ( "None" ) let simulation_block_time = lazy ( 30 ) -let simulation_ci_threshold = lazy ( 1.e-8 ) +let simulation_ci_threshold = lazy ( 1.e-8 ) let simulation_method = lazy ( "VMC" ) let simulation_sampling = lazy ( "Langevin" ) let simulation_stop_time = lazy ( 3600 ) let simulation_time_step = lazy ( 0.15 ) -let simulation_srmc_projection_time = lazy ( 1. ) +let simulation_srmc_projection_time = lazy ( 1. ) let reset_defaults () = List.iter ~f:(fun x -> Sys.remove ( (Lazy.force Qputils.ezfio_filename) ^ x)) @@ -26,9 +32,9 @@ let reset_defaults () = "/jastrow/jast_type" ; "/simulation/block_time" ; "/simulation/ci_threshold" ; - "/simulation/do_nucl_fitcusp" ; "/simulation/method" ; "/simulation/sampling" ; "/simulation/stop_time" ; - "/simulation/time_step" ] + "/simulation/time_step" ; + "/simulation/nucl_fitcusp_factor" ] diff --git a/ocaml/Input.ml b/ocaml/Input.ml index a61444b..ad9c770 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -66,81 +66,69 @@ end = struct end -module Fitcusp : sig - type t = bool + +module Fitcusp_factor : sig + + type t = float val doc : string val read : unit -> t val write : t -> unit - val to_bool : t -> bool - val of_bool : bool -> t - val to_int : t -> int - val of_int : int -> t + val to_float : t -> float + val of_float : float -> t val to_string : t -> string val of_string : string -> t end = struct - type t = bool + type t = float - let doc = "Correct wave function to verify electron-nucleus cusp condition" + let doc = "Correct wave function to verify electron-nucleus cusp condition. +Fit is done for r < r_c(f) where r_c(f) = (1s orbital radius) x f. Value of f" - let of_bool x = x + let of_float x = + if (x < 0.) then + failwith "Fitcusp_factor should be >= 0."; + if (x > 10.) then + failwith "Fitcusp_factor is too large."; + x - let to_bool x = x + let to_float x = x let read () = - let _ = - Lazy.force Qputils.ezfio_filename - in - if (not (Ezfio.has_simulation_do_nucl_fitcusp ())) then - Lazy.force Default.simulation_do_nucl_fitcusp - |> Ezfio.set_simulation_do_nucl_fitcusp ; - Ezfio.get_simulation_do_nucl_fitcusp () - |> of_bool + ignore @@ + Lazy.force Qputils.ezfio_filename ; + if (not (Ezfio.has_simulation_nucl_fitcusp_factor ())) then + begin + let factor = + Lazy.force Default.simulation_nucl_fitcusp_factor ; + in + Ezfio.set_simulation_nucl_fitcusp_factor factor + end ; + Ezfio.get_simulation_nucl_fitcusp_factor () + |> of_float let write t = let _ = Lazy.force Qputils.ezfio_filename in - let () = - match (Pseudo.read () |> Pseudo.to_bool, to_bool t) with - | (true, true) -> failwith "Pseudopotentials and Fitcusp are incompatible" - | _ -> () - in - to_bool t - |> Ezfio.set_simulation_do_nucl_fitcusp + to_float t + |> Ezfio.set_simulation_nucl_fitcusp_factor let to_string t = - to_bool t - |> Bool.to_string + to_float t + |> Float.to_string let of_string t = try - String.lowercase t - |> Bool.of_string - |> of_bool + Float.of_string t + |> of_float with | Invalid_argument msg -> failwith msg - - let to_int t = - let t = - to_bool t - in - if t then 1 - else 0 - - - let of_int = function - | 0 -> false - | 1 -> true - | _ -> failwith "Expected 0 or 1" - - end module Block_time : sig @@ -855,8 +843,6 @@ let validate () = Time_step.read () and jast_type = Jastrow_type.read () - and do_fitcusp = - Fitcusp.read () and do_pseudo = Pseudo.read () in @@ -915,13 +901,23 @@ let validate () = | _ -> () in - (* Fitcusp is not recommended with pseudo *) + (* Fitcusp is incompatible with pseudo *) let () = - match (Pseudo.to_bool do_pseudo, Fitcusp.to_bool do_fitcusp) with - | (true, true) -> warn "Fitcusp is incompatible with Pseudopotentials" + let f = + Fitcusp_factor.read () + |> Fitcusp_factor.to_float + in + match (Pseudo.to_bool do_pseudo, f > 0.) with + | (true, true) -> + begin + warn "Electron-nucleus cusp fitting is incompatible with Pseudopotentials."; + Fitcusp_factor.of_float 0. + |> Fitcusp_factor.write + end | _ -> () in + (* Other Checks *) let () = let _ = diff --git a/ocaml/Md5.ml b/ocaml/Md5.ml index 8abd5ca..48d0b9b 100644 --- a/ocaml/Md5.ml +++ b/ocaml/Md5.ml @@ -37,10 +37,9 @@ let files_to_track = [ "mo_basis/mo_tot_num" ; "nuclei/nucl_charge.gz" ; "nuclei/nucl_coord.gz" ; - "nuclei/nucl_fitcusp_radius.gz" ; "nuclei/nucl_num" ; "simulation/ci_threshold" ; - "simulation/do_nucl_fitcusp" ; + "simulation/nucl_fitcusp_factor" ; "simulation/jast_a_up_dn" ; "simulation/jast_a_up_up" ; "simulation/jast_b_up_dn" ; diff --git a/ocaml/Qmcchem_edit.ml b/ocaml/Qmcchem_edit.ml index 9b75263..2f6a6e7 100644 --- a/ocaml/Qmcchem_edit.ml +++ b/ocaml/Qmcchem_edit.ml @@ -20,7 +20,7 @@ type field = | Walk_num | Walk_num_tot | Stop_time - | Fitcusp + | Fitcusp_factor | Method | Sampling | Ref_energy @@ -54,8 +54,8 @@ let get field = option_to_string Walk_num_tot.read Walk_num_tot.to_string Walk_num_tot.doc | Stop_time -> option_to_string Stop_time.read Stop_time.to_string Stop_time.doc - | Fitcusp -> - option_to_string Fitcusp.read Fitcusp.to_string Fitcusp.doc + | Fitcusp_factor -> + option_to_string Fitcusp_factor.read Fitcusp_factor.to_string Fitcusp_factor.doc | Method -> option_to_string Method.read Method.to_string Method.doc | Sampling -> @@ -130,19 +130,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.(of_int , 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.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 = [ @@ -155,7 +155,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = Ref_energy ; Walk_num ; Walk_num_tot ; - Fitcusp ; + Fitcusp_factor ; CI_threshold ; Jastrow_type ; Properties ; @@ -214,7 +214,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = begin match f with | Stop_time -> Stop_time.(of_string s |> write) - | Fitcusp -> Fitcusp.(of_string s |> write) + | Fitcusp_factor -> Fitcusp_factor.(of_string s |> write) | Block_time -> Block_time.(of_string s |> write) | Method -> Method.(of_string s |> write) | Ref_energy -> Ref_energy.(of_string s |> write) @@ -271,8 +271,8 @@ let spec = empty +> flag "c" no_arg ~doc:(" Clear blocks") - +> flag "f" (optional int) - ~doc:("0|1 "^Input.Fitcusp.doc) + +> flag "f" (optional float) + ~doc:("float "^Input.Fitcusp_factor.doc) +> flag "t" (optional int) ~doc:("seconds "^Input.Stop_time.doc) +> flag "l" (optional int) diff --git a/src/AO/ao_full.irp.f b/src/AO/ao_full.irp.f index 3ff16ba..c59abb7 100644 --- a/src/AO/ao_full.irp.f +++ b/src/AO/ao_full.irp.f @@ -73,6 +73,7 @@ BEGIN_PROVIDER [ logical, primitives_reduced ] PROVIDE ao_power PROVIDE ao_coef PROVIDE ao_nucl + PROVIDE mo_fitcusp_normalization_before do i=1,ao_num if (ao_oned_p(i) /= 0.) then l=ao_power(i,1)+ao_power(i,2)+ao_power(i,3) diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index 32e4f6d..08c9cfe 100644 --- a/src/ezfio_interface.irp.f +++ b/src/ezfio_interface.irp.f @@ -6,7 +6,6 @@ data = [ \ ("nuclei_nucl_num" , "integer" , "" ), ("nuclei_nucl_charge" , "real" , "(nucl_num)" ), ("nuclei_nucl_coord" , "real" , "(nucl_num,3)" ), -("nuclei_nucl_fitcusp_radius" , "real" , "(nucl_num)" ), ("mo_basis_mo_coef" , "real" , "(ao_num,mo_tot_num)" ), ("electrons_elec_fitcusp_radius" , "real" , "" ), ("electrons_elec_alpha_num" , "integer" , "" ), @@ -38,9 +37,9 @@ data = [ \ ("simulation_time_step" , "real" , "" ), ("simulation_srmc_projection_time" , "real" , "" ), ("simulation_method" , "character*(32)", "" ), +("simulation_nucl_fitcusp_factor" , "real" , "" ), ("simulation_save_data" , "logical" , "" ), ("simulation_print_level" , "integer" , "" ), -("simulation_do_nucl_fitcusp" , "logical" , "" ), ("simulation_sampling" , "character*(32)", "" ), ("simulation_ci_threshold" , "double precision" , "" ), ("simulation_http_server" , "character*(128)", "" ), diff --git a/src/mo.irp.f b/src/mo.irp.f index 233137a..2dfa9c6 100644 --- a/src/mo.irp.f +++ b/src/mo.irp.f @@ -47,6 +47,7 @@ END_PROVIDER END_DOC mo_scale = 1.d0/(0.4d0*log(float(elec_num+1))) mo_norm = mo_scale*mo_scale + END_PROVIDER @@ -273,6 +274,15 @@ END_PROVIDER enddo endif + do i=1,mo_num + do j=1,elec_num + mo_value_transp(i,j) *= mo_cusp_rescale(i) + mo_grad_transp_x(i,j) *= mo_cusp_rescale(i) + mo_grad_transp_y(i,j) *= mo_cusp_rescale(i) + mo_grad_transp_z(i,j) *= mo_cusp_rescale(i) + mo_lapl_transp(i,j) *= mo_cusp_rescale(i) + enddo + enddo END_PROVIDER @@ -401,6 +411,7 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ] integer :: i, j, k, l real :: ao_value_at_nucl_no_S(ao_num) + PROVIDE mo_fitcusp_normalization_before do k=1,nucl_num point(1) = nucl_coord(k,1) point(2) = nucl_coord(k,2) @@ -466,6 +477,99 @@ END_PROVIDER FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p SOFT_TOUCH point +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_before, (mo_tot_num) ] + implicit none + BEGIN_DOC + ! Renormalization factor of MOs due to cusp fitting + END_DOC + include 'constants.F' + integer :: i,j,k,l + double precision :: dr, r, f, t + integer, save :: ifirst = 0 + + if (ifirst == 0) then + ifirst = 1 + mo_fitcusp_normalization_before = 0.d0 + do k=1,nucl_num + dr = nucl_fitcusp_radius(k)*1.d-2 + point(1) = nucl_coord(k,1) + point(2) = nucl_coord(k,2) + point(3) = nucl_coord(k,3)-dr + do l=1,101 + r = point(3) + dr + point(3) = r + TOUCH point + f = dfour_pi*r*r*dr + do i=1,mo_tot_num + mo_fitcusp_normalization_before(i) += f*mo_value_p(i)**2 + enddo + enddo + enddo + endif + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num) ] + implicit none + BEGIN_DOC + ! Renormalization factor of MOs due to cusp fitting + END_DOC + include 'constants.F' + integer :: i,j,k,l + double precision :: dr, r, f, t, t2 + integer, save :: ifirst = 0 + + PROVIDE primitives_reduced + if (ifirst == 0) then + ifirst = 1 + mo_fitcusp_normalization_after = 0.d0 + do k=1,nucl_num + dr = nucl_fitcusp_radius(k)*1.d-2 + point(1) = nucl_coord(k,1) + point(2) = nucl_coord(k,2) + point(3) = nucl_coord(k,3)- dr + do l=1,101 + point(3) = point(3)+ dr + TOUCH point nucl_fitcusp_param primitives_reduced mo_coef + r = point(3) + f = dfour_pi*r*r*dr + do i=1,mo_tot_num + t = 0.d0 + do j=1,ao_num + if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then + t = t + mo_coef(j,i) * ao_value_p(j) + endif + enddo + t = t + nucl_fitcusp_param(1,i,k) + & + r * (nucl_fitcusp_param(2,i,k) + & + r * (nucl_fitcusp_param(3,i,k) + & + r * nucl_fitcusp_param(4,i,k) )) + mo_fitcusp_normalization_after(i) += t*t*f + enddo + enddo + enddo + endif + +END_PROVIDER + +BEGIN_PROVIDER [ real, mo_cusp_rescale, (mo_tot_num) ] + implicit none + BEGIN_DOC + ! Rescaling coefficient to normalize MOs after applying fitcusp + END_DOC + integer :: i + if (do_nucl_fitcusp) then + do i=1,mo_tot_num +! mo_cusp_rescale(i) = dsqrt(mo_fitcusp_normalization_before(i) / mo_fitcusp_normalization_after(i)) + mo_cusp_rescale(i) = 1.d0/dsqrt(1.d0 - mo_fitcusp_normalization_before(i) + mo_fitcusp_normalization_after(i)) + enddo + else + mo_cusp_rescale = 1.d0 + endif + END_PROVIDER diff --git a/src/nuclei.irp.f b/src/nuclei.irp.f index 39d84a5..cdcbc8b 100644 --- a/src/nuclei.irp.f +++ b/src/nuclei.irp.f @@ -125,17 +125,20 @@ BEGIN_PROVIDER [ real, nucl_fitcusp_radius, (nucl_num) ] BEGIN_DOC ! Distance threshold for the fit END_DOC - real :: def(nucl_num) + real :: def(nucl_num), factor integer :: k + real, parameter :: a = 1.74891 + real, parameter :: b = 0.126057 - if (.not.do_nucl_fitcusp) then + if (.not. do_nucl_fitcusp) then nucl_fitcusp_radius = 0.d0 return endif + do k=1,nucl_num - nucl_fitcusp_radius(k) = .5/nucl_charge(k) + nucl_fitcusp_radius(k) = nucl_fitcusp_factor/(a*nucl_charge(k)+b) enddo - call get_nuclei_nucl_fitcusp_radius(nucl_fitcusp_radius) + ! Avoid dummy atoms do k=1,nucl_num if (nucl_charge(k) < 5.d-1) then diff --git a/src/simulation.irp.f b/src/simulation.irp.f index dfdb084..3608e08 100644 --- a/src/simulation.irp.f +++ b/src/simulation.irp.f @@ -250,16 +250,21 @@ BEGIN_PROVIDER [ character*(64), hostname] END_PROVIDER -BEGIN_PROVIDER [ logical, do_nucl_fitcusp ] - implicit none - BEGIN_DOC - ! If true, do the fit of the electron-nucleus cusp - END_DOC - do_nucl_fitcusp = .True. - call get_simulation_do_nucl_fitcusp(do_nucl_fitcusp) - call linfo(irp_here,'do_nucl_fitcusp',do_nucl_fitcusp) + BEGIN_PROVIDER [ real, nucl_fitcusp_factor ] +&BEGIN_PROVIDER [ logical, do_nucl_fitcusp ] + implicit none + BEGIN_DOC + ! The electron-nucleus cusp fitting is done between 0 and r_c, + ! where r_c is chosen as nucl_fitcusp_factor * (radius_of_1s AO) + END_DOC + nucl_fitcusp_factor = 0. + call get_simulation_nucl_fitcusp_factor(nucl_fitcusp_factor) + do_nucl_fitcusp = nucl_fitcusp_factor > 0. + call info(irp_here,'nucl_fitcusp_factor',nucl_fitcusp_factor) + END_PROVIDER + BEGIN_PROVIDER [ integer, vmc_algo ] implicit none From 756f2ccea4cf6e0685426c5c2d5be33bfdcd5800 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 17 Mar 2016 15:51:45 +0100 Subject: [PATCH 8/8] Added Rousset FKMC algorithm --- ocaml/Default.ml | 6 +- ocaml/Input.ml | 17 +- src/SAMPLING/block.irp.f | 6 + src/SAMPLING/dmc_step.irp.f | 14 +- src/SAMPLING/fkmc_step.irp.f | 375 +++++++++++++++++++++++++++++++++++ src/SAMPLING/srmc_step.irp.f | 4 +- src/mo.irp.f | 3 +- src/simulation.irp.f | 9 +- src/types.F | 3 +- 9 files changed, 417 insertions(+), 20 deletions(-) create mode 100644 src/SAMPLING/fkmc_step.irp.f diff --git a/ocaml/Default.ml b/ocaml/Default.ml index 65b542d..12f5f86 100644 --- a/ocaml/Default.ml +++ b/ocaml/Default.ml @@ -3,7 +3,7 @@ open Core.Std;; let simulation_nucl_fitcusp_factor = lazy( let default = - 0.5 + 1. in if (Ezfio.has_pseudo_do_pseudo ()) then if (Ezfio.get_pseudo_do_pseudo ()) then @@ -14,8 +14,8 @@ let simulation_nucl_fitcusp_factor = lazy( default ) -let electrons_elec_walk_num = lazy ( 30 ) -let electrons_elec_walk_num_tot = lazy ( 10000 ) +let electrons_elec_walk_num = lazy ( 100 ) +let electrons_elec_walk_num_tot = lazy ( 1000 ) let jastrow_jast_type = lazy ( "None" ) let simulation_block_time = lazy ( 30 ) let simulation_ci_threshold = lazy ( 1.e-8 ) diff --git a/ocaml/Input.ml b/ocaml/Input.ml index ad9c770..01e8f97 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -387,7 +387,7 @@ end module Method : sig - type t = VMC | DMC | SRMC + type t = VMC | DMC | SRMC | FKMC val doc : string val read : unit -> t val write : t -> unit @@ -396,21 +396,23 @@ module Method : sig end = struct - type t = VMC | DMC | SRMC + type t = VMC | DMC | SRMC | FKMC - let doc = "QMC Method : [ VMC | DMC | SRMC ]" + let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC ]" let of_string = function | "VMC" | "vmc" -> VMC | "DMC" | "dmc" -> DMC | "SRMC" | "srmc" -> SRMC - | x -> failwith ("Method should be [ VMC | DMC | SRMC ], not "^x^".") + | "FKMC" | "fkmc" -> FKMC + | x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC ], not "^x^".") let to_string = function | VMC -> "VMC" | DMC -> "DMC" | SRMC -> "SRMC" + | FKMC -> "FKMC" let read () = @@ -851,10 +853,12 @@ let validate () = let () = match (sampling, meth, Pseudo.to_bool do_pseudo) with | (Sampling.Brownian, Method.DMC, true) + | (Sampling.Brownian, Method.FKMC, true) | (Sampling.Brownian, Method.SRMC, true) -> if ( (Time_step.to_float ts) >= 0.5 ) then warn ( "Time step seems large for "^(Method.to_string meth) ) | (Sampling.Brownian, Method.SRMC, false) + | (Sampling.Brownian, Method.FKMC, false) | (Sampling.Brownian, Method.DMC, false) -> if ( (Time_step.to_float ts) >= 0.01 ) then warn ( "Time step seems large for "^(Method.to_string meth) ) @@ -865,8 +869,9 @@ let validate () = if ( (Time_step.to_float ts) <= 0.01 ) then warn "Time step seems small for Langevin sampling." | (Sampling.Langevin, Method.SRMC, _) + | (Sampling.Langevin, Method.FKMC, _) | (Sampling.Langevin, Method.DMC, _) -> - failwith "Lanvegin sampling is incompatible with DMC/SRMC" + failwith "Lanvegin sampling is incompatible with DMC" in @@ -874,6 +879,7 @@ let validate () = let () = match (meth, Ref_energy.(read () |> to_float) ) with | (Method.SRMC,0.) + | (Method.FKMC,0.) | (Method.DMC,0.) -> failwith ("E_ref should not be zero in "^(Method.to_string meth) ) | _ -> () in @@ -888,6 +894,7 @@ let validate () = let () = match (meth, Property.(calc E_loc)) with | (Method.SRMC, false) + | (Method.FKMC, false) | (Method.DMC, false) -> failwith ( "E_loc should be sampled in "^(Method.to_string meth) ) | (Method.VMC, false) -> warn "Sampling of E_loc is not activated in input" | _ -> () diff --git a/src/SAMPLING/block.irp.f b/src/SAMPLING/block.irp.f index 092eb17..a732941 100644 --- a/src/SAMPLING/block.irp.f +++ b/src/SAMPLING/block.irp.f @@ -30,6 +30,12 @@ t = """ $X_block_walk = $X_srmc_block_walk $X_2_block_walk = $X_2_srmc_block_walk endif + else if (qmc_method == t_FKMC) then + PROVIDE E_loc_fkmc_block_walk + if (calc_$X) then + $X_block_walk = $X_fkmc_block_walk + $X_2_block_walk = $X_2_fkmc_block_walk + endif endif END_PROVIDER diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f index e11662e..2764609 100644 --- a/src/SAMPLING/dmc_step.irp.f +++ b/src/SAMPLING/dmc_step.irp.f @@ -195,7 +195,7 @@ END_SHELL do k=1,walk_num_dmc sum_weight += dmc_weight(k) enddo - E0 = E_ref - log(real(walk_num_dmc)/real(walk_num)) * 0.1d0/dtime_step + E0 = E_ref - log(sum_weight/real(walk_num)) * 0.1d0 /dtime_step ! Branching integer :: ipos(walk_num_dmc_max), walk_num_dmc_new @@ -214,20 +214,22 @@ END_SHELL ipos(k) = k enddo + walk_num_dmc_new = walk_num_dmc do k=1,walk_num_dmc r = qmc_ranf() if (dmc_weight(k) > 1.d0) then if ( 1.d0+r < dmc_weight(k) ) then - walk_num_dmc = walk_num_dmc+1 - ipos(walk_num_dmc) = k + walk_num_dmc_new = walk_num_dmc_new+1 + ipos(walk_num_dmc_new) = k endif else if ( r > dmc_weight(k) ) then - ipos(k) = ipos(walk_num_dmc) - walk_num_dmc = walk_num_dmc-1 + ipos(k) = ipos(walk_num_dmc_new) + walk_num_dmc_new = walk_num_dmc_new-1 endif endif enddo + walk_num_dmc = walk_num_dmc_new integer :: ipm do k=1,walk_num_dmc @@ -328,7 +330,7 @@ BEGIN_PROVIDER [ integer, walk_num_dmc_max ] BEGIN_DOC ! Max number of walkers in DMC END_DOC - walk_num_dmc_max = 3 * walk_num + walk_num_dmc_max = max(3 * walk_num, 30) END_PROVIDER diff --git a/src/SAMPLING/fkmc_step.irp.f b/src/SAMPLING/fkmc_step.irp.f new file mode 100644 index 0000000..878af69 --- /dev/null +++ b/src/SAMPLING/fkmc_step.irp.f @@ -0,0 +1,375 @@ +! Providers of *_fkmc_block_walk +!============================== +BEGIN_SHELL [ /usr/bin/python ] +from properties import * + +t = """ + BEGIN_PROVIDER [ $T, $X_fkmc_block_walk $D1 ] +&BEGIN_PROVIDER [ $T, $X_fkmc_block_walk_kahan $D2 ] +&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk $D1 ] +&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk_kahan $D2 ] + implicit none + BEGIN_DOC +! fkMC averages of $X. Computed in E_loc_fkmc_block_walk + END_DOC + $X_fkmc_block_walk = 0.d0 + $X_fkmc_block_walk_kahan = 0.d0 + $X_2_fkmc_block_walk = 0.d0 + $X_2_fkmc_block_walk_kahan = 0.d0 +END_PROVIDER +""" +for p in properties: + if p[1] != 'e_loc': + if p[2] == "": + D1 = "" + D2 = ", (3)" + else: + D1 = ", ("+p[2][1:-1]+")" + D2 = ", ("+p[2][1:-1]+",3)" + print t.replace("$X",p[1]).replace("$T",p[0]).replace("$D1",D1).replace("$D2",D2) + +END_SHELL + + + + BEGIN_PROVIDER [ double precision, E_loc_fkmc_block_walk ] +&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk ] +&BEGIN_PROVIDER [ double precision, E_loc_fkmc_block_walk_kahan, (3) ] +&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk_kahan, (3) ] + implicit none + include '../types.F' + BEGIN_DOC +! Properties averaged over the block using the FKMC method + END_DOC + + integer, parameter :: BIRTH=1, DEATH=2 + real, allocatable :: elec_coord_tmp(:,:,:) + integer :: mod_align + double precision :: E_loc_save(walk_num_dmc_max) + double precision :: E_loc_save_tmp(walk_num_dmc_max) + double precision :: psi_value_save(walk_num) + double precision :: psi_value_save_tmp(walk_num) + double precision :: fkmc_weight(walk_num) + double precision :: delta(walk_num) + double precision, allocatable :: psi_grad_psi_inv_save(:,:,:) + double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:) + double precision, allocatable :: fkmc_clock_tmp(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save_tmp + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save_tmp + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_weight + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_clock_tmp + allocate ( psi_grad_psi_inv_save(elec_num_8,3,walk_num), & + psi_grad_psi_inv_save_tmp(elec_num_8,3,walk_num), & + elec_coord_tmp(mod_align(elec_num+1),3,walk_num), & + fkmc_clock_tmp(2,walk_num) ) + psi_value_save = 0.d0 + psi_value_save_tmp = 0.d0 + fkmc_weight = 1.d0 + +! Initialization + if (vmc_algo /= t_Brownian) then + call abrt(irp_here,'FKMC should run with Brownian algorithm') + endif + + integer :: k, i_walk, i_step + +BEGIN_SHELL [ /usr/bin/python ] +from properties import * +t = """ + if (calc_$X) then + !DIR$ VECTOR ALIGNED + $X_fkmc_block_walk = 0.d0 + !DIR$ VECTOR ALIGNED + $X_fkmc_block_walk_kahan = 0.d0 + !DIR$ VECTOR ALIGNED + $X_2_fkmc_block_walk = 0.d0 + !DIR$ VECTOR ALIGNED + $X_2_fkmc_block_walk_kahan = 0.d0 + endif +""" +for p in properties: + print t.replace("$X",p[1]) +END_SHELL + + logical :: loop + integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max + + loop = .True. + call system_clock(cpu0, count_rate, count_max) + cpu2 = cpu0 + + block_weight = 0.d0 + + real, external :: accep_rate + double precision :: thr + + thr = 2.d0/time_step_sq + + logical :: first_loop + first_loop = .True. + + do while (loop) + + ! Every walker makes a step + do i_walk=1,walk_num + + if (.not.first_loop) then + integer :: i,j,l + do l=1,3 + do i=1,elec_num+1 + elec_coord(i,l) = elec_coord_full(i,l,i_walk) + enddo + do i=1,elec_num + psi_grad_psi_inv_x(i) = psi_grad_psi_inv_save(i,1,i_walk) + psi_grad_psi_inv_y(i) = psi_grad_psi_inv_save(i,2,i_walk) + psi_grad_psi_inv_z(i) = psi_grad_psi_inv_save(i,3,i_walk) + enddo + psi_value = psi_value_save(i_walk) + E_loc = E_loc_save(i_walk) + enddo + SOFT_TOUCH elec_coord psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z psi_value E_loc + else + do l=1,3 + do i=1,elec_num+1 + elec_coord(i,l) = elec_coord_full(i,l,i_walk) + enddo + enddo + TOUCH elec_coord + E_loc_save(i_walk) = E_loc + psi_value_save(i_walk) = psi_value + endif + + double precision :: p,q + real :: delta_x + logical :: accepted + call brownian_step(p,q,accepted,delta_x) + + if ( psi_value * psi_value_save(i_walk) >= 0.d0 ) then + delta(i_walk) = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p + if ( delta(i_walk) > thr ) then + delta(i_walk) = thr + else if ( delta(i_walk) < -thr ) then + delta(i_walk) = -thr + endif + fkmc_weight(i_walk) = dexp(-dtime_step*delta(i_walk)) + elec_coord(elec_num+1,1) += p*time_step + elec_coord(elec_num+1,2) = E_loc + elec_coord(elec_num+1,3) = fkmc_weight(i_walk) + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,i_walk) = elec_coord(i,l) + enddo + enddo + do i=1,elec_num + psi_grad_psi_inv_save(i,1,i_walk) = psi_grad_psi_inv_x(i) + psi_grad_psi_inv_save(i,2,i_walk) = psi_grad_psi_inv_y(i) + psi_grad_psi_inv_save(i,3,i_walk) = psi_grad_psi_inv_z(i) + enddo + + psi_value_save(i_walk) = psi_value + E_loc_save(i_walk) = E_loc + +BEGIN_SHELL [ /usr/bin/python ] +from properties import * +t = """ + if (calc_$X) then + ! Kahan's summation algorithm to compute these sums reducing the rounding error: + ! $X_fkmc_block_walk += $X * fkmc_weight(i_walk) + ! $X_2_fkmc_block_walk += $X_2 * fkmc_weight(i_walk) + ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm + + $X_fkmc_block_walk_kahan($D2 3) = $X * fkmc_weight(i_walk) - $X_fkmc_block_walk_kahan($D2 1) + $X_fkmc_block_walk_kahan($D2 2) = $X_fkmc_block_walk $D1 + $X_fkmc_block_walk_kahan($D2 3) + $X_fkmc_block_walk_kahan($D2 1) = ($X_fkmc_block_walk_kahan($D2 2) - $X_fkmc_block_walk $D1 ) & + - $X_fkmc_block_walk_kahan($D2 3) + $X_fkmc_block_walk $D1 = $X_fkmc_block_walk_kahan($D2 2) + + + $X_2_fkmc_block_walk_kahan($D2 3) = $X_2 * fkmc_weight(i_walk) - $X_2_fkmc_block_walk_kahan($D2 1) + $X_2_fkmc_block_walk_kahan($D2 2) = $X_2_fkmc_block_walk $D1 + $X_2_fkmc_block_walk_kahan($D2 3) + $X_2_fkmc_block_walk_kahan($D2 1) = ($X_2_fkmc_block_walk_kahan($D2 2) - $X_2_fkmc_block_walk $D1 ) & + - $X_2_fkmc_block_walk_kahan($D2 3) + $X_2_fkmc_block_walk $D1 = $X_2_fkmc_block_walk_kahan($D2 2) + endif +""" +for p in properties: + if p[2] == "": + D1 = "" + D2 = "" + else: + D1 = "("+":"*(p[2].count(',')+1)+")" + D2 = ":"*(p[2].count(',')+1)+"," + print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2) + +END_SHELL + + block_weight += fkmc_weight(i_walk) + + else + fkmc_weight(i_walk) = 0.d0 + delta(i_walk) = 1.d5 + endif + + enddo + + ! Compute the new weight of the population + double precision :: sum_weight + sum_weight = 0.d0 + do k=1,walk_num + sum_weight += fkmc_weight(k) + enddo + + do k=1,walk_num + do l=1,3 + do i=1,elec_num+1 + elec_coord_tmp(i,l,k) = elec_coord_full(i,l,k) + enddo + do i=1,elec_num + psi_grad_psi_inv_save_tmp(i,l,k) = psi_grad_psi_inv_save(i,l,k) + enddo + enddo + psi_value_save_tmp(k) = psi_value_save(k) + E_loc_save_tmp(k) = E_loc_save(k) + if (fkmc_weight(k) == 0.d0) then + fkmc_clock(DEATH,k) = -1.d0 + endif + if ( delta(k) <= 0.d0 ) then + fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k) +time_step * delta(k) + fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k) + else + fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k) + fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k) -time_step * delta(k) + endif + enddo + +! Reconfiguration +! =============== + + ! Identify first which walkers will be killed to place branched walkers there + ! later + + double precision, external :: qmc_ranf + integer :: ipm, m + integer :: killed(walk_num) + + m=1 + do k=1,walk_num + fkmc_clock(DEATH,k) = fkmc_clock_tmp(DEATH,k) + if (fkmc_clock_tmp(DEATH,k) <= 0.d0) then + killed(m) = k + m += 1 + fkmc_clock(DEATH,k) = -dlog(qmc_ranf()) + fkmc_clock(BIRTH,k) = -dlog(qmc_ranf()) + ipm = k + do while (ipm == k) + ipm = 1 + int (walk_num*qmc_ranf()) + enddo + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm) + enddo + do i=1,elec_num + psi_grad_psi_inv_save(i,l,k) = psi_grad_psi_inv_save_tmp(i,l,ipm) + enddo + enddo + psi_value_save(k) = psi_value_save_tmp(ipm) + E_loc_save(k) = E_loc_save_tmp(ipm) + endif + enddo + killed(m) = 0 + + m=1 + do k=1,walk_num + fkmc_clock(BIRTH,k) = fkmc_clock_tmp(BIRTH,k) + if (fkmc_clock_tmp(BIRTH,k) <= 0.d0) then + fkmc_clock(BIRTH,k) = -dlog(qmc_ranf()) + if (killed(m) == 0) then + ipm = k + do while (ipm == k) + ipm = 1 + int (walk_num*qmc_ranf()) + enddo + else + ipm = killed(m) + m +=1 + endif + fkmc_clock(BIRTH,ipm) = -dlog(qmc_ranf()) + fkmc_clock(DEATH,ipm) = -dlog(qmc_ranf()) + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,ipm) = elec_coord_tmp(i,l,k) + enddo + do i=1,elec_num + psi_grad_psi_inv_save(i,l,ipm) = psi_grad_psi_inv_save_tmp(i,l,k) + enddo + enddo + psi_value_save(ipm) = psi_value_save_tmp(k) + E_loc_save(ipm) = E_loc_save_tmp(k) + + endif + + enddo + + + + call system_clock(cpu1, count_rate, count_max) + if (cpu1 < cpu0) then + cpu1 = cpu1+cpu0 + endif + loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate) + if (cpu1-cpu2 > count_rate) then + integer :: do_run + call get_running(do_run) + loop = do_run == t_Running + cpu2 = cpu1 + endif + +! Update E_ref to take into account the weight of the population + E_ref -= dlog(sum_weight / dble(walk_num) ) / time_step + SOFT_TOUCH elec_coord_full E_ref + + first_loop = .False. + + enddo + + double precision :: factor + factor = 1.d0/block_weight + SOFT_TOUCH block_weight + +BEGIN_SHELL [ /usr/bin/python ] +from properties import * +t = """ + if (calc_$X) then + $X_fkmc_block_walk *= factor + $X_2_fkmc_block_walk *= factor + endif +""" +for p in properties: + print t.replace("$X",p[1]) +END_SHELL + + deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp, & + fkmc_clock_tmp ) + +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, fkmc_clock, (2,walk_num) ] + implicit none + BEGIN_DOC + ! Branching clocks for the FKMC algotithm. (1,:) is the birth clock and + ! (2,:) is the death clock. + END_DOC + integer :: i + double precision, external :: qmc_ranf + do i=1, walk_num + fkmc_clock(1,i) = -dlog(qmc_ranf()) + fkmc_clock(2,i) = -dlog(qmc_ranf()) + enddo + +END_PROVIDER + diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f index 3afaab9..fc98ab5 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -133,7 +133,9 @@ END_SHELL elec_coord(i,l) = elec_coord_full(i,l,i_walk) enddo enddo - SOFT_TOUCH elec_coord + TOUCH elec_coord + psi_value_save(i_walk) = psi_value + E_loc_save(i_walk) = E_loc endif double precision :: p,q diff --git a/src/mo.irp.f b/src/mo.irp.f index 2dfa9c6..5e43659 100644 --- a/src/mo.irp.f +++ b/src/mo.irp.f @@ -47,7 +47,6 @@ END_PROVIDER END_DOC mo_scale = 1.d0/(0.4d0*log(float(elec_num+1))) mo_norm = mo_scale*mo_scale - END_PROVIDER @@ -536,7 +535,7 @@ BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num) TOUCH point nucl_fitcusp_param primitives_reduced mo_coef r = point(3) f = dfour_pi*r*r*dr - do i=1,mo_tot_num + do i=1,mo_num t = 0.d0 do j=1,ao_num if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then diff --git a/src/simulation.irp.f b/src/simulation.irp.f index 3608e08..0995e2d 100644 --- a/src/simulation.irp.f +++ b/src/simulation.irp.f @@ -148,7 +148,7 @@ BEGIN_PROVIDER [ integer, qmc_method ] implicit none include 'types.F' BEGIN_DOC - ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC + ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC, t_FKMC END_DOC character*(32) :: method method = types(t_VMC) @@ -160,8 +160,10 @@ BEGIN_PROVIDER [ integer, qmc_method ] qmc_method = t_DMC else if (method == types(t_SRMC)) then qmc_method = t_SRMC + else if (method == types(t_FKMC)) then + qmc_method = t_FKMC else - call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC )') + call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC )') endif call cinfo(irp_here,'qmc_method',trim(method)) @@ -288,6 +290,9 @@ BEGIN_PROVIDER [ integer, vmc_algo ] if (qmc_method == t_SRMC) then stop 'Langevin incompatible with SRMC' endif + if (qmc_method == t_FKMC) then + stop 'Langevin incompatible with FKMC' + endif else if (Sampling == types(t_MTM)) then vmc_algo = t_MTM else diff --git a/src/types.F b/src/types.F index 26ce25d..f2684ec 100644 --- a/src/types.F +++ b/src/types.F @@ -6,6 +6,7 @@ integer, parameter :: t_VMC = 7 integer, parameter :: t_DMC = 8 integer, parameter :: t_SRMC = 9 + integer, parameter :: t_FKMC = 10 integer, parameter :: t_Simple = 11 integer, parameter :: t_None = 12 @@ -26,7 +27,7 @@ 'VMC ', & 'DMC ', & 'SRMC ', & - ' ', & + 'FKMC ', & 'Simple ', & 'None ', & ' ', &