LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-12-21 20:03:31 +01:00
Code Releases Activity

Merge branch 'master' of github.com:scemama/qmcchem

Browse Source
This commit is contained in:
Anthony Scemama 2017-03-26 21:37:20 +02:00
commit b3ca85eb2d
22 changed files with 988 additions and 120 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ezfio_config/qmc.config
View File

@ -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

5
install/build.ninja
View File

@ -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

1
install/scripts/install_zmq.sh
View File

@ -3,7 +3,6 @@
TARGET=zmq
function _install()
{
# LIBVERSION=4
LIBVERSION=5
cd .. ; QMCCHEM_PATH="$PWD" ; cd -
set +u

106
ocaml/Input.ml
View File

@ -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 =

38
ocaml/Qmcchem_edit.ml
View File

@ -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 )

9
ocaml/Qmcchem_result.ml
View 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

72
ocaml/Random_variable.ml
View File

@ -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) *)

271
promela/qmcchem.pml Normal file
View File

@ -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
}

31
src/PROPERTIES/properties_energy.irp.f
View File

@ -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

16
src/PROPERTIES/properties_general.irp.f
View File

@ -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

6
src/SAMPLING/block.irp.f
View File

@ -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

5
src/SAMPLING/dmc_step.irp.f
View File

@ -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

4
src/SAMPLING/fkmc_step.irp.f
View File

@ -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

380
src/SAMPLING/pdmc_step.irp.f Normal file
View 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(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

5
src/SAMPLING/srmc_step.irp.f
View File

@ -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

4
src/SAMPLING/vmc_step.irp.f
View File

@ -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

7
src/det.irp.f
View File

@ -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 ]

1
src/ezfio_interface.irp.f
View File

@ -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 " , "" ),

28
src/mo.irp.f
View File

@ -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

57
src/pseudo.irp.f
View File

@ -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

12
src/simulation.irp.f
View File

@ -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
View File

@ -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 ', &
' '/)