LCPQ
/
qmcchem
mirror of https://gitlab.com/scemama/qmcchem.git
10
1
Fork 1
Code Releases Activity

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

This commit is contained in:
Anthony Scemama 2017-10-20 17:15:55 +02:00
parent 7806f6df53 7883be32c1
commit c10217157c
45 changed files with 546 additions and 189 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
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.3.tar.gz"
URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/v4.1.3.tar.gz"
URL_ZMQ ="http://download.zeromq.org/zeromq-4.1.4.tar.gz"
URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/4.1.4.tar.gz"
# Rules
#######
rule download
command = [[ -e ${out} ]] || (wget --no-check-certificate ${url} -O ${out}.tmp -o /dev/null && mv ${out}.tmp ${out})
command = [ -e ${out} ] || (wget --no-check-certificate ${url} -O ${out}.tmp -o /dev/null && mv ${out}.tmp ${out})
description = Downloading ${descr}
rule install

40
install/scripts/install_ocaml.sh
View File

@ -6,9 +6,43 @@ set -e
cd .. ; QMCCHEM_PATH="$PWD" ; cd -
PACKAGES="core cryptokit ocamlfind sexplib" # ppx_sexp_conv"
declare -i i
i=$(gcc -dumpversion | cut -d '.' -f 2)
if [[ i -lt 6 ]]
# return 0 if program version is equal or greater than check version
check_version () {
if [[ $1 == $2 ]]
then
return 0
fi
local IFS=.
local i ver1=($1) ver2=($2)
# fill empty fields in ver1 with zeros
for ((i=${#ver1[@]}; i<${#ver2[@]}; i++))
do
ver1[i]=0
done
for ((i=0; i<${#ver1[@]}; i++))
do
if [[ -z ${ver2[i]} ]]
then
# fill empty fields in ver2 with zeros
ver2[i]=0
fi
if ((10#${ver1[i]} > 10#${ver2[i]}))
then
return 1
fi
if ((10#${ver1[i]} < 10#${ver2[i]}))
then
return 0
fi
done
return 0
}
i=$(gcc -dumpversion)
check_version 4.6 $i
if [[ $? == 1 ]]
then
echo "GCC version $(gcc -dumpversion) too old. GCC >= 4.6 required."
exit 1

2
install/scripts/install_ocaml_zmq.sh
View File

@ -32,7 +32,7 @@ export C_INCLUDE_PATH="${QMCCHEM_PATH}/lib":$C_INCLUDE_PATH
export LIBRARY_PATH="${QMCCHEM_PATH}/lib":$LIBRARY_PATH
export LD_LIBRARY_PATH="${QMCCHEM_PATH}/lib":$LD_LIBRARY_PATH
set -u
opam install zmq
opam install zmq conf-zmq
rm -f _build/ocaml_zmq.log
exit 0

10
install/scripts/install_zmq.sh
View File

@ -3,7 +3,7 @@
TARGET=zmq
function _install()
{
LIBVERSION=4
LIBVERSION=5
cd .. ; QMCCHEM_PATH="$PWD" ; cd -
set +u
export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./
@ -14,10 +14,10 @@ function _install()
make -j 8
cd -
rm -f -- "${QMCCHEM_PATH}"/lib/libzmq.{a,so,so.$LIBVERSION}
# cp "${BUILD}"/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/
# cp "${BUILD}"/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION
cp "${BUILD}"/src/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/
cp "${BUILD}"/src/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION
cp "${BUILD}"/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/
cp "${BUILD}"/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION
# cp "${BUILD}"/src/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/
# cp "${BUILD}"/src/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION
cp "${BUILD}"/include/{zmq,zmq_utils}.h "${QMCCHEM_PATH}"/lib/
cd "${QMCCHEM_PATH}"/lib
ln libzmq.so.$LIBVERSION libzmq.so || cp libzmq.so.$LIBVERSION libzmq.so

4
ocaml/Block.ml
View File

@ -1,5 +1,5 @@
open Core.Std;;
open Qptypes;;
open Core
open Qptypes
type t =
{ property : Property.t ;

2
ocaml/Default.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core
let simulation_nucl_fitcusp_factor = lazy(

2
ocaml/Input.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
open Qptypes
open Qputils

2
ocaml/Launcher.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core
type t =
| Srun

2
ocaml/Md5.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
(** Directory containing the list of input files. The directory is created is inexistant. *)
let input_directory = lazy (

2
ocaml/Message.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
open Qptypes
type t =

2
ocaml/Qmcchem_config.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core
(** QMC=Chem installation directory *)

7
ocaml/Qmcchem_dataserver.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
open Qptypes
(** Data server of QMC=Chem.
@ -36,7 +36,7 @@ let run ?(daemon=true) ezfio_filename =
begin
Printf.printf "Generating initial walkers...\n%!";
Unix.fork_exec ~prog:(Lazy.force Qmcchem_config.qmc_create_walkers)
~args:["qmc_create_walkers" ; ezfio_filename] ()
~argv:["qmc_create_walkers" ; ezfio_filename] ()
|> Unix.waitpid_exn ;
Printf.printf "Initial walkers ready\n%!"
end ;
@ -98,8 +98,7 @@ let run ?(daemon=true) ezfio_filename =
let result =
try
ZMQ.Socket.bind socket address;
ZMQ.Socket.unbind socket address;
accu;
accu
with
| _ -> false;
in

4
ocaml/Qmcchem_debug.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let run ~t ezfio_filename=
@ -49,7 +49,7 @@ let run ~t ezfio_filename=
in
tot_size := Byte_units.create `Bytes ((Byte_units.bytes !tot_size) +. (Byte_units.bytes bytes));
Printf.printf "%s\n%!" (Byte_units.to_string !tot_size);
Time.pause (Time.Span.of_float 1.)
Time.pause (Time.Span.of_sec 1.)
done
end
else

2
ocaml/Qmcchem_edit.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let file_header filename = Printf.sprintf
"

12
ocaml/Qmcchem_forwarder.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core
let bind_socket ~socket_type ~socket ~address =
let rec loop = function
@ -11,7 +11,7 @@ let bind_socket ~socket_type ~socket ~address =
ZMQ.Socket.bind socket address;
loop (-1)
with
| Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) )
| Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_sec 1. ; loop (i-1) )
| other_exception -> raise other_exception
in loop 10
@ -40,7 +40,7 @@ let run ezfio_filename dataserver =
in
(* Build qmc executable command *)
let prog, args =
let prog, argv =
qmc,
[ qmc ; ezfio_filename ;
Printf.sprintf "ipc://%s:%d" Qmcchem_config.dev_shm port ];
@ -57,7 +57,7 @@ let run ezfio_filename dataserver =
| Unix.Unix_error _ ->
begin
Unix.chdir tmpdir;
Time.pause @@ Time.Span.of_float 0.1;
Time.pause @@ Time.Span.of_sec 0.1;
match (Sys.file_exists "PID") with
| `No
| `Unknown -> ()
@ -75,7 +75,7 @@ let run ezfio_filename dataserver =
begin
match Signal.send (Signal.of_system_int 0) (`Pid (Pid.of_int pid)) with
| `No_such_process -> ()
| _ -> ignore @@ Unix.exec ~prog ~args ()
| _ -> ignore @@ Unix.exec ~prog ~argv ()
end
end
in
@ -89,7 +89,7 @@ let run ezfio_filename dataserver =
(* Fork a qmc *)
ignore @@
Watchdog.fork_exec ~prog ~args ();
Watchdog.fork_exec ~prog ~argv ();
(* If there are MICs, use them here (TODO) *)

6
ocaml/Qmcchem_info.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let run ezfio_filename =
@ -6,12 +6,12 @@ let run ezfio_filename =
let qmcchem_info =
Lazy.force Qmcchem_config.qmcchem_info
in
let prog, args =
let prog, argv =
qmcchem_info,
[ qmcchem_info ; ezfio_filename ]
in
ignore @@
Unix.exec ~prog ~args ()
Unix.exec ~prog ~argv ()
let spec =

2
ocaml/Qmcchem_md5.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let run ?c ?d ~l ~update ezfio_filename =

4
ocaml/Qmcchem_result.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
open Qptypes
(** Display a table that can be plotted by gnuplot *)
@ -65,7 +65,7 @@ let display_cumulants ~range property =
Printf.printf "Variance = %16.10f\n" cum.(1);
Printf.printf "Centered k3 = %16.10f\n" cum.(2);
Printf.printf "Centered k4 = %16.10f\n" cum.(3);
print_newline ();
Printf.printf "\n%!";
let n = 1. /. 12. *. cum.(2) *. cum.(2) +.
1. /. 48. *. cum.(3) *. cum.(3)
in

18
ocaml/Qmcchem_run.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let full_run ?(start_dataserver=true) ezfio_filename =
(* Identify the job scheduler *)
@ -36,13 +36,13 @@ let full_run ?(start_dataserver=true) ezfio_filename =
(* Start the data server *)
let prog, args =
let prog, argv =
qmcchem, [ qmcchem; "run" ; "-d" ; ezfio_filename]
in
let pid_dataserver =
Watchdog.fork_exec ~prog ~args ()
Watchdog.fork_exec ~prog ~argv ()
in
Printf.printf "%7d : %s\n%!" (Pid.to_int pid_dataserver) (String.concat ~sep:" " args)
Printf.printf "%7d : %s\n%!" (Pid.to_int pid_dataserver) (String.concat ~sep:" " argv)
end;
@ -83,7 +83,7 @@ let full_run ?(start_dataserver=true) ezfio_filename =
| n ->
if (not (test_open_rep_socket ())) then
begin
Time.pause (Time.Span.of_float 0.5);
Time.pause (Time.Span.of_sec 0.5);
count (n-1);
end
else
@ -94,7 +94,7 @@ let full_run ?(start_dataserver=true) ezfio_filename =
(* Start the qmc processes *)
let prog, args =
let prog, argv =
let launcher =
Launcher.(find () |> to_string)
in
@ -110,12 +110,12 @@ let full_run ?(start_dataserver=true) ezfio_filename =
in
let pid_qmc =
try
Watchdog.fork_exec ~prog ~args ()
Watchdog.fork_exec ~prog ~argv ()
with
| Unix.Unix_error _ ->
begin
let command =
String.concat ~sep:" " args
String.concat ~sep:" " argv
in
Printf.printf "
============================================================
@ -126,7 +126,7 @@ Error: Unable to run the following command
Watchdog.kill ()
end
in
Printf.printf "%7d : %s\n%!" (Pid.to_int pid_qmc) (String.concat ~sep:" " args);
Printf.printf "%7d : %s\n%!" (Pid.to_int pid_qmc) (String.concat ~sep:" " argv);
(* Wait for processes to finish *)
Watchdog.join ()

2
ocaml/Qmcchem_stop.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let run ezfio_filename =

2
ocaml/Qputils.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let split_re =
Str.regexp " +"

4
ocaml/Random_variable.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
open Qptypes
type t =
@ -64,7 +64,7 @@ end = struct
(x -. mu) *. ( x -. mu) /. sigma2
in
let pi =
acos (-1.)
Float.acos (-1.)
in
let c =
1. /. (sqrt (sigma2 *. (pi +. pi)))

4
ocaml/Sample.ml
View File

@ -1,9 +1,9 @@
open Core.Std
open Core
type t =
| One_dimensional of float
| Multidimensional of (float array * int)
with sexp
[@ deriving sexp]
let dimension = function
| One_dimensional _ -> 1

4
ocaml/Sample.mli
View File

@ -1,4 +1,6 @@
type t with sexp
open Core
type t [@@ deriving sexp]
val to_float : ?idx:int -> t -> float
val to_float_array : t -> float array
val of_float : float -> t

2
ocaml/Scheduler.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core
type t =
| SGE

6
ocaml/Watchdog.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core
let _list = ref [] ;;
let _running = ref false;;
@ -90,9 +90,9 @@ let del pid =
;;
(** Fork and exec a new process *)
let fork_exec ~prog ~args () =
let fork_exec ~prog ~argv () =
let pid =
Unix.fork_exec ~prog ~args ()
Unix.fork_exec ~prog ~argv ()
in
let f () =

4
ocaml/build.ninja
View File

@ -1,7 +1,7 @@
MAIN=qmcchem
# Main program to build
PACKAGES=-package core,cryptokit,str,ZMQ,sexplib.syntax
PACKAGES=-package core,cryptokit,str,ZMQ
#,ppx_sexp_conv
# Required opam packages, for example:
# PACKAGES=-package core,sexplib.syntax
@ -10,7 +10,7 @@ THREAD=-thread
# If you need threding support, use:
# THREAD=-thread
SYNTAX=-syntax camlp4o
SYNTAX=
# If you need pre-processing, use:
# SYNTAX=-syntax camlp4o

2
ocaml/ninja_ocaml.py
View File

@ -196,7 +196,7 @@ MAIN=
PACKAGES=
# Required opam packages, for example:
# PACKAGES=-package core,sexplib.syntax
# PACKAGES=-package core
THREAD=
# If you need threding support, use:

2
ocaml/qmcchem.ml
View File

@ -1,4 +1,4 @@
open Core.Std
open Core
let command =

12
ocaml/qptypes_generator.ml
View File

@ -1,4 +1,4 @@
open Core.Std;;
open Core
let input_data = "
* Positive_float : float
@ -156,12 +156,12 @@ let untouched = "
let template = format_of_string "
module %s : sig
type t with sexp
type t [@@ deriving sexp]
val to_%s : t -> %s
val of_%s : %s %s -> t
val to_string : t -> string
end = struct
type t = %s with sexp
type t = %s [@@ deriving sexp]
let to_%s x = x
let of_%s %s x = ( %s x )
let to_string x = %s.to_string x
@ -199,13 +199,13 @@ let parse_input input=
let ezfio_template = format_of_string "
module %s : sig
type t with sexp
type t [@@ deriving sexp]
val to_%s : t -> %s
val get_max : unit -> %s
val of_%s : ?min:%s -> ?max:%s -> %s -> t
val to_string : t -> string
end = struct
type t = %s with sexp
type t = %s [@@ deriving sexp]
let to_string x = %s.to_string x
let get_max () =
if (Ezfio.has_%s ()) then
@ -312,7 +312,7 @@ match msg with " ] @
let () =
let input =
String.concat ~sep:"\n"
[ "open Core.Std\nlet warning = print_string\n\n" ;
[ "open Core\nlet warning = print_string\n\n" ;
parse_input input_data ;
parse_input_ezfio input_ezfio ;
create_ezfio_handler ();

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
}

4
scripts/create_properties_ezfio.py
View File

@ -67,7 +67,7 @@ file = open(tmp_filename,'w')
# ----
print >>file, """
(* File generated by ${QMCCHEM_PATH}/src/create_properties.py. Do not
(* File generated by ${QMCCHEM_PATH}/scripts/create_properties.py. Do not
modify here
*)
@ -125,7 +125,7 @@ for p in properties:
print >>file, """;;
let of_string s =
match (String.lowercase s) with
match (String.lowercase_ascii s) with
| "cpu" -> Cpu
| "wall" -> Wall
| "accep" -> Accep"""

2
src/AO/ao_axis.irp.f
View File

@ -45,7 +45,7 @@ subroutine pow_l(r,a,x1,x2,x3)
x3 = 0.
return
end select
end function
end
BEGIN_PROVIDER [ real, ao_axis_block, (ao_block_num_8) ]

2
src/JASTROW/jastrow_core.irp.f
View File

@ -14,7 +14,7 @@
jast_elec_Core_range(i) = 0.d0
else
double precision :: rc
double precision, parameter :: thresh = 0.5 ! function = thresh at rc
double precision, parameter :: thresh = 0.5d0 ! function = thresh at rc
rc = min(0.8d0,max(4.0d0/nucl_charge(i), 0.25d0))
jast_elec_Core_expo(i) = -1.d0/rc**2 * log(thresh)
jast_elec_Core_range(i) = dsqrt(15.d0/jast_elec_Core_expo(i))

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

@ -270,21 +270,14 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
BEGIN_DOC
! Zero-variance parameter on E_loc
END_DOC
E_loc_zv = E_loc + (E_trial-E_loc) * dmc_zv_weight
E_loc_zv = E_loc
E_loc_zv += (E_trial-E_loc) * dmc_zv_weight
! E_loc_zv += - time_step*(E_trial**2 + 1.44341217940434 - E_loc**2)*dmc_zv_weight
! E_loc_zv(3) = dmc_zv_weight_half
! E_loc_zv(:) = 0.d0
END_PROVIDER
BEGIN_PROVIDER [ double precision, E_loc_zv_diag ]
implicit none
BEGIN_DOC
! Zero-variance parameter on E_loc
END_DOC
E_loc_zv_diag = E_trial
END_PROVIDER

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

24
src/SAMPLING/pdmc_step.irp.f
View File

@ -109,9 +109,9 @@ END_SHELL
endif
integer :: info
double precision :: H(0:pdmc_n_diag/2,0:pdmc_n_diag/2), S(0:pdmc_n_diag/2,0:pdmc_n_diag/2), w(0:pdmc_n_diag/2), work(3*pdmc_n_diag+1)
H = 0.d0
S = 0.d0
! double precision :: H(0:pdmc_n_diag/2,0:pdmc_n_diag/2), S(0:pdmc_n_diag/2,0:pdmc_n_diag/2), w(0:pdmc_n_diag/2), work(3*pdmc_n_diag+1)
! H = 0.d0
! S = 0.d0
do while (loop)
@ -234,13 +234,13 @@ END_SHELL
block_weight += pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk)
pdmc_pop_weight_mult(0) = 1.d0/pdmc_weight(i_walk)
do k=0,pdmc_n_diag/2
do l=0,pdmc_n_diag/2
H(k,l) += E_loc*pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
S(k,l) += pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
enddo
enddo
H = H + (E_trial - E_loc)
! do k=0,pdmc_n_diag/2
! do l=0,pdmc_n_diag/2
! H(k,l) += E_loc*pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
! S(k,l) += pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
! enddo
! enddo
! H = H + (E_trial - E_loc)
! else
! pdmc_weight(i_walk) = 1.d0
@ -280,11 +280,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(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

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

@ -90,8 +90,8 @@ for p in properties:
print t.replace("$X",p[1])
END_SHELL
logical :: loop
integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max
logical :: loop
integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max
loop = .True.
call system_clock(cpu0, count_rate, count_max)
@ -320,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(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

49
src/TOOLS/invert.irp.f
View File

@ -1,3 +1,45 @@
! This file contains the fast inversion routines of QMC=Chem for
! small matrices. It may be downloaded here:
! https://raw.githubusercontent.com/scemama/qmcchem/master/src/TOOLS/invert.irp.f
!
! To use it in your Fortran code, you will need to~:
! 1) rename it inverse.f90
! 2) replace all $IRP_ALIGN occurences by
! a) 16 for SSE4.2
! b) 32 for AVX or AVX2
! c) 64 for AVX-512
!
!
! GPL license :
! -------------
!
! QMC=Chem : Quantum Monte Carlo for Chemistry
! Copyright (C) 2009 Anthony SCEMAMA
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along
! with this program; if not, write to the Free Software Foundation, Inc.,
! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
!
! Anthony Scemama
! LCPQ - IRSAMC - CNRS
! Universite Paul Sabatier
! 118, route de Narbonne
! 31062 Toulouse Cedex 4
! scemama@irsamc.ups-tlse.fr
!
!
subroutine invert(a,LDA,na,det_l)
implicit none
double precision, intent(inout) :: a (LDA,na)
@ -43,7 +85,6 @@ subroutine invert_general(a,LDA,na,det_l)
integer :: ipiv(LDA)
!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
integer :: lwork
double precision :: f
integer :: i, j
call dgetrf(na, na, a, LDA, ipiv, inf )
det_l = 1.d0
@ -74,7 +115,6 @@ subroutine sinvert(a,LDA,na,det_l)
integer :: ipiv(LDA)
!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
integer :: lwork
real :: f
integer :: i, j
call sgetrf(na, na, a, LDA, ipiv, inf )
det_l = 1.d0
@ -113,7 +153,6 @@ subroutine invert2(a,LDA,na,det_l)
double precision :: det_l
double precision :: b(2,2)
double precision :: f
b(1,1) = a(1,1)
b(2,1) = a(2,1)
b(1,2) = a(1,2)
@ -134,7 +173,6 @@ subroutine invert3(a,LDA,na,det_l)
double precision :: b(4,3)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
integer :: i
double precision :: f
det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) &
-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) &
+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
@ -166,7 +204,6 @@ subroutine invert4(a,LDA,na,det_l)
double precision :: b(4,4)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
integer :: i,j
double precision :: f
det_l = a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
@ -217,7 +254,6 @@ subroutine invert5(a,LDA,na,det_l)
double precision :: b(5,5)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
integer :: i,j
double precision :: f
det_l = a(1,1)*(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))- &
a(2,3)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)- &
@ -410,7 +446,6 @@ subroutine invert_update(a,LDA,na,det_l,b)
integer :: ipiv(LDA)
!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
integer :: lwork
double precision :: f
integer :: i, j
double precision :: bold(LDA,na)
double precision :: ba(LDA,na)

130
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 ]
@ -1401,15 +1408,10 @@ END_PROVIDER
endif
det_alpha_value(det_i) = det_alpha_value_curr
do i=1,elec_alpha_num
!DIR$ VECTOR ALIGNED
do k=1,4
det_alpha_grad_lapl(k,i,det_i) = det_alpha_grad_lapl_curr(k,i)
enddo
if (do_pseudo) then
det_alpha_pseudo(i,det_i) = det_alpha_pseudo_curr(i)
endif
enddo
det_alpha_grad_lapl(:,:,det_i) = det_alpha_grad_lapl_curr(:,:)
if (do_pseudo) then
det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:)
endif
enddo
@ -1453,15 +1455,10 @@ END_PROVIDER
endif
det_beta_value(det_j) = det_beta_value_curr
!DIR$ LOOP COUNT (200)
do i=elec_alpha_num+1,elec_num
do k=1,4
det_beta_grad_lapl(k,i,det_j) = det_beta_grad_lapl_curr(k,i)
enddo
if (do_pseudo) then
det_beta_pseudo(i,det_j) = det_beta_pseudo_curr(i)
endif
enddo
det_beta_grad_lapl(:,:,det_j) = det_beta_grad_lapl_curr(:,:)
if (do_pseudo) then
det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:)
endif
enddo
@ -1536,53 +1533,64 @@ END_PROVIDER
DaC = 0.d0
CDb = 0.d0
det_num4 = iand(det_num,not(3))
!DIR$ VECTOR ALIGNED
do k=1,det_num4,4
i1 = det_coef_matrix_rows(k )
i2 = det_coef_matrix_rows(k+1)
i3 = det_coef_matrix_rows(k+2)
i4 = det_coef_matrix_rows(k+3)
j1 = det_coef_matrix_columns(k )
j2 = det_coef_matrix_columns(k+1)
j3 = det_coef_matrix_columns(k+2)
j4 = det_coef_matrix_columns(k+3)
if ( (j1 == j2).and.(j1 == j3).and.(j1 == j4) ) then
f = det_beta_value (j1)
CDb(i1) = CDb(i1) + det_coef_matrix_values(k )*f
CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f
if (det_num < ishft(det_alpha_num*det_beta_num,-2)) then
if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then
DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+ det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
+ det_coef_matrix_values(k+2)*det_alpha_value(i1+2) &
+ det_coef_matrix_values(k+3)*det_alpha_value(i1+3)
det_num4 = iand(det_num,not(3))
!DIR$ VECTOR ALIGNED
do k=1,det_num4,4
i1 = det_coef_matrix_rows(k )
i2 = det_coef_matrix_rows(k+1)
i3 = det_coef_matrix_rows(k+2)
i4 = det_coef_matrix_rows(k+3)
j1 = det_coef_matrix_columns(k )
j2 = det_coef_matrix_columns(k+1)
j3 = det_coef_matrix_columns(k+2)
j4 = det_coef_matrix_columns(k+3)
if ( (j1 == j2).and.(j1 == j3).and.(j1 == j4) ) then
f = det_beta_value (j1)
CDb(i1) = CDb(i1) + det_coef_matrix_values(k )*f
CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f
if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then
DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+ det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
+ det_coef_matrix_values(k+2)*det_alpha_value(i1+2) &
+ det_coef_matrix_values(k+3)*det_alpha_value(i1+3)
else
DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+ det_coef_matrix_values(k+1)*det_alpha_value(i2) &
+ det_coef_matrix_values(k+2)*det_alpha_value(i3) &
+ det_coef_matrix_values(k+3)*det_alpha_value(i4)
endif
else
DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+ det_coef_matrix_values(k+1)*det_alpha_value(i2) &
+ det_coef_matrix_values(k+2)*det_alpha_value(i3) &
+ det_coef_matrix_values(k+3)*det_alpha_value(i4)
DaC(j1) = DaC(j1) + det_coef_matrix_values(k )*det_alpha_value(i1)
DaC(j2) = DaC(j2) + det_coef_matrix_values(k+1)*det_alpha_value(i2)
DaC(j3) = DaC(j3) + det_coef_matrix_values(k+2)*det_alpha_value(i3)
DaC(j4) = DaC(j4) + det_coef_matrix_values(k+3)*det_alpha_value(i4)
CDb(i1) = CDb(i1) + det_coef_matrix_values(k )*det_beta_value (j1)
CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*det_beta_value (j2)
CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*det_beta_value (j3)
CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
endif
else
DaC(j1) = DaC(j1) + det_coef_matrix_values(k )*det_alpha_value(i1)
DaC(j2) = DaC(j2) + det_coef_matrix_values(k+1)*det_alpha_value(i2)
DaC(j3) = DaC(j3) + det_coef_matrix_values(k+2)*det_alpha_value(i3)
DaC(j4) = DaC(j4) + det_coef_matrix_values(k+3)*det_alpha_value(i4)
CDb(i1) = CDb(i1) + det_coef_matrix_values(k )*det_beta_value (j1)
CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*det_beta_value (j2)
CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*det_beta_value (j3)
CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
endif
enddo
enddo
do k=det_num4+1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
enddo
do k=det_num4+1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
enddo
else
call dgemv('T',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
size(det_coef_matrix_dense,1), det_alpha_value, 1, 0.d0, DaC, 1)
call dgemv('N',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1)
endif
! Value
! -----

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

2
src/simulation.irp.f
View File

@ -264,7 +264,7 @@ END_PROVIDER
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)
call rinfo(irp_here,'nucl_fitcusp_factor',nucl_fitcusp_factor)
END_PROVIDER

15
src/wf.irp.f
View File

@ -80,6 +80,21 @@ END_PROVIDER
deallocate(buffer)
END_PROVIDER
BEGIN_PROVIDER [ double precision, det_coef_matrix_dense, (det_alpha_num, det_beta_num) ]
implicit none
BEGIN_DOC
! Dense version of det_coef_matrix
END_DOC
integer :: i,j,k
det_coef_matrix_dense = 0.d0
do k=1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
det_coef_matrix_dense(i,j) = det_coef_matrix_values(k)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, det_num ]
implicit none
BEGIN_DOC