LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-17 10:45:28 +02:00
Code Issues Releases Wiki Activity

Merge branch 'scemama-master'

Browse Source
This commit is contained in:
Emmanuel Giner 2017-03-27 15:15:56 +02:00
commit e6505dfb98
91 changed files with 6335 additions and 3780 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
README.md
View File

@ -82,11 +82,11 @@ If you have set the `--developement` flag you can go in any module directory and
### 4) Compiling the OCaml
make -C ocaml
make -C $QP_ROOT/ocaml
### 5) Testing if all is ok
cd tests ; bats bats/qp.bats
cd tests ; ./run_tests.sh
@ -137,10 +137,6 @@ interface: ezfio
#FAQ
### Opam error: cryptokit
You need to install `gmp-dev`.
### Error: ezfio_* is already defined.
#### Why ?
@ -166,5 +162,5 @@ It's caused when we call the DGEMM routine of LAPACK.
##### Fix
Set `ulimit -s unlimited`, before runing `qp_run`. It seem to fix the problem.
Set `ulimit -s unlimited`, before runing `qp_run`. It seems to fix the problem.

2
configure vendored
View File

@ -102,7 +102,7 @@ curl = Info(
default_path=join(QP_ROOT_BIN, "curl"))
zlib = Info(
url='http://www.zlib.net/zlib-1.2.11.tar.gz',
url='http://www.zlib.net/fossils/zlib-1.2.10.tar.gz',
description=' zlib',
default_path=join(QP_ROOT_LIB, "libz.a"))

14
ocaml/Basis.ml
View File

@ -36,9 +36,11 @@ let read_element in_channel at_number element =
let to_string_general ~fmt ~atom_sep b =
let to_string_general ~fmt ~atom_sep ?ele_array b =
let new_nucleus n =
Printf.sprintf "Atom %d" n
match ele_array with
| None -> Printf.sprintf "Atom %d" n
| Some x -> Printf.sprintf "%s" (Element.to_string x.(n-1))
in
let rec do_work accu current_nucleus = function
| [] -> List.rev accu
@ -56,12 +58,12 @@ let to_string_general ~fmt ~atom_sep b =
do_work [new_nucleus 1] 1 b
|> String.concat ~sep:"\n"
let to_string_gamess =
to_string_general ~fmt:Gto.Gamess ~atom_sep:""
let to_string_gamess ?ele_array =
to_string_general ?ele_array ~fmt:Gto.Gamess ~atom_sep:""
let to_string_gaussian b =
let to_string_gaussian ?ele_array b =
String.concat ~sep:"\n"
[ to_string_general ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ]
[ to_string_general ?ele_array ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ]
let to_string ?(fmt=Gto.Gamess) =
match fmt with

2
ocaml/Basis.mli
View File

@ -14,7 +14,7 @@ val read_element :
in_channel -> Nucl_number.t -> Element.t -> (Gto.t * Nucl_number.t) list
(** Convert the basis to a string *)
val to_string : ?fmt:Gto.fmt -> (Gto.t * Nucl_number.t) list -> string
val to_string : ?fmt:Gto.fmt -> ?ele_array:Element.t array -> (Gto.t * Nucl_number.t) list -> string
(** Convert the basis to an MD5 hash *)
val to_md5 : (Gto.t * Nucl_number.t) list -> MD5.t

1
ocaml/Makefile
View File

@ -13,6 +13,7 @@ LIBS=
PKGS=
OCAMLCFLAGS="-g -warn-error A"
OCAMLBUILD=ocamlbuild -j 0 -syntax camlp4o -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS)
MLLFILES=$(wildcard *.mll)
MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml
MLIFILES=$(wildcard *.mli) git
ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml))

175
ocaml/Message.ml
View File

@ -110,7 +110,7 @@ module Disconnect_msg : sig
{ client_id: Id.Client.t ;
state: State.t ;
}
val create : state:string -> client_id:string -> t
val create : state:string -> client_id:int -> t
val to_string : t -> string
end = struct
type t =
@ -118,7 +118,7 @@ end = struct
state: State.t ;
}
let create ~state ~client_id =
{ client_id = Id.Client.of_string client_id ; state = State.of_string state }
{ client_id = Id.Client.of_int client_id ; state = State.of_string state }
let to_string x =
Printf.sprintf "disconnect %s %d"
(State.to_string x.state)
@ -150,18 +150,18 @@ end
module AddTask_msg : sig
type t =
{ state: State.t;
task: string;
tasks: string list;
}
val create : state:string -> task:string -> t
val create : state:string -> tasks:string list -> t
val to_string : t -> string
end = struct
type t =
{ state: State.t;
task: string;
tasks: string list;
}
let create ~state ~task = { state = State.of_string state ; task }
let create ~state ~tasks = { state = State.of_string state ; tasks }
let to_string x =
Printf.sprintf "add_task %s %s" (State.to_string x.state) x.task
Printf.sprintf "add_task %s %s" (State.to_string x.state) (String.concat ~sep:"|" x.tasks)
end
@ -182,44 +182,44 @@ end
module DelTask_msg : sig
type t =
{ state: State.t;
task_id: Id.Task.t
task_ids: Id.Task.t list
}
val create : state:string -> task_id:string -> t
val create : state:string -> task_ids:int list -> t
val to_string : t -> string
end = struct
type t =
{ state: State.t;
task_id: Id.Task.t
task_ids: Id.Task.t list
}
let create ~state ~task_id =
let create ~state ~task_ids =
{ state = State.of_string state ;
task_id = Id.Task.of_string task_id
task_ids = List.map ~f:Id.Task.of_int task_ids
}
let to_string x =
Printf.sprintf "del_task %s %d"
Printf.sprintf "del_task %s %s"
(State.to_string x.state)
(Id.Task.to_int x.task_id)
(String.concat ~sep:"|" @@ List.map ~f:Id.Task.to_string x.task_ids)
end
(** DelTaskReply : Reply to the DelTask message *)
module DelTaskReply_msg : sig
type t
val create : task_id:Id.Task.t -> more:bool -> t
val create : task_ids:Id.Task.t list -> more:bool -> t
val to_string : t -> string
end = struct
type t = {
task_id : Id.Task.t ;
task_ids : Id.Task.t list;
more : bool;
}
let create ~task_id ~more = { task_id ; more }
let create ~task_ids ~more = { task_ids ; more }
let to_string x =
let more =
if x.more then "more"
else "done"
in
Printf.sprintf "del_task_reply %s %d"
more (Id.Task.to_int x.task_id)
Printf.sprintf "del_task_reply %s %s"
more (String.concat ~sep:"|" @@ List.map ~f:Id.Task.to_string x.task_ids)
end
@ -230,7 +230,7 @@ module GetTask_msg : sig
{ client_id: Id.Client.t ;
state: State.t ;
}
val create : state:string -> client_id:string -> t
val create : state:string -> client_id:int -> t
val to_string : t -> string
end = struct
type t =
@ -238,7 +238,7 @@ end = struct
state: State.t ;
}
let create ~state ~client_id =
{ client_id = Id.Client.of_string client_id ; state = State.of_string state }
{ client_id = Id.Client.of_int client_id ; state = State.of_string state }
let to_string x =
Printf.sprintf "get_task %s %d"
(State.to_string x.state)
@ -269,14 +269,14 @@ module GetPsi_msg : sig
type t =
{ client_id: Id.Client.t ;
}
val create : client_id:string -> t
val create : client_id:int -> t
val to_string : t -> string
end = struct
type t =
{ client_id: Id.Client.t ;
}
let create ~client_id =
{ client_id = Id.Client.of_string client_id }
{ client_id = Id.Client.of_int client_id }
let to_string x =
Printf.sprintf "get_psi %d"
(Id.Client.to_int x.client_id)
@ -365,14 +365,14 @@ module PutPsi_msg : sig
n_det_selectors : Strictly_positive_int.t option;
psi : Psi.t option }
val create :
client_id:string ->
n_state:string ->
n_det:string ->
psi_det_size:string ->
client_id:int ->
n_state:int ->
n_det:int ->
psi_det_size:int ->
psi_det:string option ->
psi_coef:string option ->
n_det_generators: string option ->
n_det_selectors:string option ->
n_det_generators: int option ->
n_det_selectors:int option ->
energy:string option -> t
val to_string_list : t -> string list
val to_string : t -> string
@ -388,20 +388,17 @@ end = struct
let create ~client_id ~n_state ~n_det ~psi_det_size ~psi_det ~psi_coef
~n_det_generators ~n_det_selectors ~energy =
let n_state, n_det, psi_det_size =
Int.of_string n_state
|> Strictly_positive_int.of_int ,
Int.of_string n_det
|> Strictly_positive_int.of_int ,
Int.of_string psi_det_size
|> Strictly_positive_int.of_int
Strictly_positive_int.of_int n_state,
Strictly_positive_int.of_int n_det,
Strictly_positive_int.of_int psi_det_size
in
assert (Strictly_positive_int.to_int psi_det_size >=
Strictly_positive_int.to_int n_det);
let n_det_generators, n_det_selectors =
match n_det_generators, n_det_selectors with
| Some x, Some y ->
Some (Strictly_positive_int.of_int @@ Int.of_string x),
Some (Strictly_positive_int.of_int @@ Int.of_string y)
Some (Strictly_positive_int.of_int x),
Some (Strictly_positive_int.of_int y)
| _ -> None, None
in
let psi =
@ -411,7 +408,7 @@ end = struct
~psi_coef ~n_det_generators ~n_det_selectors ~energy)
| _ -> None
in
{ client_id = Id.Client.of_string client_id ;
{ client_id = Id.Client.of_int client_id ;
n_state ; n_det ; psi_det_size ; n_det_generators ;
n_det_selectors ; psi }
@ -463,48 +460,48 @@ module TaskDone_msg : sig
type t =
{ client_id: Id.Client.t ;
state: State.t ;
task_id: Id.Task.t ;
task_ids: Id.Task.t list ;
}
val create : state:string -> client_id:string -> task_id:string -> t
val create : state:string -> client_id:int -> task_ids:int list -> t
val to_string : t -> string
end = struct
type t =
{ client_id: Id.Client.t ;
state: State.t ;
task_id: Id.Task.t;
task_ids: Id.Task.t list;
}
let create ~state ~client_id ~task_id =
{ client_id = Id.Client.of_string client_id ;
let create ~state ~client_id ~task_ids =
{ client_id = Id.Client.of_int client_id ;
state = State.of_string state ;
task_id = Id.Task.of_string task_id;
task_ids = List.map ~f:Id.Task.of_int task_ids;
}
let to_string x =
Printf.sprintf "task_done %s %d %d"
Printf.sprintf "task_done %s %d %s"
(State.to_string x.state)
(Id.Client.to_int x.client_id)
(Id.Task.to_int x.task_id)
(String.concat ~sep:"|" @@ List.map ~f:Id.Task.to_string x.task_ids)
end
(** Terminate *)
module Terminate_msg : sig
type t
val create : unit -> t
val create : t
val to_string : t -> string
end = struct
type t = Terminate
let create () = Terminate
let create = Terminate
let to_string x = "terminate"
end
(** OK *)
module Ok_msg : sig
type t
val create : unit -> t
val create : t
val to_string : t -> string
end = struct
type t = Ok
let create () = Ok
let create = Ok
let to_string x = "ok"
end
@ -551,45 +548,45 @@ type t =
let of_string s =
let l =
String.split ~on:' ' s
|> List.filter ~f:(fun x -> (String.strip x) <> "")
|> List.map ~f:String.lowercase
in
match l with
| "add_task" :: state :: task ->
AddTask (AddTask_msg.create ~state ~task:(String.concat ~sep:" " task) )
| "del_task" :: state :: task_id :: [] ->
DelTask (DelTask_msg.create ~state ~task_id)
| "get_task" :: state :: client_id :: [] ->
GetTask (GetTask_msg.create ~state ~client_id)
| "task_done" :: state :: client_id :: task_id :: [] ->
TaskDone (TaskDone_msg.create ~state ~client_id ~task_id)
| "disconnect" :: state :: client_id :: [] ->
Disconnect (Disconnect_msg.create ~state ~client_id)
| "connect" :: t :: [] ->
Connect (Connect_msg.create t)
| "new_job" :: state :: push_address_tcp :: push_address_inproc :: [] ->
Newjob (Newjob_msg.create push_address_tcp push_address_inproc state)
| "end_job" :: state :: [] ->
Endjob (Endjob_msg.create state)
| "terminate" :: [] ->
Terminate (Terminate_msg.create () )
| "get_psi" :: client_id :: [] ->
GetPsi (GetPsi_msg.create ~client_id)
| "put_psi" :: client_id :: n_state :: n_det :: psi_det_size :: n_det_generators :: n_det_selectors :: [] ->
PutPsi (PutPsi_msg.create ~client_id ~n_state ~n_det ~psi_det_size
~n_det_generators:(Some n_det_generators) ~n_det_selectors:(Some n_det_selectors)
~psi_det:None ~psi_coef:None ~energy:None )
| "put_psi" :: client_id :: n_state :: n_det :: psi_det_size :: [] ->
PutPsi (PutPsi_msg.create ~client_id ~n_state ~n_det ~psi_det_size ~n_det_generators:None
~n_det_selectors:None ~psi_det:None ~psi_coef:None ~energy:None)
| "ok" :: [] -> Ok (Ok_msg.create ())
| "error" :: rest -> Error (Error_msg.create (String.concat ~sep:" " rest))
| "set_stopped" :: [] -> SetStopped
| "set_running" :: [] -> SetRunning
| "set_waiting" :: [] -> SetWaiting
| _ -> failwith "Message not understood"
let open Message_lexer in
match parse s with
| AddTask_ { state ; tasks } ->
AddTask (AddTask_msg.create ~state ~tasks)
| DelTask_ { state ; task_ids } ->
DelTask (DelTask_msg.create ~state ~task_ids)
| GetTask_ { state ; client_id } ->
GetTask (GetTask_msg.create ~state ~client_id)
| TaskDone_ { state ; task_ids ; client_id } ->
TaskDone (TaskDone_msg.create ~state ~client_id ~task_ids)
| Disconnect_ { state ; client_id } ->
Disconnect (Disconnect_msg.create ~state ~client_id)
| Connect_ socket ->
Connect (Connect_msg.create socket)
| NewJob_ { state ; push_address_tcp ; push_address_inproc } ->
Newjob (Newjob_msg.create push_address_tcp push_address_inproc state)
| EndJob_ state ->
Endjob (Endjob_msg.create state)
| GetPsi_ client_id ->
GetPsi (GetPsi_msg.create ~client_id)
| PutPsi_ { client_id ; n_state ; n_det ; psi_det_size ; n_det_generators ; n_det_selectors } ->
begin
match n_det_selectors, n_det_generators with
| Some s, Some g ->
PutPsi (PutPsi_msg.create ~client_id ~n_state ~n_det ~psi_det_size
~n_det_generators:(Some g) ~n_det_selectors:(Some s)
~psi_det:None ~psi_coef:None ~energy:None )
| _ ->
PutPsi (PutPsi_msg.create ~client_id ~n_state ~n_det ~psi_det_size
~n_det_generators:None ~n_det_selectors:None
~psi_det:None ~psi_coef:None ~energy:None )
end
| Terminate_ -> Terminate (Terminate_msg.create )
| SetWaiting_ -> SetWaiting
| SetStopped_ -> SetStopped
| SetRunning_ -> SetRunning
| Ok_ -> Ok (Ok_msg.create)
| Error_ m -> Error (Error_msg.create m)
let to_string = function

265
ocaml/Message_lexer.mll Normal file
View File

@ -0,0 +1,265 @@
{
type kw_type =
| TEXT of string
| WORD of string
| INTEGER of int
| FLOAT of float
| NONE
| ADD_TASK
| DEL_TASK
| GET_TASK
| TASK_DONE
| DISCONNECT
| CONNECT
| NEW_JOB
| END_JOB
| TERMINATE
| GET_PSI
| PUT_PSI
| OK
| ERROR
| SET_STOPPED
| SET_RUNNING
| SET_WAITING
type state_tasks = { state : string ; tasks : string list ; }
type state_taskids = { state : string ; task_ids : int list ; }
type state_taskids_clientid = { state : string ; task_ids : int list ; client_id : int ; }
type state_clientid = { state : string ; client_id : int ; }
type state_tcp_inproc = { state : string ; push_address_tcp : string ; push_address_inproc : string ; }
type psi = { client_id: int ; n_state: int ; n_det: int ; psi_det_size: int ;
n_det_generators: int option ; n_det_selectors: int option }
type msg =
| AddTask_ of state_tasks
| DelTask_ of state_taskids
| GetTask_ of state_clientid
| TaskDone_ of state_taskids_clientid
| Disconnect_ of state_clientid
| Connect_ of string
| NewJob_ of state_tcp_inproc
| EndJob_ of string
| Terminate_
| GetPsi_ of int
| PutPsi_ of psi
| Ok_
| Error_ of string
| SetStopped_
| SetRunning_
| SetWaiting_
}
let word = [^' ' '\t' '\n']+
let text = [^ ' ' '|']+[^ '|']+
let integer = ['0'-'9']+
let real = '-'? integer '.' integer (['e' 'E'] '-'? integer)?
let white = [' ' '\t']+
rule get_text = parse
| text as t { TEXT t }
| eof { TERMINATE }
| _ { NONE }
and get_int = parse
| integer as i { INTEGER (int_of_string i) }
| eof { TERMINATE }
| _ { NONE }
and get_word = parse
| word as w { WORD w }
| eof { TERMINATE }
| _ { NONE }
and kw = parse
| "add_task" { ADD_TASK }
| "del_task" { DEL_TASK }
| "get_task" { GET_TASK }
| "task_done" { TASK_DONE }
| "disconnect" { DISCONNECT }
| "connect" { CONNECT }
| "new_job" { NEW_JOB }
| "end_job" { END_JOB }
| "terminate" { TERMINATE }
| "get_psi" { GET_PSI }
| "put_psi" { PUT_PSI }
| "ok" { OK }
| "error" { ERROR }
| "set_stopped" { SET_STOPPED }
| "set_running" { SET_RUNNING }
| "set_waiting" { SET_WAITING }
| _ { NONE }
{
let rec read_text ?(accu=[]) lexbuf =
let token =
get_text lexbuf
in
match token with
| TEXT t -> read_text ~accu:(t::accu) lexbuf
| TERMINATE -> List.rev accu
| NONE -> read_text ~accu lexbuf
| _ -> failwith "Error in MessageLexer (2)"
and read_word lexbuf =
let token =
get_word lexbuf
in
match token with
| WORD w -> w
| NONE -> read_word lexbuf
| _ -> failwith "Error in MessageLexer (3)"
and read_int lexbuf =
let token =
get_int lexbuf
in
match token with
| INTEGER i -> i
| NONE -> read_int lexbuf
| _ -> failwith "Error in MessageLexer (4)"
and read_ints ?(accu=[]) lexbuf =
let token =
get_int lexbuf
in
match token with
| INTEGER i -> read_ints ~accu:(i::accu) lexbuf
| TERMINATE -> List.rev accu
| NONE -> read_ints ~accu lexbuf
| _ -> failwith "Error in MessageLexer (4)"
and parse_rec lexbuf =
let token =
kw lexbuf
in
match token with
| ADD_TASK ->
let state = read_word lexbuf in
let tasks = read_text lexbuf in
AddTask_ { state ; tasks }
| DEL_TASK ->
let state = read_word lexbuf in
let task_ids = read_ints lexbuf in
DelTask_ { state ; task_ids }
| GET_TASK ->
let state = read_word lexbuf in
let client_id = read_int lexbuf in
GetTask_ { state ; client_id }
| TASK_DONE ->
let state = read_word lexbuf in
let client_id = read_int lexbuf in
let task_ids = read_ints lexbuf in
TaskDone_ { state ; task_ids ; client_id }
| DISCONNECT ->
let state = read_word lexbuf in
let client_id = read_int lexbuf in
Disconnect_ { state ; client_id }
| GET_PSI ->
let client_id = read_int lexbuf in
GetPsi_ client_id
| PUT_PSI ->
let client_id = read_int lexbuf in
let n_state = read_int lexbuf in
let n_det = read_int lexbuf in
let psi_det_size = read_int lexbuf in
let n_det_generators, n_det_selectors =
try
(Some (read_int lexbuf), Some (read_int lexbuf))
with (Failure _) -> (None, None)
in
PutPsi_ { client_id ; n_state ; n_det ; psi_det_size ; n_det_generators ; n_det_selectors }
| CONNECT ->
let socket = read_word lexbuf in
Connect_ socket
| NEW_JOB ->
let state = read_word lexbuf in
let push_address_tcp = read_word lexbuf in
let push_address_inproc = read_word lexbuf in
NewJob_ { state ; push_address_tcp ; push_address_inproc }
| END_JOB ->
let state = read_word lexbuf in
EndJob_ state
| ERROR ->
let message = List.hd (read_text lexbuf) in
Error_ message
| OK -> Ok_
| SET_WAITING -> SetWaiting_
| SET_RUNNING -> SetRunning_
| SET_STOPPED -> SetStopped_
| TERMINATE -> Terminate_
| NONE -> parse_rec lexbuf
| _ -> failwith "Error in MessageLexer"
let parse message =
let lexbuf =
Lexing.from_string message
in
parse_rec lexbuf
let debug () =
let l = [
"add_task state_pouet Task pouet zob" ;
"add_task state_pouet Task pouet zob |Task2 zob | Task3 prout" ;
"del_task state_pouet 12345" ;
"del_task state_pouet 12345 | 6789 | 10 | 11" ;
"get_task state_pouet 12" ;
"task_done state_pouet 12 12345";
"task_done state_pouet 12 12345 | 678 | 91011";
"connect tcp";
"disconnect state_pouet 12";
"new_job state_pouet tcp://test.com:12345 ipc:///dev/shm/x.socket";
"end_job state_pouet";
"terminate" ;
"set_running" ;
"set_stopped" ;
"set_waiting" ;
"ok" ;
"error my_error" ;
"get_psi 12" ;
"put_psi 12 2 1000 10000 800 900" ;
"put_psi 12 2 1000 10000"
]
|> List.map parse
in
List.map (function
| AddTask_ { state ; tasks } -> Printf.sprintf "ADD_TASK state:\"%s\" tasks:{\"%s\"}" state (String.concat "\"}|{\"" tasks)
| DelTask_ { state ; task_ids } -> Printf.sprintf "DEL_TASK state:\"%s\" task_ids:{%s}" state (String.concat "|" @@ List.map string_of_int task_ids)
| GetTask_ { state ; client_id } -> Printf.sprintf "GET_TASK state:\"%s\" task_id:%d" state client_id
| TaskDone_ { state ; task_ids ; client_id } -> Printf.sprintf "TASK_DONE state:\"%s\" task_ids:{%s} client_id:%d" state (String.concat "|" @@ List.map string_of_int task_ids) client_id
| Disconnect_ { state ; client_id } -> Printf.sprintf "DISCONNECT state:\"%s\" client_id:%d" state client_id
| Connect_ socket -> Printf.sprintf "CONNECT socket:\"%s\"" socket
| NewJob_ { state ; push_address_tcp ; push_address_inproc } -> Printf.sprintf "NEW_JOB state:\"%s\" tcp:\"%s\" inproc:\"%s\"" state push_address_tcp push_address_inproc
| EndJob_ state -> Printf.sprintf "END_JOB state:\"%s\"" state
| GetPsi_ client_id -> Printf.sprintf "GET_PSI client_id:%d" client_id
| PutPsi_ { client_id ; n_state ; n_det ; psi_det_size ; n_det_generators ; n_det_selectors } ->
begin
match n_det_selectors, n_det_generators with
| Some s, Some g -> Printf.sprintf "PUT_PSI client_id:%d n_state:%d n_det:%d psi_det_size:%d n_det_generators:%d n_det_selectors:%d" client_id n_state n_det psi_det_size g s
| _ -> Printf.sprintf "PUT_PSI client_id:%d n_state:%d n_det:%d psi_det_size:%d" client_id n_state n_det psi_det_size
end
| Terminate_ -> "TERMINATE"
| SetWaiting_ -> "SET_WAITING"
| SetStopped_ -> "SET_STOPPED"
| SetRunning_ -> "SET_RUNNING"
| Ok_ -> "OK"
| Error_ s -> Printf.sprintf "ERROR: \"%s\"" s
) l
|> List.iter print_endline
}

103
ocaml/TaskServer.ml
View File

@ -62,7 +62,15 @@ let bind_socket ~socket_type ~socket ~port =
| Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) )
| other_exception -> raise other_exception
in loop 60;
ZMQ.Socket.bind socket @@ Printf.sprintf "ipc:///tmp/qp_run:%d" port
let filename =
Printf.sprintf "/tmp/qp_run:%d" port
in
begin
match Sys.file_exists filename with
| `Yes -> Sys.remove filename
| _ -> ()
end;
ZMQ.Socket.bind socket ("ipc://"^filename)
let hostname = lazy (
@ -99,7 +107,7 @@ let ip_address = lazy (
let reply_ok rep_socket =
Message.Ok_msg.create ()
Message.Ok_msg.create
|> Message.Ok_msg.to_string
|> ZMQ.Socket.send rep_socket
@ -121,7 +129,7 @@ let stop ~port =
ZMQ.Socket.set_linger_period req_socket 1_000_000;
ZMQ.Socket.connect req_socket address;
Message.Terminate (Message.Terminate_msg.create ())
Message.Terminate (Message.Terminate_msg.create)
|> Message.to_string
|> ZMQ.Socket.send req_socket ;
@ -289,9 +297,9 @@ let disconnect msg program_state rep_socket =
let del_task msg program_state rep_socket =
let state, task_id =
let state, task_ids =
msg.Message.DelTask_msg.state,
msg.Message.DelTask_msg.task_id
msg.Message.DelTask_msg.task_ids
in
let failure () =
@ -302,13 +310,14 @@ let del_task msg program_state rep_socket =
let new_program_state =
{ program_state with
queue = Queuing_system.del_task ~task_id program_state.queue
queue = List.fold ~f:(fun queue task_id -> Queuing_system.del_task ~task_id queue)
~init:program_state.queue task_ids
}
in
let more =
(Queuing_system.number_of_tasks new_program_state.queue > 0)
in
Message.DelTaskReply (Message.DelTaskReply_msg.create ~task_id ~more)
Message.DelTaskReply (Message.DelTaskReply_msg.create ~task_ids ~more)
|> Message.to_string
|> ZMQ.Socket.send ~block:true rep_socket ; (** /!\ Has to be blocking *)
new_program_state
@ -329,9 +338,9 @@ let del_task msg program_state rep_socket =
let add_task msg program_state rep_socket =
let state, task =
let state, tasks =
msg.Message.AddTask_msg.state,
msg.Message.AddTask_msg.task
msg.Message.AddTask_msg.tasks
in
let increment_progress_bar = function
@ -339,59 +348,17 @@ let add_task msg program_state rep_socket =
| None -> None
in
let rec add_task_triangle program_state imax = function
| 0 -> program_state
| i ->
let task =
Printf.sprintf "%d %d" i imax
in
let new_program_state =
{ program_state with
queue = Queuing_system.add_task ~task program_state.queue ;
progress_bar = increment_progress_bar program_state.progress_bar ;
}
in
add_task_triangle new_program_state imax (i-1)
in
let rec add_task_range program_state i = function
| j when (j < i) -> program_state
| j ->
let task =
Printf.sprintf "%d" j
in
let new_program_state =
{ program_state with
queue = Queuing_system.add_task ~task program_state.queue ;
progress_bar = increment_progress_bar program_state.progress_bar ;
}
in
add_task_range new_program_state i (j-1)
in
let new_program_state = function
| "triangle" :: i_str :: [] ->
let imax =
Int.of_string i_str
in
add_task_triangle program_state imax imax
| "range" :: i_str :: j_str :: [] ->
let i, j =
Int.of_string i_str,
Int.of_string j_str
in
add_task_range program_state i j
| _ ->
{ program_state with
queue = Queuing_system.add_task ~task program_state.queue ;
progress_bar = increment_progress_bar program_state.progress_bar ;
}
in
let result =
String.split ~on:' ' task
|> List.filter ~f:(fun x -> x <> "")
|> new_program_state
let new_queue, new_bar =
List.fold ~f:(fun (queue, bar) task ->
Queuing_system.add_task ~task queue,
increment_progress_bar bar)
~init:(program_state.queue, program_state.progress_bar) tasks
in
{ program_state with
queue = new_queue;
progress_bar = new_bar
}
in
reply_ok rep_socket;
result
@ -448,10 +415,10 @@ let get_task msg program_state rep_socket pair_socket =
let task_done msg program_state rep_socket =
let state, client_id, task_id =
let state, client_id, task_ids =
msg.Message.TaskDone_msg.state,
msg.Message.TaskDone_msg.client_id,
msg.Message.TaskDone_msg.task_id
msg.Message.TaskDone_msg.task_ids
in
let increment_progress_bar = function
@ -464,10 +431,16 @@ let task_done msg program_state rep_socket =
program_state
and success () =
let new_queue, new_bar =
List.fold ~f:(fun (queue, bar) task_id ->
Queuing_system.end_task ~task_id ~client_id queue,
increment_progress_bar bar)
~init:(program_state.queue, program_state.progress_bar) task_ids
in
let result =
{ program_state with
queue = Queuing_system.end_task ~task_id ~client_id program_state.queue ;
progress_bar = increment_progress_bar program_state.progress_bar ;
queue = new_queue;
progress_bar = new_bar
}
in
reply_ok rep_socket;

8
ocaml/qptypes_generator.ml
View File

@ -42,8 +42,8 @@ let input_data = "
* Det_number_max : int
assert (x > 0) ;
if (x > 100000000) then
warning \"More than 100 million determinants\";
if (x > 10000000000) then
warning \"More than 10 billion determinants\";
* States_number : int
assert (x > 0) ;
@ -140,8 +140,8 @@ let input_ezfio = "
* Det_number : int
determinants_n_det
1 : 100000000
More than 100 million of determinants
1 : 10000000000
More than 10 billion of determinants
"
;;

9
plugins/CAS_SD_ZMQ/EZFIO.cfg
View File

@ -1,10 +1,15 @@
[energy]
type: double precision
doc: "Calculated CAS-SD energy"
doc: Calculated CAS-SD energy
interface: ezfio
[energy_pt2]
type: double precision
doc: "Calculated selected CAS-SD energy with PT2 correction"
doc: Calculated selected CAS-SD energy with PT2 correction
interface: ezfio
[do_ddci]
type: logical
doc: If true, remove purely inactive double excitations
interface: ezfio,provider,ocaml
default: False

121
plugins/CAS_SD_ZMQ/cassd_zmq.irp.f
View File

@ -132,124 +132,3 @@ program fci_zmq
call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before(1)+pt2(1))
end
subroutine ZMQ_selection(N_in, pt2)
use f77_zmq
use selection_types
implicit none
character*(512) :: task
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer, intent(in) :: N_in
type(selection_buffer) :: b
integer :: i, N
integer, external :: omp_get_thread_num
double precision, intent(out) :: pt2(N_states)
if (.True.) then
PROVIDE pt2_e0_denominator
N = max(N_in,1)
provide nproc
call new_parallel_job(zmq_to_qp_run_socket,"selection")
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
call zmq_set_running(zmq_to_qp_run_socket)
call create_selection_buffer(N, N*2, b)
endif
integer :: i_generator, i_generator_start, i_generator_max, step
! step = int(max(1.,10*elec_num/mo_tot_num)
step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num ))
step = max(1,step)
do i= 1, N_det_generators,step
i_generator_start = i
i_generator_max = min(i+step-1,N_det_generators)
write(task,*) i_generator_start, i_generator_max, 1, N
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
end do
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call selection_collector(b, pt2)
else
call selection_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'selection')
if (N_in > 0) then
call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN
call copy_H_apply_buffer_to_wf()
if (s2_eig) then
call make_s2_eigenfunction
endif
endif
end subroutine
subroutine selection_slave_inproc(i)
implicit none
integer, intent(in) :: i
call run_selection_slave(1,i,pt2_e0_denominator)
end
subroutine selection_collector(b, pt2)
use f77_zmq
use selection_types
use bitmasks
implicit none
type(selection_buffer), intent(inout) :: b
double precision, intent(out) :: pt2(N_states)
double precision :: pt2_mwen(N_states)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_pull_socket
integer(ZMQ_PTR) :: zmq_socket_pull
integer :: msg_size, rc, more
integer :: acc, i, j, robin, N, ntask
double precision, allocatable :: val(:)
integer(bit_kind), allocatable :: det(:,:,:)
integer, allocatable :: task_id(:)
integer :: done
real :: time, time0
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket()
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det))
done = 0
more = 1
pt2(:) = 0d0
call CPU_TIME(time0)
do while (more == 1)
call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask)
pt2 += pt2_mwen
do i=1, N
call add_to_selection_buffer(b, det(1,1,i), val(i))
end do
do i=1, ntask
if(task_id(i) == 0) then
print *, "Error in collector"
endif
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
end do
done += ntask
call CPU_TIME(time)
! print *, "DONE" , done, time - time0
end do
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull)
call sort_selection_buffer(b)
end subroutine

4
plugins/CAS_SD_ZMQ/ezfio_interface.irp.f
View File

@ -1,4 +0,0 @@
! DO NOT MODIFY BY HAND
! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
! from file /home/scemama/quantum_package/src/CAS_SD_ZMQ/EZFIO.cfg

6
plugins/CAS_SD_ZMQ/run_selection_slave.irp.f
View File

@ -50,8 +50,6 @@ subroutine run_selection_slave(thread,iproc,energy)
else
if(N /= buf%N) stop "N changed... wtf man??"
end if
!print *, "psi_selectors_coef ", psi_selectors_coef(N_det_selectors-5:N_det_selectors, 1)
!call debug_det(psi_selectors(1,1,N_det_selectors), N_int)
do i_generator=i_generator_start,i_generator_max,step
call select_connected(i_generator,energy,pt2,buf)
enddo
@ -115,7 +113,7 @@ subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ
! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
end subroutine
@ -149,7 +147,7 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, nt
if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP
! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
end subroutine

2539
plugins/CAS_SD_ZMQ/selection.irp.f
View File

File diff suppressed because it is too large Load Diff

109
plugins/CAS_SD_ZMQ/target_pt2_ratio_cassd.irp.f Normal file
View File

@ -0,0 +1,109 @@
program fci_zmq
implicit none
integer :: i,j,k
logical, external :: detEq
double precision, allocatable :: pt2(:)
integer :: Nmin, Nmax
integer :: n_det_before, to_select
double precision :: threshold_davidson_in, ratio, E_ref
double precision, allocatable :: psi_coef_ref(:,:)
integer(bit_kind), allocatable :: psi_det_ref(:,:,:)
allocate (pt2(N_states))
pt2 = 1.d0
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
! Stopping criterion is the PT2max
double precision :: E_CI_before(N_states)
do while (dabs(pt2(1)) > pt2_max)
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k=1, N_states
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k)
print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k)
enddo
print *, '-----'
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_cas_sd_zmq_energy(CI_energy(1))
n_det_before = N_det
to_select = N_det
to_select = max(64-to_select, to_select)
call ZMQ_selection(to_select, pt2)
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
call diagonalize_CI
call save_wavefunction
call ezfio_set_cas_sd_zmq_energy(CI_energy(1))
enddo
threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
threshold_generators = max(threshold_generators,threshold_generators_pt2)
threshold_davidson = threshold_davidson_in
TOUCH threshold_selectors threshold_generators threshold_davidson
call diagonalize_CI
call ZMQ_selection(0, pt2)
E_ref = CI_energy(1) + pt2(1)
print *, 'Est FCI = ', E_ref
Nmax = N_det
Nmin = 2
allocate (psi_coef_ref(size(psi_coef_sorted,1),size(psi_coef_sorted,2)))
allocate (psi_det_ref(N_int,2,size(psi_det_sorted,3)))
psi_coef_ref = psi_coef_sorted
psi_det_ref = psi_det_sorted
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
TOUCH psi_coef psi_det
do while (Nmax-Nmin > 1)
psi_coef = psi_coef_ref
psi_det = psi_det_ref
TOUCH psi_det psi_coef
call diagonalize_CI
ratio = (CI_energy(1) - HF_energy) / (E_ref - HF_energy)
if (ratio < var_pt2_ratio) then
Nmin = N_det
else
Nmax = N_det
psi_coef_ref = psi_coef
psi_det_ref = psi_det
TOUCH psi_det psi_coef
endif
N_det = Nmin + (Nmax-Nmin)/2
print *, '-----'
print *, 'Det min, Det max: ', Nmin, Nmax
print *, 'Ratio : ', ratio, ' ~ ', var_pt2_ratio
print *, 'N_det = ', N_det
print *, 'E = ', CI_energy(1)
call save_wavefunction
enddo
call ZMQ_selection(0, pt2)
print *, '------'
print *, 'HF_energy = ', HF_energy
print *, 'Est FCI = ', E_ref
print *, 'E = ', CI_energy(1)
print *, 'PT2 = ', pt2(1)
print *, 'E+PT2 = ', CI_energy(1)+pt2(1)
E_CI_before(1:N_states) = CI_energy(1:N_states)
call save_wavefunction
call ezfio_set_cas_sd_zmq_energy(CI_energy(1))
call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before(1)+pt2(1))
end

16
plugins/Full_CI_ZMQ/energy.irp.f
View File

@ -1,11 +1,23 @@
BEGIN_PROVIDER [ logical, initialize_pt2_E0_denominator ]
implicit none
BEGIN_DOC
! If true, initialize pt2_E0_denominator
END_DOC
initialize_pt2_E0_denominator = .True.
END_PROVIDER
BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ]
implicit none
BEGIN_DOC
! E0 in the denominator of the PT2
END_DOC
pt2_E0_denominator(1:N_states) = CI_electronic_energy(1:N_states)
if (initialize_pt2_E0_denominator) then
pt2_E0_denominator(1:N_states) = CI_electronic_energy(1:N_states)
! pt2_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
! pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
call write_double(6,pt2_E0_denominator(1)+nuclear_repulsion, 'PT2 Energy denominator')
call write_double(6,pt2_E0_denominator(1)+nuclear_repulsion, 'PT2 Energy denominator')
else
pt2_E0_denominator = -huge(1.d0)
endif
END_PROVIDER

137
plugins/Full_CI_ZMQ/fci_zmq.irp.f
View File

@ -68,8 +68,8 @@ program fci_zmq
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
n_det_before = N_det
to_select = 2*N_det
to_select = max(64-to_select, to_select)
to_select = N_det
to_select = max(N_det, to_select)
to_select = min(to_select, N_det_max-n_det_before)
call ZMQ_selection(to_select, pt2)
@ -96,11 +96,17 @@ program fci_zmq
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
threshold_generators = max(threshold_generators,threshold_generators_pt2)
TOUCH threshold_selectors threshold_generators
!threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
!threshold_generators = max(threshold_generators,threshold_generators_pt2)
!TOUCH threshold_selectors threshold_generators
threshold_selectors = 1.d0
threshold_generators = 1d0
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ZMQ_selection(0, pt2)
double precision :: relative_error
relative_error=1.d-3
pt2 = 0.d0
call ZMQ_pt2(pt2,relative_error)
!call ZMQ_selection(0, pt2)! pour non-stochastic
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
@ -119,122 +125,3 @@ program fci_zmq
end
subroutine ZMQ_selection(N_in, pt2)
use f77_zmq
use selection_types
implicit none
character*(512) :: task
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer, intent(in) :: N_in
type(selection_buffer) :: b
integer :: i, N
integer, external :: omp_get_thread_num
double precision, intent(out) :: pt2(N_states)
if (.True.) then
PROVIDE pt2_e0_denominator
N = max(N_in,1)
provide nproc
call new_parallel_job(zmq_to_qp_run_socket,"selection")
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
call zmq_set_running(zmq_to_qp_run_socket)
call create_selection_buffer(N, N*2, b)
endif
integer :: i_generator, i_generator_start, i_generator_max, step
! step = int(max(1.,10*elec_num/mo_tot_num)
step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num ))
step = max(1,step)
do i= 1, N_det_generators,step
i_generator_start = i
i_generator_max = min(i+step-1,N_det_generators)
write(task,*) i_generator_start, i_generator_max, 1, N
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
end do
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call selection_collector(b, pt2)
else
call selection_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'selection')
if (N_in > 0) then
call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN
call copy_H_apply_buffer_to_wf()
if (s2_eig) then
call make_s2_eigenfunction
endif
endif
end subroutine
subroutine selection_slave_inproc(i)
implicit none
integer, intent(in) :: i
call run_selection_slave(1,i,pt2_e0_denominator)
end
subroutine selection_collector(b, pt2)
use f77_zmq
use selection_types
use bitmasks
implicit none
type(selection_buffer), intent(inout) :: b
double precision, intent(out) :: pt2(N_states)
double precision :: pt2_mwen(N_states)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_pull_socket
integer(ZMQ_PTR) :: zmq_socket_pull
integer :: msg_size, rc, more
integer :: acc, i, j, robin, N, ntask
double precision, allocatable :: val(:)
integer(bit_kind), allocatable :: det(:,:,:)
integer, allocatable :: task_id(:)
integer :: done
real :: time, time0
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket()
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det))
done = 0
more = 1
pt2(:) = 0d0
call CPU_TIME(time0)
do while (more == 1)
call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask)
pt2 += pt2_mwen
do i=1, N
call add_to_selection_buffer(b, det(1,1,i), val(i))
end do
do i=1, ntask
if(task_id(i) == 0) then
print *, "Error in collector"
endif
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
end do
done += ntask
call CPU_TIME(time)
! print *, "DONE" , done, time - time0
end do
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull)
call sort_selection_buffer(b)
end subroutine

70
plugins/Full_CI_ZMQ/pt2_slave.irp.f Normal file
View File

@ -0,0 +1,70 @@
program pt2_slave
implicit none
BEGIN_DOC
! Helper program to compute the PT2 in distributed mode.
END_DOC
read_wf = .False.
SOFT_TOUCH read_wf
call provide_everything
call switch_qp_run_to_master
call run_wf
end
subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
end
subroutine run_wf
use f77_zmq
implicit none
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states_diag)
character*(64) :: states(1)
integer :: rc, i
call provide_everything
zmq_context = f77_zmq_ctx_new ()
states(1) = 'pt2'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_states(states,zmq_state,1)
if(trim(zmq_state) == 'Stopped') then
exit
else if (trim(zmq_state) == 'pt2') then
! Selection
! ---------
print *, 'PT2'
call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call pt2_slave_tcp(i, energy)
!$OMP END PARALLEL
print *, 'PT2 done'
endif
end do
end
subroutine pt2_slave_tcp(i,energy)
implicit none
double precision, intent(in) :: energy(N_states_diag)
integer, intent(in) :: i
logical :: lstop
lstop = .False.
call run_pt2_slave(0,i,energy,lstop)
end

47
plugins/Full_CI_ZMQ/pt2_stoch.irp.f Normal file
View File

@ -0,0 +1,47 @@
program pt2_stoch
implicit none
initialize_pt2_E0_denominator = .False.
read_wf = .True.
SOFT_TOUCH initialize_pt2_E0_denominator read_wf
PROVIDE mo_bielec_integrals_in_map
call run
end
subroutine run
implicit none
integer :: i,j,k
logical, external :: detEq
double precision, allocatable :: pt2(:)
integer :: degree
integer :: n_det_before, to_select
double precision :: threshold_davidson_in
double precision :: E_CI_before(N_states), relative_error
if (.true.) then
call ezfio_get_full_ci_zmq_energy(E_CI_before(1))
pt2_e0_denominator(:) = E_CI_before(1) - nuclear_repulsion
SOFT_TOUCH pt2_e0_denominator read_wf
endif
allocate (pt2(N_states))
pt2 = 0.d0
threshold_selectors = 1.d0
threshold_generators = 1d0
relative_error = 1.d-6
call ZMQ_pt2(pt2, relative_error)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k=1,N_states
print *, 'State', k
print *, 'PT2 = ', pt2
print *, 'E = ', E_CI_before
print *, 'E+PT2 = ', E_CI_before+pt2
print *, '-----'
enddo
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
end

567
plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f Normal file
View File

@ -0,0 +1,567 @@
BEGIN_PROVIDER [ integer, fragment_first ]
implicit none
fragment_first = first_det_of_teeth(1)
END_PROVIDER
subroutine ZMQ_pt2(pt2,relative_error)
use f77_zmq
use selection_types
implicit none
character(len=64000) :: task
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_to_qp_run_socket2
type(selection_buffer) :: b
integer, external :: omp_get_thread_num
double precision, intent(in) :: relative_error
double precision, intent(out) :: pt2(N_states)
double precision, allocatable :: pt2_detail(:,:), comb(:)
logical, allocatable :: computed(:)
integer, allocatable :: tbc(:)
integer :: i, j, Ncomb, generator_per_task, i_generator_end
integer, external :: pt2_find
double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
double precision, external :: omp_get_wtime
double precision :: time0, time
allocate(pt2_detail(N_states, N_det_generators), comb(N_det_generators/2), computed(N_det_generators), tbc(0:size_tbc))
sumabove = 0d0
sum2above = 0d0
Nabove = 0d0
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral
!call random_seed()
computed = .false.
tbc(0) = first_det_of_comb - 1
do i=1, tbc(0)
tbc(i) = i
computed(i) = .true.
end do
pt2_detail = 0d0
time0 = omp_get_wtime()
print *, "time - avg - err - n_combs"
generator_per_task = 1
do while(.true.)
call write_time(6)
call new_parallel_job(zmq_to_qp_run_socket,"pt2")
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
call create_selection_buffer(1, 1*2, b)
Ncomb=size(comb)
call get_carlo_workbatch(computed, comb, Ncomb, tbc)
call write_time(6)
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer :: ipos
logical :: tasks
tasks = .False.
ipos=1
do i=1,tbc(0)
if(tbc(i) > fragment_first) then
write(task(ipos:ipos+20),'(I9,X,I9,''|'')') 0, tbc(i)
ipos += 20
if (ipos > 64000) then
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
ipos=1
tasks = .True.
endif
else
do j=1,fragment_count
write(task(ipos:ipos+20),'(I9,X,I9,''|'')') j, tbc(i)
ipos += 20
if (ipos > 64000) then
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
ipos=1
tasks = .True.
endif
end do
end if
end do
if (ipos > 1) then
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
tasks = .True.
endif
if (tasks) then
call zmq_set_running(zmq_to_qp_run_socket)
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc+1) &
!$OMP PRIVATE(i)
i = omp_get_thread_num()
if (i==0) then
call pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2)
else
call pt2_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'pt2')
else
pt2(1) = sum(pt2_detail(1,:))
endif
tbc(0) = 0
if (pt2(1) /= 0.d0) then
exit
endif
end do
end subroutine
subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, computed, sumabove, sum2above, Nabove)
integer, intent(in) :: tbc(0:size_tbc), Ncomb
logical, intent(in) :: computed(N_det_generators)
double precision, intent(in) :: comb(Ncomb), pt2_detail(N_states, N_det_generators)
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
integer :: i, dets(comb_teeth)
double precision :: myVal, myVal2
mainLoop : do i=1,Ncomb
call get_comb(comb(i), dets, comb_teeth)
do j=1,comb_teeth
if(.not.(computed(dets(j)))) then
exit mainLoop
end if
end do
myVal = 0d0
myVal2 = 0d0
do j=comb_teeth,1,-1
myVal += pt2_detail(1, dets(j)) * pt2_weight_inv(dets(j)) * comb_step
sumabove(j) += myVal
sum2above(j) += myVal*myVal
Nabove(j) += 1
end do
end do mainLoop
end subroutine
subroutine pt2_slave_inproc(i)
implicit none
integer, intent(in) :: i
call run_pt2_slave(1,i,pt2_e0_denominator)
end
subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2)
use f77_zmq
use selection_types
use bitmasks
implicit none
integer, intent(in) :: Ncomb
double precision, intent(inout) :: pt2_detail(N_states, N_det_generators)
double precision, intent(in) :: comb(Ncomb), relative_error
logical, intent(inout) :: computed(N_det_generators)
integer, intent(in) :: tbc(0:size_tbc)
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
double precision, intent(out) :: pt2(N_states)
type(selection_buffer), intent(inout) :: b
double precision, allocatable :: pt2_mwen(:,:)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_pull_socket
integer(ZMQ_PTR) :: zmq_socket_pull
integer :: msg_size, rc, more
integer :: acc, i, j, robin, N, ntask
double precision, allocatable :: val(:)
integer(bit_kind), allocatable :: det(:,:,:)
integer, allocatable :: task_id(:)
integer :: done, Nindex
integer, allocatable :: index(:)
double precision, save :: time0 = -1.d0
double precision :: time, timeLast
double precision, external :: omp_get_wtime
integer :: tooth, firstTBDcomb, orgTBDcomb
integer, allocatable :: parts_to_get(:)
logical, allocatable :: actually_computed(:)
allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators), &
pt2_mwen(N_states, N_det_generators) )
actually_computed(:) = computed(:)
parts_to_get(:) = 1
if(fragment_first > 0) then
parts_to_get(1:fragment_first) = fragment_count
endif
do i=1,tbc(0)
actually_computed(tbc(i)) = .false.
end do
orgTBDcomb = Nabove(1)
firstTBDcomb = 1
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket()
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det_generators), index(1))
more = 1
if (time0 < 0.d0) then
time0 = omp_get_wtime()
endif
timeLast = time0
print *, 'N_deterministic = ', first_det_of_teeth(1)-1
pullLoop : do while (more == 1)
call pull_pt2_results(zmq_socket_pull, Nindex, index, pt2_mwen, task_id, ntask)
do i=1,Nindex
pt2_detail(:, index(i)) += pt2_mwen(:,i)
parts_to_get(index(i)) -= 1
if(parts_to_get(index(i)) < 0) then
print *, i, index(i), parts_to_get(index(i)), Nindex
print *, "PARTS ??"
print *, parts_to_get
stop "PARTS ??"
end if
if(parts_to_get(index(i)) == 0) actually_computed(index(i)) = .true.
end do
do i=1, ntask
if(task_id(i) == 0) then
print *, "Error in collector"
endif
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
end do
time = omp_get_wtime()
if(time - timeLast > 1d1 .or. more /= 1) then
timeLast = time
do i=1, first_det_of_teeth(1)-1
if(.not.(actually_computed(i))) then
print *, "PT2 : deterministic part not finished"
cycle pullLoop
end if
end do
double precision :: E0, avg, eqt, prop
call do_carlo(tbc, Ncomb+1-firstTBDcomb, comb(firstTBDcomb), pt2_detail, actually_computed, sumabove, sum2above, Nabove)
firstTBDcomb = Nabove(1) - orgTBDcomb + 1
if(Nabove(1) < 2d0) cycle
call get_first_tooth(actually_computed, tooth)
done = 0
do i=first_det_of_teeth(tooth), first_det_of_teeth(tooth+1)-1
if(actually_computed(i)) done = done + 1
end do
E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1))
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop
avg = E0 + (sumabove(tooth) / Nabove(tooth))
eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
time = omp_get_wtime()
if (dabs(eqt/avg) < relative_error) then
pt2(1) = avg
! exit pullLoop
endif
print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
end if
end do pullLoop
print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull)
call sort_selection_buffer(b)
end subroutine
integer function pt2_find(v, w, sze, imin, imax)
implicit none
integer, intent(in) :: sze, imin, imax
double precision, intent(in) :: v, w(sze)
integer :: i,l,h
integer, parameter :: block=64
l = imin
h = imax-1
do while(h-l >= block)
i = ishft(h+l,-1)
if(w(i+1) > v) then
h = i-1
else
l = i+1
end if
end do
!DIR$ LOOP COUNT (64)
do pt2_find=l,h
if(w(pt2_find) >= v) then
exit
end if
end do
end function
BEGIN_PROVIDER [ integer, comb_teeth ]
implicit none
comb_teeth = 100
END_PROVIDER
subroutine get_first_tooth(computed, first_teeth)
implicit none
logical, intent(in) :: computed(N_det_generators)
integer, intent(out) :: first_teeth
integer :: i, first_det
first_det = 1
first_teeth = 1
do i=first_det_of_comb, N_det_generators
if(.not.(computed(i))) then
first_det = i
exit
end if
end do
do i=comb_teeth, 1, -1
if(first_det_of_teeth(i) < first_det) then
first_teeth = i
exit
end if
end do
end subroutine
subroutine get_last_full_tooth(computed, last_tooth)
implicit none
logical, intent(in) :: computed(N_det_generators)
integer, intent(out) :: last_tooth
integer :: i, j, missing
last_tooth = 0
combLoop : do i=comb_teeth, 1, -1
missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-12) ! /4096
do j=first_det_of_teeth(i), first_det_of_teeth(i+1)-1
if(.not.computed(j)) then
missing -= 1
if(missing < 0) cycle combLoop
end if
end do
last_tooth = i
exit
end do combLoop
end subroutine
BEGIN_PROVIDER [ integer, size_tbc ]
implicit none
BEGIN_DOC
! Size of the tbc array
END_DOC
size_tbc = N_det_generators + fragment_count*fragment_first
END_PROVIDER
subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
implicit none
double precision, intent(out) :: comb(Ncomb)
integer, intent(inout) :: tbc(0:size_tbc)
integer, intent(inout) :: Ncomb
logical, intent(inout) :: computed(N_det_generators)
integer :: i, j, last_full, dets(comb_teeth), tbc_save
integer :: icount, n
n = tbc(0)
icount = 0
call RANDOM_NUMBER(comb)
do i=1,size(comb)
comb(i) = comb(i) * comb_step
tbc_save = tbc(0)
!DIR$ FORCEINLINE
call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
if (tbc(0) < size(tbc)) then
Ncomb = i
else
tbc(0) = tbc_save
return
endif
icount = icount + tbc(0) - tbc_save
if (icount > n) then
call get_filling_teeth(computed, tbc)
icount = 0
n = ishft(tbc_save,-4)
endif
enddo
end subroutine
subroutine get_filling_teeth(computed, tbc)
implicit none
integer, intent(inout) :: tbc(0:size_tbc)
logical, intent(inout) :: computed(N_det_generators)
integer :: i, j, k, last_full, dets(comb_teeth)
call get_last_full_tooth(computed, last_full)
if(last_full /= 0) then
if (tbc(0) > size(tbc) - first_det_of_teeth(last_full+1) -2) then
return
endif
k = tbc(0)+1
do j=1,first_det_of_teeth(last_full+1)-1
if(.not.(computed(j))) then
tbc(k) = j
k=k+1
computed(j) = .true.
end if
end do
tbc(0) = k-1
end if
end subroutine
subroutine reorder_tbc(tbc)
implicit none
integer, intent(inout) :: tbc(0:size_tbc)
logical, allocatable :: ltbc(:)
integer :: i, ci
allocate(ltbc(size_tbc))
ltbc(:) = .false.
do i=1,tbc(0)
ltbc(tbc(i)) = .true.
end do
ci = 0
do i=1,size_tbc
if(ltbc(i)) then
ci = ci+1
tbc(ci) = i
end if
end do
end subroutine
subroutine get_comb(stato, dets, ct)
implicit none
integer, intent(in) :: ct
double precision, intent(in) :: stato
integer, intent(out) :: dets(ct)
double precision :: curs
integer :: j
integer, external :: pt2_find
curs = 1d0 - stato
do j = comb_teeth, 1, -1
!DIR$ FORCEINLINE
dets(j) = pt2_find(curs, pt2_cweight,size(pt2_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
curs -= comb_step
end do
end subroutine
subroutine add_comb(comb, computed, tbc, stbc, ct)
implicit none
integer, intent(in) :: stbc, ct
double precision, intent(in) :: comb
logical, intent(inout) :: computed(N_det_generators)
integer, intent(inout) :: tbc(0:stbc)
integer :: i, k, l, dets(ct)
!DIR$ FORCEINLINE
call get_comb(comb, dets, ct)
k=tbc(0)+1
do i = 1, ct
l = dets(i)
if(.not.(computed(l))) then
tbc(k) = l
k = k+1
computed(l) = .true.
end if
end do
tbc(0) = k-1
end subroutine
BEGIN_PROVIDER [ double precision, pt2_weight, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, pt2_cweight, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, pt2_cweight_cache, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, comb_step ]
&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ]
&BEGIN_PROVIDER [ integer, first_det_of_comb ]
implicit none
integer :: i
double precision :: norm_left, stato
integer, external :: pt2_find
pt2_weight(1) = psi_coef_generators(1,1)**2
pt2_cweight(1) = psi_coef_generators(1,1)**2
do i=2,N_det_generators
pt2_weight(i) = psi_coef_generators(i,1)**2
pt2_cweight(i) = pt2_cweight(i-1) + psi_coef_generators(i,1)**2
end do
pt2_weight = pt2_weight / pt2_cweight(N_det_generators)
pt2_cweight = pt2_cweight / pt2_cweight(N_det_generators)
norm_left = 1d0
comb_step = 1d0/dfloat(comb_teeth)
do i=1,N_det_generators
if(pt2_weight(i)/norm_left < comb_step*.5d0) then
first_det_of_comb = i
exit
end if
norm_left -= pt2_weight(i)
end do
comb_step = (1d0 - pt2_cweight(first_det_of_comb-1)) * comb_step
stato = 1d0 - comb_step
iloc = N_det_generators
do i=comb_teeth, 1, -1
integer :: iloc
iloc = pt2_find(stato, pt2_cweight, N_det_generators, 1, iloc)
first_det_of_teeth(i) = iloc
stato -= comb_step
end do
first_det_of_teeth(comb_teeth+1) = N_det_generators + 1
first_det_of_teeth(1) = first_det_of_comb
if(first_det_of_teeth(1) /= first_det_of_comb) then
print *, 'Error in ', irp_here
stop "comb provider"
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, pt2_weight_inv, (N_det_generators) ]
implicit none
BEGIN_DOC
! Inverse of pt2_weight array
END_DOC
integer :: i
do i=1,N_det_generators
pt2_weight_inv(i) = 1.d0/pt2_weight(i)
enddo
END_PROVIDER

168
plugins/Full_CI_ZMQ/run_pt2_slave.irp.f Normal file
View File

@ -0,0 +1,168 @@
subroutine run_pt2_slave(thread,iproc,energy)
use f77_zmq
use selection_types
implicit none
double precision, intent(in) :: energy(N_states_diag)
integer, intent(in) :: thread, iproc
integer :: rc, i
integer :: worker_id, task_id(1), ctask, ltask
character*(512) :: task
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_push_socket
integer(ZMQ_PTR) :: zmq_socket_push
type(selection_buffer) :: buf, buf2
logical :: done
double precision :: pt2(N_states)
double precision,allocatable :: pt2_detail(:,:)
integer :: index
integer :: Nindex
allocate(pt2_detail(N_states, N_det))
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_push = new_zmq_push_socket(thread)
call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
if(worker_id == -1) then
print *, "WORKER -1"
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_push_socket(zmq_socket_push,thread)
return
end if
buf%N = 0
ctask = 1
Nindex=1
pt2 = 0d0
pt2_detail = 0d0
do
call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task)
done = task_id(ctask) == 0
if (done) then
ctask = ctask - 1
else
integer :: i_generator, i_i_generator, N, subset
read (task,*) subset, index
!!!!!
N=1
!!!!!
if(buf%N == 0) then
! Only first time
call create_selection_buffer(N, N*2, buf)
call create_selection_buffer(N, N*3, buf2)
else
if(N /= buf%N) stop "N changed... wtf man??"
end if
do i_i_generator=1, Nindex
i_generator = index
call select_connected(i_generator,energy,pt2_detail(1, i_i_generator),buf,subset)
pt2(:) += pt2_detail(:, i_generator)
enddo
endif
if(done .or. ctask == size(task_id)) then
if(buf%N == 0 .and. ctask > 0) stop "uninitialized selection_buffer"
do i=1, ctask
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
end do
if(ctask > 0) then
call push_pt2_results(zmq_socket_push, Nindex, index, pt2_detail, task_id(1), ctask)
do i=1,buf%cur
call add_to_selection_buffer(buf2, buf%det(1,1,i), buf%val(i))
enddo
call sort_selection_buffer(buf2)
buf%mini = buf2%mini
pt2 = 0d0
pt2_detail(:,:Nindex) = 0d0
buf%cur = 0
end if
ctask = 0
end if
if(done) exit
ctask = ctask + 1
end do
call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_push_socket(zmq_socket_push,thread)
end subroutine
subroutine push_pt2_results(zmq_socket_push, N, index, pt2_detail, task_id, ntask)
use f77_zmq
use selection_types
implicit none
integer(ZMQ_PTR), intent(in) :: zmq_socket_push
double precision, intent(in) :: pt2_detail(N_states, N_det)
integer, intent(in) :: ntask, N, index, task_id(*)
integer :: rc
rc = f77_zmq_send( zmq_socket_push, N, 4, ZMQ_SNDMORE)
if(rc /= 4) stop "push"
rc = f77_zmq_send( zmq_socket_push, index, 4, ZMQ_SNDMORE)
if(rc /= 4*N) stop "push"
rc = f77_zmq_send( zmq_socket_push, pt2_detail, 8*N_states*N, ZMQ_SNDMORE)
if(rc /= 8*N_states*N) stop "push"
rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE)
if(rc /= 4) stop "push"
rc = f77_zmq_send( zmq_socket_push, task_id, ntask*4, 0)
if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ
character*(2) :: ok
rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0)
end subroutine
subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntask)
use f77_zmq
use selection_types
implicit none
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
double precision, intent(inout) :: pt2_detail(N_states, N_det)
integer, intent(out) :: index
integer, intent(out) :: N, ntask, task_id(*)
integer :: rc, rn, i
rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0)
if(rc /= 4) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, index, 4, 0)
if(rc /= 4*N) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, pt2_detail, N_states*8*N, 0)
if(rc /= 8*N_states*N) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
if(rc /= 4) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)
if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0)
end subroutine
BEGIN_PROVIDER [ double precision, pt2_workload, (N_det) ]
integer :: i
do i=1,N_det
pt2_workload(:) = dfloat(N_det - i + 1)**2
end do
pt2_workload = pt2_workload / sum(pt2_workload)
END_PROVIDER

14
plugins/Full_CI_ZMQ/run_selection_slave.irp.f
View File

@ -41,8 +41,8 @@ subroutine run_selection_slave(thread,iproc,energy)
if (done) then
ctask = ctask - 1
else
integer :: i_generator, i_generator_start, i_generator_max, step, N
read (task,*) i_generator_start, i_generator_max, step, N
integer :: i_generator, N
read(task,*) i_generator, N
if(buf%N == 0) then
! Only first time
call create_selection_buffer(N, N*2, buf)
@ -50,11 +50,7 @@ subroutine run_selection_slave(thread,iproc,energy)
else
if(N /= buf%N) stop "N changed... wtf man??"
end if
!print *, "psi_selectors_coef ", psi_selectors_coef(N_det_selectors-5:N_det_selectors, 1)
!call debug_det(psi_selectors(1,1,N_det_selectors), N_int)
do i_generator=i_generator_start,i_generator_max,step
call select_connected(i_generator,energy,pt2,buf)
enddo
call select_connected(i_generator,energy,pt2,buf,0)
endif
if(done .or. ctask == size(task_id)) then
@ -115,7 +111,7 @@ subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ
! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
end subroutine
@ -149,7 +145,7 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, nt
if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP
! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
end subroutine

2314
plugins/Full_CI_ZMQ/selection.irp.f
View File

File diff suppressed because it is too large Load Diff

14
plugins/Full_CI_ZMQ/selection_buffer.irp.f
View File

@ -27,7 +27,7 @@ subroutine add_to_selection_buffer(b, det, val)
if(dabs(val) >= b%mini) then
b%cur += 1
b%det(:,:,b%cur) = det(:,:)
b%det(1:N_int,1:2,b%cur) = det(1:N_int,1:2)
b%val(b%cur) = val
if(b%cur == size(b%val)) then
call sort_selection_buffer(b)
@ -57,13 +57,15 @@ subroutine sort_selection_buffer(b)
call dsort(absval, iorder, b%cur)
do i=1, nmwen
detmp(:,:,i) = b%det(:,:,iorder(i))
detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i))
detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i))
vals(i) = b%val(iorder(i))
end do
b%det(:,:,:nmwen) = detmp(:,:,:)
b%det(:,:,nmwen+1:) = 0_bit_kind
b%val(:nmwen) = vals(:)
b%val(nmwen+1:) = 0d0
b%det = 0_bit_kind
b%val = 0d0
b%det(1:N_int,1,1:nmwen) = detmp(1:N_int,1,1:nmwen)
b%det(1:N_int,2,1:nmwen) = detmp(1:N_int,2,1:nmwen)
b%val(1:nmwen) = vals(1:nmwen)
b%mini = max(b%mini,dabs(b%val(b%N)))
b%cur = nmwen
end subroutine

59
plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f
View File

@ -12,8 +12,8 @@ program selection_slave
end
subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context mo_mono_elec_integral
! PROVIDE pt2_e0_denominator mo_tot_num N_int
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count
end
subroutine run_wf
@ -23,16 +23,19 @@ subroutine run_wf
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states)
character*(64) :: states(2)
character*(64) :: states(4)
integer :: rc, i
logical :: force_update
call provide_everything
zmq_context = f77_zmq_ctx_new ()
states(1) = 'selection'
states(2) = 'davidson'
states(3) = 'pt2'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
force_update = .True.
do
@ -52,7 +55,7 @@ subroutine run_wf
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call selection_slave_tcp(i, energy)
call run_selection_slave(0,i,energy)
!$OMP END PARALLEL
print *, 'Selection done'
@ -62,46 +65,34 @@ subroutine run_wf
! --------
print *, 'Davidson'
call davidson_miniserver_get()
call davidson_miniserver_get(force_update)
force_update = .False.
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call davidson_slave_tcp(i)
!$OMP END PARALLEL
print *, 'Davidson done'
else if (trim(zmq_state) == 'pt2') then
! PT2
! ---
print *, 'PT2'
call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
logical :: lstop
lstop = .False.
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call run_pt2_slave(0,i,energy,lstop)
!$OMP END PARALLEL
print *, 'PT2 done'
endif
end do
end
subroutine update_energy(energy)
implicit none
double precision, intent(in) :: energy(N_states)
BEGIN_DOC
! Update energy when it is received from ZMQ
END_DOC
integer :: j,k
do j=1,N_states
do k=1,N_det
CI_eigenvectors(k,j) = psi_coef(k,j)
enddo
enddo
call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int)
if (.True.) then
do k=1,N_states
ci_electronic_energy(k) = energy(k)
enddo
TOUCH ci_electronic_energy CI_eigenvectors_s2 CI_eigenvectors
endif
call write_double(6,ci_energy,'Energy')
end
subroutine selection_slave_tcp(i,energy)
implicit none
double precision, intent(in) :: energy(N_states)
integer, intent(in) :: i
call run_selection_slave(0,i,energy)
end

24
plugins/Full_CI_ZMQ/selection_slave.irp.f
View File

@ -13,7 +13,7 @@ end
subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
PROVIDE pt2_e0_denominator mo_tot_num N_int
PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count
end
subroutine run_wf
@ -60,28 +60,6 @@ subroutine run_wf
end do
end
subroutine update_energy(energy)
implicit none
double precision, intent(in) :: energy(N_states)
BEGIN_DOC
! Update energy when it is received from ZMQ
END_DOC
integer :: j,k
do j=1,N_states
do k=1,N_det
CI_eigenvectors(k,j) = psi_coef(k,j)
enddo
enddo
call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int)
if (.True.) then
do k=1,N_states
ci_electronic_energy(k) = energy(k)
enddo
TOUCH ci_electronic_energy CI_eigenvectors_s2 CI_eigenvectors
endif
call write_double(6,ci_energy,'Energy')
end
subroutine selection_slave_tcp(i,energy)
implicit none

109
plugins/Full_CI_ZMQ/target_pt2_ratio_zmq.irp.f Normal file
View File

@ -0,0 +1,109 @@
program fci_zmq
implicit none
integer :: i,j,k
logical, external :: detEq
double precision, allocatable :: pt2(:)
integer :: Nmin, Nmax
integer :: n_det_before, to_select
double precision :: threshold_davidson_in, ratio, E_ref
double precision, allocatable :: psi_coef_ref(:,:)
integer(bit_kind), allocatable :: psi_det_ref(:,:,:)
allocate (pt2(N_states))
pt2 = 1.d0
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
! Stopping criterion is the PT2max
double precision :: E_CI_before(N_states)
do while (dabs(pt2(1)) > pt2_max)
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k=1, N_states
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k)
print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k)
enddo
print *, '-----'
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
n_det_before = N_det
to_select = N_det
to_select = max(64-to_select, to_select)
call ZMQ_selection(to_select, pt2)
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
call diagonalize_CI
call save_wavefunction
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
enddo
threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
threshold_generators = max(threshold_generators,threshold_generators_pt2)
threshold_davidson = threshold_davidson_in
TOUCH threshold_selectors threshold_generators threshold_davidson
call diagonalize_CI
call ZMQ_selection(0, pt2)
E_ref = CI_energy(1) + pt2(1)
print *, 'Est FCI = ', E_ref
Nmax = N_det
Nmin = 2
allocate (psi_coef_ref(size(psi_coef_sorted,1),size(psi_coef_sorted,2)))
allocate (psi_det_ref(N_int,2,size(psi_det_sorted,3)))
psi_coef_ref = psi_coef_sorted
psi_det_ref = psi_det_sorted
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
TOUCH psi_coef psi_det
do while (Nmax-Nmin > 1)
psi_coef = psi_coef_ref
psi_det = psi_det_ref
TOUCH psi_det psi_coef
call diagonalize_CI
ratio = (CI_energy(1) - HF_energy) / (E_ref - HF_energy)
if (ratio < var_pt2_ratio) then
Nmin = N_det
else
Nmax = N_det
psi_coef_ref = psi_coef
psi_det_ref = psi_det
TOUCH psi_det psi_coef
endif
N_det = Nmin + (Nmax-Nmin)/2
print *, '-----'
print *, 'Det min, Det max: ', Nmin, Nmax
print *, 'Ratio : ', ratio, ' ~ ', var_pt2_ratio
print *, 'N_det = ', N_det
print *, 'E = ', CI_energy(1)
call save_wavefunction
enddo
call ZMQ_selection(0, pt2)
print *, '------'
print *, 'HF_energy = ', HF_energy
print *, 'Est FCI = ', E_ref
print *, 'E = ', CI_energy(1)
print *, 'PT2 = ', pt2(1)
print *, 'E+PT2 = ', CI_energy(1)+pt2(1)
E_CI_before(1:N_states) = CI_energy(1:N_states)
call save_wavefunction
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
end

95
plugins/Full_CI_ZMQ/target_pt2_zmq.irp.f Normal file
View File

@ -0,0 +1,95 @@
program fci_zmq
implicit none
integer :: i,j,k
logical, external :: detEq
double precision, allocatable :: pt2(:)
integer :: Nmin, Nmax
integer :: n_det_before, to_select
double precision :: threshold_davidson_in, ratio, E_ref, pt2_ratio
allocate (pt2(N_states))
pt2 = 1.d0
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
double precision :: E_CI_before(N_states)
do while (dabs(pt2(1)) > pt2_max)
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k=1, N_states
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k)
print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k)
enddo
print *, '-----'
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
n_det_before = N_det
to_select = N_det
to_select = max(64-to_select, to_select)
call ZMQ_selection(to_select, pt2)
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
call diagonalize_CI
call save_wavefunction
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
enddo
threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
threshold_generators = max(threshold_generators,threshold_generators_pt2)
threshold_davidson = threshold_davidson_in
TOUCH threshold_selectors threshold_generators threshold_davidson
call diagonalize_CI
call ZMQ_selection(0, pt2)
E_ref = CI_energy(1) + pt2(1)
pt2_ratio = (E_ref + pt2_max - HF_energy) / (E_ref - HF_energy)
print *, 'Est FCI = ', E_ref
Nmax = N_det
Nmin = N_det/8
do while (Nmax-Nmin > 1)
call diagonalize_CI
ratio = (CI_energy(1) - HF_energy) / (E_ref - HF_energy)
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
TOUCH psi_coef psi_det
if (ratio < pt2_ratio) then
Nmin = N_det
to_select = (Nmax-Nmin)/2
call ZMQ_selection(to_select, pt2)
else
Nmax = N_det
N_det = Nmin + (Nmax-Nmin)/2
endif
print *, '-----'
print *, 'Det min, Det max: ', Nmin, Nmax
print *, 'Ratio : ', ratio, ' ~ ', pt2_ratio
print *, 'HF_energy = ', HF_energy
print *, 'Est FCI = ', E_ref
print *, 'N_det = ', N_det
print *, 'E = ', CI_energy(1)
print *, 'PT2 = ', pt2(1)
enddo
call ZMQ_selection(0, pt2)
print *, '------'
print *, 'E = ', CI_energy(1)
print *, 'PT2 = ', pt2(1)
E_CI_before(1:N_states) = CI_energy(1:N_states)
call save_wavefunction
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
end

114
plugins/Full_CI_ZMQ/zmq_selection.irp.f Normal file
View File

@ -0,0 +1,114 @@
subroutine ZMQ_selection(N_in, pt2)
use f77_zmq
use selection_types
implicit none
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer, intent(in) :: N_in
type(selection_buffer) :: b
integer :: i, N
integer, external :: omp_get_thread_num
double precision, intent(out) :: pt2(N_states)
PROVIDE fragment_count
if (.True.) then
PROVIDE pt2_e0_denominator
N = max(N_in,1)
provide nproc
call new_parallel_job(zmq_to_qp_run_socket,"selection")
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
call zmq_set_running(zmq_to_qp_run_socket)
call create_selection_buffer(N, N*2, b)
endif
character(len=:), allocatable :: task
task = repeat(' ',20*N_det_generators)
do i= 1, N_det_generators
write(task(20*(i-1)+1:20*i),'(I9,X,I9,''|'')') i, N
end do
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call selection_collector(b, pt2)
else
call selection_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'selection')
if (N_in > 0) then
call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0)
call copy_H_apply_buffer_to_wf()
if (s2_eig) then
call make_s2_eigenfunction
endif
call save_wavefunction
endif
end subroutine
subroutine selection_slave_inproc(i)
implicit none
integer, intent(in) :: i
call run_selection_slave(1,i,pt2_e0_denominator)
end
subroutine selection_collector(b, pt2)
use f77_zmq
use selection_types
use bitmasks
implicit none
type(selection_buffer), intent(inout) :: b
double precision, intent(out) :: pt2(N_states)
double precision :: pt2_mwen(N_states)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_pull_socket
integer(ZMQ_PTR) :: zmq_socket_pull
integer :: msg_size, rc, more
integer :: acc, i, j, robin, N, ntask
double precision, allocatable :: val(:)
integer(bit_kind), allocatable :: det(:,:,:)
integer, allocatable :: task_id(:)
integer :: done
real :: time, time0
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket()
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det))
done = 0
more = 1
pt2(:) = 0d0
call CPU_TIME(time0)
do while (more == 1)
call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask)
pt2 += pt2_mwen
do i=1, N
call add_to_selection_buffer(b, det(1,1,i), val(i))
end do
do i=1, ntask
if(task_id(i) == 0) then
print *, "Error in collector"
endif
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
end do
done += ntask
call CPU_TIME(time)
! print *, "DONE" , done, time - time0
end do
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull)
call sort_selection_buffer(b)
end subroutine

75
plugins/Hartree_Fock/localize_mos.irp.f Normal file
View File

@ -0,0 +1,75 @@
program localize_mos
implicit none
integer :: rank, i,j,k
double precision, allocatable :: W(:,:)
double precision :: f, f_incr
allocate (W(ao_num,ao_num))
W = 0.d0
do k=1,elec_beta_num
do j=1,ao_num
do i=1,ao_num
W(i,j) = W(i,j) + mo_coef(i,k) * mo_coef(j,k)
enddo
enddo
enddo
! call svd_mo(ao_num,elec_beta_num,W, size(W,1), &
! mo_coef(1,1),size(mo_coef,1))
call cholesky_mo(ao_num,elec_beta_num,W, size(W,1), &
mo_coef(1,1),size(mo_coef,1),1.d-6,rank)
print *, rank
if (elec_alpha_num>elec_alpha_num) then
W = 0.d0
do k=elec_beta_num+1,elec_alpha_num
do j=1,ao_num
do i=1,ao_num
W(i,j) = W(i,j) + mo_coef(i,k) * mo_coef(j,k)
enddo
enddo
enddo
! call svd_mo(ao_num,elec_alpha_num-elec_beta_num,W, size(W,1), &
! mo_coef(1,1),size(mo_coef,1))
call cholesky_mo(ao_num,elec_alpha_num-elec_beta_num,W, size(W,1), &
mo_coef(1,elec_beta_num+1),size(mo_coef,1),1.d-6,rank)
print *, rank
endif
W = 0.d0
do k=elec_alpha_num+1,mo_tot_num
do j=1,ao_num
do i=1,ao_num
W(i,j) = W(i,j) + mo_coef(i,k) * mo_coef(j,k)
enddo
enddo
enddo
! call svd_mo(ao_num,mo_tot_num-elec_alpha_num,W, size(W,1), &
! mo_coef(1,1),size(mo_coef,1))
call cholesky_mo(ao_num,mo_tot_num-elec_alpha_num,W, size(W,1), &
mo_coef(1,elec_alpha_num+1),size(mo_coef,1),1.d-6,rank)
print *, rank
mo_label = "Localized"
TOUCH mo_coef
W(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num)
integer :: iorder(mo_tot_num)
double precision :: s(mo_tot_num), swap(ao_num)
do k=1,mo_tot_num
iorder(k) = k
s(k) = Fock_matrix_diag_mo(k)
enddo
call dsort(s(1),iorder(1),elec_beta_num)
call dsort(s(elec_beta_num+1),iorder(elec_beta_num+1),elec_alpha_num-elec_beta_num)
call dsort(s(elec_alpha_num+1),iorder(elec_alpha_num+1),mo_tot_num-elec_alpha_num)
do k=1,mo_tot_num
mo_coef(1:ao_num,k) = W(1:ao_num,iorder(k))
print *, k, s(k)
enddo
call save_mos
end

2
plugins/MRCC_Utils/H_apply.irp.f
View File

@ -31,7 +31,7 @@ s.set_perturbation("epstein_nesbet_2x2")
s.unset_openmp()
print s
s = H_apply_zmq("mrcepa_PT2")
s = H_apply("mrcepa_PT2")
s.energy = "psi_energy"
s.set_perturbation("epstein_nesbet_2x2")
s.unset_openmp()

71
plugins/MRCC_Utils/davidson.irp.f
View File

@ -35,21 +35,20 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,N_st_diag,Ni
PROVIDE mo_bielec_integrals_in_map
allocate(H_jj(sze))
H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(sze,H_jj,N_det_ref,dets_in,Nint,istate,delta_ii,idx_ref) &
!$OMP PRIVATE(i)
!$OMP DO SCHEDULE(guided)
do i=1,sze
!$OMP DO
do i=2,sze
H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint)
enddo
!$OMP END DO
!$OMP DO SCHEDULE(guided)
do i=1,N_det_ref
H_jj(idx_ref(i)) += delta_ii(istate,i)
enddo
!$OMP END DO
!$OMP END PARALLEL
do i=1,N_det_ref
H_jj(idx_ref(i)) += delta_ii(istate,i)
enddo
call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate)
deallocate (H_jj)
end
@ -224,17 +223,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s
W(i,k,iter+1) = 0.d0
enddo
enddo
! do k=1,N_st_diag
! do iter2=1,iter
! do l=1,N_st_diag
! do i=1,sze
! U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
! W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
! enddo
! enddo
! enddo
! enddo
!
!
call dgemm('N','N', sze, N_st_diag, N_st_diag*iter, &
1.d0, U, size(U,1), y, size(y,1)*size(y,2), 0.d0, U(1,1,iter+1), size(U,1))
@ -276,27 +264,11 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s
do k=1,N_st_diag
! do iter2=1,iter
! do l=1,N_st_diag
! c(1) = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze)
! do i=1,sze
! U(i,k,iter+1) = U(i,k,iter+1) - c(1) * U(i,l,iter2)
! enddo
! enddo
! enddo
!
call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), &
U(1,k,iter+1),1,0.d0,c,1)
call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), &
c,1,1.d0,U(1,k,iter+1),1)
!
! do l=1,k-1
! c(1) = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze)
! do i=1,sze
! U(i,k,iter+1) = U(i,k,iter+1) - c(1) * U(i,l,iter+1)
! enddo
! enddo
!
call dgemv('T',sze,k-1,1.d0,U(1,1,iter+1),size(U,1), &
U(1,k,iter+1),1,0.d0,c,1)
call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1), &
@ -429,7 +401,7 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8)
allocate(vt(sze_8,N_st))
Vt = 0.d0
!$OMP DO SCHEDULE(dynamic)
!$OMP DO SCHEDULE(static,1)
do sh=1,shortcut(0,1)
do sh2=sh,shortcut(0,1)
exa = 0
@ -468,9 +440,9 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8)
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP END DO
!$OMP DO SCHEDULE(dynamic)
!$OMP DO SCHEDULE(static,1)
do sh=1,shortcut(0,2)
do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2)
@ -490,7 +462,7 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8)
end do
end do
enddo
!$OMP END DO NOWAIT
!$OMP END DO
!$OMP DO
do ii=1,n_det_ref
@ -505,13 +477,12 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8)
enddo
!$OMP END DO
!$OMP CRITICAL
do istate=1,N_st
do i=n,1,-1
!$OMP ATOMIC
v_0(i,istate) = v_0(i,istate) + vt(i,istate)
enddo
enddo
!$OMP END CRITICAL
deallocate(vt)
!$OMP END PARALLEL
@ -559,25 +530,26 @@ subroutine davidson_diag_mrcc_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_dia
ASSERT (sze > 0)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_bielec_integrals_in_map
allocate(H_jj(sze), S2_jj(sze))
H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
call get_s2(dets_in(1,1,1),dets_in(1,1,1),Nint,S2_jj(1))
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(sze,H_jj,S2_jj, dets_in,Nint,N_det_ref,delta_ii, &
!$OMP idx_ref, istate) &
!$OMP PRIVATE(i)
!$OMP DO SCHEDULE(guided)
do i=1,sze
!$OMP DO
do i=2,sze
H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint)
call get_s2(dets_in(1,1,i),dets_in(1,1,i),Nint,S2_jj(i))
enddo
!$OMP END DO
!$OMP DO SCHEDULE(guided)
!$OMP END PARALLEL
do i=1,N_det_ref
H_jj(idx_ref(i)) += delta_ii(istate,i)
enddo
!$OMP END DO
!$OMP END PARALLEL
call davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate)
deallocate (H_jj,S2_jj)
@ -1094,6 +1066,7 @@ subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_i
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(guided)
do sh=1,shortcut(0,2)
do i=shortcut(sh,2),shortcut(sh+1,2)-1
@ -1142,14 +1115,14 @@ subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_i
! End Specific to dressing
! ------------------------
!$OMP CRITICAL
do istate=1,N_st
do i=n,1,-1
!$OMP ATOMIC
v_0(i,istate) = v_0(i,istate) + vt(istate,i)
!$OMP ATOMIC
s_0(i,istate) = s_0(i,istate) + st(istate,i)
enddo
enddo
!$OMP END CRITICAL
deallocate(vt,st)
!$OMP END PARALLEL

174
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -5,6 +5,7 @@ use bitmasks
END_PROVIDER
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states, N_det_non_ref) ]
&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ]
&BEGIN_PROVIDER [ integer, lambda_mrcc_kept, (0:psi_det_size) ]
@ -62,6 +63,65 @@ END_PROVIDER
END_PROVIDER
! BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states, N_det_non_ref) ]
!&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ]
!&BEGIN_PROVIDER [ integer, lambda_mrcc_kept, (0:psi_det_size) ]
!&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states, N_det_non_ref) ]
! implicit none
! BEGIN_DOC
! ! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
! END_DOC
! integer :: i,k
! double precision :: ihpsi_current(N_states)
! integer :: i_pert_count
! double precision :: hii, E2(N_states), E2var(N_states)
! integer :: N_lambda_mrcc_pt2, N_lambda_mrcc_pt3
!
! i_pert_count = 0
! lambda_mrcc = 0.d0
! N_lambda_mrcc_pt2 = 0
! N_lambda_mrcc_pt3 = 0
! lambda_mrcc_pt2(0) = 0
! lambda_mrcc_kept(0) = 0
!
! E2 = 0.d0
! E2var = 0.d0
! do i=1,N_det_non_ref
! call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
! size(psi_ref_coef,1), N_states,ihpsi_current)
! call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
! do k=1,N_states
! if (ihpsi_current(k) == 0.d0) then
! ihpsi_current(k) = 1.d-32
! endif
! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
! lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
! E2(k) += ihpsi_current(k)*ihpsi_current(k) / (psi_ref_energy_diagonalized(k)-hii)
! E2var(k) += ihpsi_current(k) * psi_non_ref_coef(i,k)
! enddo
! enddo
!
! do i=1,N_det_non_ref
! call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
! size(psi_ref_coef,1), N_states,ihpsi_current)
! call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
! do k=1,N_states
! if (ihpsi_current(k) == 0.d0) then
! ihpsi_current(k) = 1.d-32
! endif
! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
! lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) * E2var(k)/E2(k)
! enddo
! enddo
! lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2
! lambda_mrcc_kept(0) = N_lambda_mrcc_pt3
! print*,'N_det_non_ref = ',N_det_non_ref
! print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
! print*,'lambda max = ',maxval(dabs(lambda_mrcc))
! print*,'Number of ignored determinants = ',i_pert_count
!
!END_PROVIDER
BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
@ -758,10 +818,8 @@ END_PROVIDER
factor = 1.d0
resold = huge(1.d0)
do k=0,10*hh_nex
do k=0,hh_nex/4
res = 0.d0
!$OMP PARALLEL default(shared) private(cx, i, a_col, a_coll) reduction(+:res)
!$OMP DO
do a_coll = 1, n_exc_active
a_col = active_pp_idx(a_coll)
cx = 0.d0
@ -772,21 +830,23 @@ END_PROVIDER
res = res + (X_new(a_col) - X(a_col))*(X_new(a_col) - X(a_col))
X(a_col) = X_new(a_col)
end do
!$OMP END DO
!$OMP END PARALLEL
if (res > resold) then
factor = factor * 0.5d0
endif
if(iand(k, 127) == 0) then
print *, k, res, 1.d0 - res/resold
endif
if ( (res < 1d-10).or.(res/resold > 0.99d0) ) then
exit
endif
resold = res
if(iand(k, 4095) == 0) then
print *, "res ", k, res
end if
if(res < 1d-10) exit
end do
dIj_unique(1:size(X), s) = X(1:size(X))
print *, k, res, 1.d0 - res/resold
enddo
@ -818,21 +878,23 @@ END_PROVIDER
do s=1,N_states
norm = 0.d0
double precision :: f
double precision :: f, g, gmax
gmax = 1.d0*maxval(dabs(psi_non_ref_coef(:,s)))
do i=1,N_det_non_ref
if (rho_mrcc(i,s) == 0.d0) then
rho_mrcc(i,s) = 1.d-32
endif
if (lambda_type == 2) then
f = 1.d0
else
if (rho_mrcc(i,s) == 0.d0) then
cycle
endif
! f is such that f.\tilde{c_i} = c_i
f = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
! Avoid numerical instabilities
f = min(f,2.d0)
f = max(f,-2.d0)
! g = 1.d0+dabs(gmax / psi_non_ref_coef(i,s) )
g = 2.d0+100.d0*exp(-20.d0*dabs(psi_non_ref_coef(i,s)/gmax))
f = min(f, g)
f = max(f,-g)
endif
norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s)
@ -844,24 +906,20 @@ END_PROVIDER
f = 1.d0/norm
! f now contains 1/ <T.Psi_0|T.Psi_0>
norm = 1.d0
do i=1,N_det_ref
norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s)
norm = 0.d0
do i=1,N_det_non_ref
norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s)
enddo
! norm now contains <Psi_SD|Psi_SD>
f = dsqrt(f*norm)
! f normalises T.Psi_0 such that (1+T)|Psi> is normalized
norm = norm*f
print *, 'norm of |T Psi_0> = ', dsqrt(norm)
norm = norm*f
if (dsqrt(norm) > 1.d0) then
stop 'Error : Norm of the SD larger than the norm of the reference.'
endif
do i=1,N_det_ref
norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
enddo
do i=1,N_det_non_ref
rho_mrcc(i,s) = rho_mrcc(i,s) * f
enddo
@ -892,6 +950,53 @@ END_PROVIDER
!double precision function f_fit(x)
! implicit none
! double precision :: x
! f_fit = 0.d0
! return
! if (x < 0.d0) then
! f_fit = 0.d0
! else if (x < 1.d0) then
! f_fit = 1.d0/0.367879441171442 * ( x**2 * exp(-x**2))
! else
! f_fit = 1.d0
! endif
!end
!
!double precision function get_dij_index(II, i, s, Nint)
! integer, intent(in) :: II, i, s, Nint
! double precision, external :: get_dij
! double precision :: HIi, phase, c, a, b, d
!
! call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi)
! call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int)
!
! a = lambda_pert(s,i)
! b = lambda_mrcc(s,i)
! c = f_fit(a/b)
!
! d = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase* rho_mrcc(i,s)
!
! c = f_fit(a*HIi/d)
!
! get_dij_index = HIi * a * c + (1.d0 - c) * d
! get_dij_index = d
! return
!
! if(lambda_type == 0) then
! call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int)
! get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase
! get_dij_index = get_dij_index * rho_mrcc(i,s)
! else if(lambda_type == 1) then
! call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi)
! get_dij_index = HIi * lambda_mrcc(s, i)
! else if(lambda_type == 2) then
! call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int)
! get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase
! get_dij_index = get_dij_index * rho_mrcc(i,s)
! end if
!end function
double precision function get_dij_index(II, i, s, Nint)
integer, intent(in) :: II, i, s, Nint
@ -1031,6 +1136,22 @@ end function
end do
hh_shortcut(hh_shortcut(0)+1) = s+1
if (hh_shortcut(0) > N_hh_exists) then
print *, 'Error in ', irp_here
print *, 'hh_shortcut(0) :', hh_shortcut(0)
print *, 'N_hh_exists : ', N_hh_exists
print *, 'Is your active space defined?'
stop
endif
if (hh_shortcut(hh_shortcut(0)+1)-1 > N_pp_exists) then
print *, 'Error 1 in ', irp_here
print *, 'hh_shortcut(hh_shortcut(0)+1)-1 :', hh_shortcut(hh_shortcut(0)+1)-1
print *, 'N_pp_exists : ', N_pp_exists
print *, 'Is your active space defined?'
stop
endif
do s=2,4,2
do i=1,hh_shortcut(0)
if(hh_exists(s, i) == 0) then
@ -1041,6 +1162,7 @@ end function
end if
end do
do i=1,hh_shortcut(hh_shortcut(0)+1)-1
if(pp_exists(s, i) == 0) then
pp_exists(s-1, i) = 0

34
plugins/Psiref_CAS/psi_ref.irp.f
View File

@ -78,3 +78,37 @@ END_PROVIDER
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, norm_psi_ref, (N_states)]
&BEGIN_PROVIDER [double precision, inv_norm_psi_ref, (N_states)]
implicit none
integer :: i,j
norm_psi_ref = 0.d0
do j = 1, N_states
do i = 1, N_det_ref
norm_psi_ref(j) += psi_ref_coef(i,j) * psi_ref_coef(i,j)
enddo
inv_norm_psi_ref(j) = 1.d0/(dsqrt(norm_psi_Ref(j)))
print *, inv_norm_psi_ref(j)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, psi_ref_coef_interm_norm, (N_det_ref,N_states)]
implicit none
integer :: i,j
do j = 1, N_states
do i = 1, N_det_ref
psi_ref_coef_interm_norm(i,j) = inv_norm_psi_ref(j) * psi_ref_coef(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, psi_non_ref_coef_interm_norm, (N_det_non_ref,N_states)]
implicit none
integer :: i,j
do j = 1, N_states
do i = 1, N_det_non_ref
psi_non_ref_coef_interm_norm(i,j) = psi_non_ref_coef(i,j) * inv_norm_psi_ref(j)
enddo
enddo
END_PROVIDER

3
plugins/Psiref_Utils/psi_ref_utils.irp.f
View File

@ -98,8 +98,7 @@ END_PROVIDER
enddo
N_det_non_ref = i_non_ref
if (N_det_non_ref < 1) then
print *, 'Error : All determinants are in the reference'
stop -1
print *, 'Warning : All determinants are in the reference'
endif
END_PROVIDER

71
plugins/Psiref_threshold/psi_ref.irp.f
View File

@ -1,5 +1,44 @@
use bitmasks
! BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ]
!&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ]
!&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ]
!&BEGIN_PROVIDER [ integer, N_det_ref ]
! implicit none
! BEGIN_DOC
! ! Reference wave function, defined as determinants with amplitudes > 0.05
! ! idx_ref gives the indice of the ref determinant in psi_det.
! END_DOC
! integer :: i, k, l
! logical :: good
! double precision, parameter :: threshold=0.01d0
! double precision :: t(N_states)
! N_det_ref = 0
! do l = 1, N_states
! t(l) = threshold * abs_psi_coef_max(l)
! enddo
! do i=1,N_det
! good = .False.
! do l=1, N_states
! psi_ref_coef(i,l) = 0.d0
! good = good.or.(dabs(psi_coef(i,l)) > t(l))
! enddo
! if (good) then
! N_det_ref = N_det_ref+1
! do k=1,N_int
! psi_ref(k,1,N_det_ref) = psi_det(k,1,i)
! psi_ref(k,2,N_det_ref) = psi_det(k,2,i)
! enddo
! idx_ref(N_det_ref) = i
! do k=1,N_states
! psi_ref_coef(N_det_ref,k) = psi_coef(i,k)
! enddo
! endif
! enddo
! call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference')
!
!END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ]
&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ]
@ -10,30 +49,16 @@ use bitmasks
! idx_ref gives the indice of the ref determinant in psi_det.
END_DOC
integer :: i, k, l
logical :: good
double precision, parameter :: threshold=0.05d0
double precision :: t(N_states)
N_det_ref = 0
do l = 1, N_states
t(l) = threshold * abs_psi_coef_max(l)
enddo
do i=1,N_det
good = .False.
do l=1, N_states
psi_ref_coef(i,l) = 0.d0
good = good.or.(dabs(psi_coef(i,l)) > t(l))
double precision, parameter :: threshold=0.01d0
call find_reference(threshold, N_det_ref, idx_ref)
do l=1,N_states
do i=1,N_det_ref
psi_ref_coef(i,l) = psi_coef(idx_ref(i), l)
enddo
if (good) then
N_det_ref = N_det_ref+1
do k=1,N_int
psi_ref(k,1,N_det_ref) = psi_det(k,1,i)
psi_ref(k,2,N_det_ref) = psi_det(k,2,i)
enddo
idx_ref(N_det_ref) = i
do k=1,N_states
psi_ref_coef(N_det_ref,k) = psi_coef(i,k)
enddo
endif
enddo
do i=1,N_det_ref
psi_ref(:,:,i) = psi_det(:,:,idx_ref(i))
enddo
call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference')

1
plugins/analyze_wf/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Determinants

12
plugins/analyze_wf/README.rst Normal file
View File

@ -0,0 +1,12 @@
==========
analyze_wf
==========
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

70
plugins/analyze_wf/analyze_wf.irp.f Normal file
View File

@ -0,0 +1,70 @@
program analyze_wf
implicit none
BEGIN_DOC
! Wave function analyzis
END_DOC
read_wf = .True.
SOFT_TOUCH read_wf
call run()
end
subroutine run
implicit none
integer :: istate, i
integer :: class(0:mo_tot_num,5)
double precision :: occupation(mo_tot_num)
write(*,'(A)') 'MO Occupation'
write(*,'(A)') '============='
write(*,'(A)') ''
do istate=1,N_states
call get_occupation_from_dets(occupation,1)
write(*,'(A)') ''
write(*,'(A,I3)'), 'State ', istate
write(*,'(A)') '---------------'
write(*,'(A)') ''
write (*,'(A)') '======== ================'
class = 0
do i=1,mo_tot_num
write (*,'(I8,X,F16.10)') i, occupation(i)
if (occupation(i) > 1.999d0) then
class(0,1) += 1
class( class(0,1), 1) = i
else if (occupation(i) > 1.95d0) then
class(0,2) += 1
class( class(0,2), 2) = i
else if (occupation(i) < 0.001d0) then
class(0,5) += 1
class( class(0,5), 5) = i
else if (occupation(i) < 0.01d0) then
class(0,4) += 1
class( class(0,4), 4) = i
else
class(0,3) += 1
class( class(0,3), 3) = i
endif
enddo
write (*,'(A)') '======== ================'
write (*,'(A)') ''
write (*,'(A)') 'Suggested classes'
write (*,'(A)') '-----------------'
write (*,'(A)') ''
write (*,'(A)') 'Core :'
write (*,*) (class(i,1), ',', i=1,class(0,1))
write (*,*) ''
write (*,'(A)') 'Inactive :'
write (*,*) (class(i,2), ',', i=1,class(0,2))
write (*,'(A)') ''
write (*,'(A)') 'Active :'
write (*,*) (class(i,3), ',', i=1,class(0,3))
write (*,'(A)') ''
write (*,'(A)') 'Virtual :'
write (*,*) (class(i,4), ',', i=1,class(0,4))
write (*,'(A)') ''
write (*,'(A)') 'Deleted :'
write (*,*) (class(i,5), ',', i=1,class(0,5))
write (*,'(A)') ''
enddo
end

23
plugins/analyze_wf/occupation.irp.f Normal file
View File

@ -0,0 +1,23 @@
subroutine get_occupation_from_dets(occupation, istate)
implicit none
double precision, intent(out) :: occupation(mo_tot_num)
integer, intent(in) :: istate
BEGIN_DOC
! Returns the average occupation of the MOs
END_DOC
integer :: i,j, ispin
integer :: list(N_int*bit_kind_size,2)
integer :: n_elements(2)
double precision :: c
occupation = 0.d0
do i=1,N_det
c = psi_coef(i,istate)*psi_coef(i,istate)
call bitstring_to_list_ab(psi_det(1,1,i), list, n_elements, N_int)
do ispin=1,2
do j=1,n_elements(ispin)
occupation( list(j,ispin) ) += c
enddo
enddo
enddo
end

33
plugins/mrcc_selected/EZFIO.cfg Normal file
View File

@ -0,0 +1,33 @@
[lambda_type]
type: Positive_int
doc: lambda type
interface: ezfio,provider,ocaml
default: 0
[energy]
type: double precision
doc: Calculated energy
interface: ezfio
[energy_pt2]
type: double precision
doc: Calculated energy with PT2 contribution
interface: ezfio
[energy]
type: double precision
doc: Calculated energy
interface: ezfio
[thresh_dressed_ci]
type: Threshold
doc: Threshold on the convergence of the dressed CI energy
interface: ezfio,provider,ocaml
default: 1.e-5
[n_it_max_dressed_ci]
type: Strictly_positive_int
doc: Maximum number of dressed CI iterations
interface: ezfio,provider,ocaml
default: 10

1
plugins/mrcc_selected/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Perturbation Selectors_full Generators_full Psiref_threshold MRCC_Utils ZMQ

12
plugins/mrcc_selected/README.rst Normal file
View File

@ -0,0 +1,12 @@
=======
mrcepa0
=======
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

80
plugins/mrcc_selected/dressing.irp.f
View File

@ -534,63 +534,9 @@ END_PROVIDER
END_PROVIDER
! BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ]
! use bitmasks
! implicit none
! integer :: i,j,k
! double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall
! integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref)
!
! ! provide lambda_mrcc
! npres = 0
! delta_cas = 0d0
! call wall_time(wall)
! print *, "dcas ", wall
! do i_state = 1, N_states
! !!$OMP PARALLEL DO default(none) schedule(dynamic) private(pre,npre,ipre,j,k,Hjk,Hki,degree) shared(npres,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref)
! do k=1,N_det_non_ref
! if(lambda_mrcc(i_state, k) == 0d0) cycle
! npre = 0
! do i=1,N_det_ref
! call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki)
! if(Hki /= 0d0) then
! !!$OMP ATOMIC
! npres(i) += 1
! npre += 1
! ipre(npre) = i
! pre(npre) = Hki
! end if
! end do
!
!
! do i=1,npre
! do j=1,i
! !!$OMP ATOMIC
! delta_cas(ipre(i),ipre(j),i_state) += pre(i) * pre(j) * lambda_mrcc(i_state, k)
! end do
! end do
! end do
! !!$OMP END PARALLEL DO
! npre=0
! do i=1,N_det_ref
! npre += npres(i)
! end do
! !stop
! do i=1,N_det_ref
! do j=1,i
! delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
! end do
! end do
! end do
!
! call wall_time(wall)
! print *, "dcas", wall
! ! stop
! END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ]
&BEGIN_PROVIDER [ double precision, delta_cas_s2, (N_det_ref, N_det_ref, N_states) ]
BEGIN_PROVIDER [ double precision, delta_ref, (N_det_ref, N_det_ref, N_states) ]
&BEGIN_PROVIDER [ double precision, delta_ref_s2, (N_det_ref, N_det_ref, N_states) ]
use bitmasks
implicit none
integer :: i,j,k
@ -600,22 +546,22 @@ END_PROVIDER
provide lambda_mrcc dIj
do i_state = 1, N_states
!$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Sjk,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,delta_cas_s2,N_det_ref,dij)
!$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Sjk,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_ref,delta_ref_s2,N_det_ref,dij)
do i=1,N_det_ref
do j=1,i
call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int)
delta_cas(i,j,i_state) = 0d0
delta_cas_s2(i,j,i_state) = 0d0
delta_ref(i,j,i_state) = 0d0
delta_ref_s2(i,j,i_state) = 0d0
do k=1,N_det_non_ref
call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk)
call get_s2(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Sjk)
delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k)
delta_cas_s2(i,j,i_state) += Sjk * dij(i, k, i_state) ! * Ski * lambda_mrcc(i_state, k)
delta_ref(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k)
delta_ref_s2(i,j,i_state) += Sjk * dij(i, k, i_state) ! * Ski * lambda_mrcc(i_state, k)
end do
delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
delta_cas_s2(j,i,i_state) = delta_cas_s2(i,j,i_state)
delta_ref(j,i,i_state) = delta_ref(i,j,i_state)
delta_ref_s2(j,i,i_state) = delta_ref_s2(i,j,i_state)
end do
end do
!$OMP END PARALLEL DO
@ -739,7 +685,7 @@ end function
!$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii, delta_mrcepa0_ij_s2, delta_mrcepa0_ii_s2) &
!$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2,contrib_s2,contrib2_s2) &
!$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) &
!$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas, delta_cas_s2) &
!$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_ref, delta_ref_s2) &
!$OMP shared(notf,i_state, sortRef, sortRefIdx, dij)
do blok=1,cepa0_shortcut(0)
do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
@ -781,8 +727,8 @@ end function
notf = notf+1
! call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
contrib_s2 = delta_cas_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
contrib = delta_ref(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
contrib_s2 = delta_ref_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
if(dabs(psi_ref_coef(J,i_state)).ge.1.d-3) then
contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
@ -828,7 +774,7 @@ END_PROVIDER
integer :: II, blok
provide delta_cas lambda_mrcc
provide delta_ref lambda_mrcc
allocate(idx_sorted_bit(N_det))
idx_sorted_bit(:) = -1
do i=1,N_det_non_ref

24
plugins/mrcc_selected/dressing_slave.irp.f
View File

@ -294,12 +294,12 @@ subroutine push_mrsc2_results(zmq_socket_push, I_i, J, delta, delta_s2, task_id)
endif
! ! Activate is zmq_socket_push is a REQ
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
end
@ -368,12 +368,12 @@ subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, delta_s2,
! ! Activate is zmq_socket_pull is a REP
! integer :: idummy
! rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)'
! stop 'error'
! endif
integer :: idummy
rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)'
stop 'error'
endif
end

2
plugins/mrcc_selected/ezfio_interface.irp.f
View File

@ -1,6 +1,6 @@
! DO NOT MODIFY BY HAND
! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
! from file /home/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg
! from file /ccc/work/cont003/gen1738/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg
BEGIN_PROVIDER [ double precision, thresh_dressed_ci ]

1
plugins/mrcc_selected/mrcc_selected.irp.f
View File

@ -8,7 +8,6 @@ program mrsc2sub
read_wf = .True.
SOFT_TOUCH read_wf
call print_cas_coefs
call set_generators_bitmasks_as_holes_and_particles
call run(N_states,energy)
if(do_pt2_end)then

15
plugins/mrcc_selected/mrcepa0_general.irp.f
View File

@ -60,16 +60,17 @@ subroutine run(N_st,energy)
end
subroutine print_cas_coefs
subroutine print_ref_coefs
implicit none
integer :: i,j
print *, 'CAS'
print *, '==='
do i=1,N_det_cas
print *, (psi_cas_coef(i,j), j=1,N_states)
call debug_det(psi_cas(1,1,i),N_int)
print *, 'Reference'
print *, '========='
do i=1,N_det_ref
print *, (psi_ref_coef(i,j), j=1,N_states)
call debug_det(psi_ref(1,1,i),N_int)
enddo
print *, ''
call write_double(6,ci_energy(1),"Initial CI energy")
end
@ -202,7 +203,7 @@ subroutine run_pt2(N_st,energy)
print*,'Last iteration only to compute the PT2'
N_det_generators = N_det_cas
N_det_generators = N_det_ref
N_det_selectors = N_det_non_ref
do i=1,N_det_generators

6
plugins/mrcepa0/EZFIO.cfg
View File

@ -14,6 +14,12 @@ type: double precision
doc: Calculated energy with PT2 contribution
interface: ezfio
[perturbative_triples]
type: logical
doc: Compute perturbative contribution of the Triples
interface: ezfio,provider,ocaml
default: true
[energy]
type: double precision
doc: Calculated energy

103
plugins/mrcepa0/dressing.irp.f
View File

@ -13,6 +13,7 @@ use bitmasks
integer(bit_kind),allocatable :: buf(:,:,:)
logical :: ok
logical, external :: detEq
integer, external :: omp_get_thread_num
delta_ij_mrcc = 0d0
delta_ii_mrcc = 0d0
@ -23,7 +24,7 @@ use bitmasks
!$OMP PARALLEL DO default(none) schedule(dynamic) &
!$OMP shared(psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) &
!$OMP shared(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc, delta_ii_s2_mrcc, delta_ij_s2_mrcc) &
!$OMP private(h, n, mask, omask, buf, ok, iproc)
!$OMP private(h, n, mask, omask, buf, ok, iproc)
do gen= 1, N_det_generators
allocate(buf(N_int, 2, N_det_non_ref))
iproc = omp_get_thread_num() + 1
@ -74,9 +75,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
logical :: good, fullMatch
integer(bit_kind),allocatable :: tq(:,:,:)
integer :: N_tq, c_ref ,degree
integer :: N_tq, c_ref ,degree1, degree2, degree
double precision :: hIk, hla, hIl, sla, dIk(N_states), dka(N_states), dIa(N_states)
double precision :: hIk, hla, hIl, sla, dIk(N_states), dka(N_states), dIa(N_states), hka
double precision, allocatable :: dIa_hla(:,:), dIa_sla(:,:)
double precision :: haj, phase, phase2
double precision :: f(N_states), ci_inv(N_states)
@ -99,6 +100,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
!double precision, external :: get_dij, get_dij_index
leng = max(N_det_generators, N_det_non_ref)
allocate(miniList(Nint, 2, leng), tq(Nint,2,n_selected), idx_minilist(leng), hij_cache(N_det_non_ref), sij_cache(N_det_non_ref))
allocate(idx_alpha(0:psi_det_size), degree_alpha(psi_det_size))
@ -189,17 +191,25 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
end do
end if
if (perturbative_triples) then
double precision :: Delta_E_inv(N_states)
double precision, external :: diag_H_mat_elem
do i_state=1,N_states
Delta_E_inv(i_state) = 1.d0 / (psi_ref_energy_diagonalized(i_state) - diag_H_mat_elem(tq(1,1,i_alpha),N_int) )
enddo
endif
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd))
call get_s2(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,sij_cache(k_sd))
enddo
! |I>
do i_I=1,N_det_ref
! Find triples and quadruple grand parents
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint)
if (degree > 4) then
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree1,Nint)
if (degree1 > 4) then
cycle
endif
@ -209,77 +219,57 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
! <I| <> |alpha>
do k_sd=1,idx_alpha(0)
! Loop if lambda == 0
logical :: loop
! loop = .True.
! do i_state=1,N_states
! if (lambda_mrcc(i_state,idx_alpha(k_sd)) /= 0.d0) then
! loop = .False.
! exit
! endif
! enddo
! if (loop) then
! cycle
! endif
call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
if (degree > 2) then
cycle
endif
! <I| /k\ |alpha>
! <I|H|k>
!hIk = hij_mrcc(idx_alpha(k_sd),i_I)
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
do i_state=1,N_states
dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state)
!dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
!dIk(i_state) = psi_non_ref_coef(idx_alpha(k_sd), i_state) / psi_ref_coef(i_I, i_state)
enddo
! |l> = Exc(k -> alpha) |I>
call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree2,phase,Nint)
call decode_exc(exc,degree2,h1,p1,h2,p2,s1,s2)
do k=1,N_int
tmp_det(k,1) = psi_ref(k,1,i_I)
tmp_det(k,2) = psi_ref(k,2,i_I)
enddo
logical :: ok
call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
if(.not. ok) cycle
if (perturbative_triples) then
ok = ok .and. ( (degree2 /= 1).and.(degree /=1) )
endif
do i_state=1,N_states
dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state)
enddo
! <I| \l/ |alpha>
do i_state=1,N_states
dka(i_state) = 0.d0
enddo
do l_sd=k_sd+1,idx_alpha(0)
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
if (degree == 0) then
! loop = .True.
! do i_state=1,N_states
! if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then
! loop = .False.
! exit
! endif
! enddo
loop = .false.
if (.not.loop) then
if (ok) then
do l_sd=k_sd+1,idx_alpha(0)
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
if (degree == 0) then
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
hIl = hij_mrcc(idx_alpha(l_sd),i_I)
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
do i_state=1,N_states
dka(i_state) = dij(i_I, idx_alpha(l_sd), i_state) * phase * phase2
!dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
!dka(i_state) = psi_non_ref_coef(idx_alpha(l_sd), i_state) / psi_ref_coef(i_I, i_state) * phase * phase2
enddo
exit
endif
enddo
else if (perturbative_triples) then
hka = hij_cache(idx_alpha(k_sd))
do i_state=1,N_states
dka(i_state) = hka * Delta_E_inv(i_state)
enddo
endif
exit
endif
enddo
do i_state=1,N_states
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
enddo
@ -292,32 +282,35 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
k_sd = idx_alpha(l_sd)
hla = hij_cache(k_sd)
sla = sij_cache(k_sd)
! call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla)
do i_state=1,N_states
dIa_hla(i_state,k_sd) = dIa(i_state) * hla
dIa_sla(i_state,k_sd) = dIa(i_state) * sla
enddo
enddo
call omp_set_lock( psi_ref_lock(i_I) )
do i_state=1,N_states
if(dabs(psi_ref_coef(i_I,i_state)).ge.1.d-3)then
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
!$OMP ATOMIC
delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
!$OMP ATOMIC
delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
!$OMP ATOMIC
delta_ij_s2_(i_state,k_sd,i_I) = delta_ij_s2_(i_state,k_sd,i_I) + dIa_sla(i_state,k_sd)
!$OMP ATOMIC
delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
enddo
else
delta_ii_(i_state,i_I) = 0.d0
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
!$OMP ATOMIC
delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0*dIa_hla(i_state,k_sd)
!$OMP ATOMIC
delta_ij_s2_(i_state,k_sd,i_I) = delta_ij_s2_(i_state,k_sd,i_I) + 0.5d0*dIa_sla(i_state,k_sd)
enddo
endif
enddo
call omp_unset_lock( psi_ref_lock(i_I) )
enddo
enddo
deallocate (dIa_hla,dIa_sla,hij_cache,sij_cache)

28
plugins/mrcepa0/dressing_slave.irp.f
View File

@ -315,13 +315,13 @@ subroutine push_mrsc2_results(zmq_socket_push, I_i, J, delta, delta_s2, task_id)
stop 'error'
endif
! ! Activate is zmq_socket_push is a REQ
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
! Activate is zmq_socket_push is a REQ
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
end
@ -389,13 +389,13 @@ subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, delta_s2,
endif
! ! Activate is zmq_socket_pull is a REP
! integer :: idummy
! rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)'
! stop 'error'
! endif
! Activate is zmq_socket_pull is a REP
integer :: idummy
rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)'
stop 'error'
endif
end

2
plugins/mrcepa0/mrcc.irp.f
View File

@ -5,7 +5,7 @@ program mrsc2sub
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
mrmode = 3
read_wf = .True.
SOFT_TOUCH read_wf
call set_generators_bitmasks_as_holes_and_particles

1
plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Psiref_CAS Determinants Davidson

12
plugins/mrsc2_no_amp/README.rst Normal file
View File

@ -0,0 +1,12 @@
============
mrsc2_no_amp
============
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

129
plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f Normal file
View File

@ -0,0 +1,129 @@
BEGIN_PROVIDER [double precision, CI_eigenvectors_sc2_no_amp, (N_det,N_states_diag)]
&BEGIN_PROVIDER [double precision, CI_eigenvectors_s2_sc2_no_amp, (N_states_diag)]
&BEGIN_PROVIDER [double precision, CI_electronic_energy_sc2_no_amp, (N_states_diag)]
implicit none
integer :: i,j,k,l
integer, allocatable :: idx(:)
integer, allocatable :: holes_part(:,:)
double precision, allocatable :: e_corr(:,:)
double precision, allocatable :: accu(:)
double precision, allocatable :: ihpsi_current(:)
double precision, allocatable :: H_jj(:),H_jj_total(:),S2_jj(:)
integer :: number_of_particles, number_of_holes, n_h,n_p
allocate(e_corr(N_det_non_ref,N_states),ihpsi_current(N_states),accu(N_states),H_jj(N_det_non_ref),idx(0:N_det_non_ref))
allocate(H_jj_total(N_det),S2_jj(N_det))
allocate(holes_part(N_det,2))
accu = 0.d0
do i = 1, N_det_non_ref
holes_part(i,1) = number_of_holes(psi_non_ref(1,1,i))
holes_part(i,2) = number_of_particles(psi_non_ref(1,1,i))
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
size(psi_ref_coef_interm_norm,1), N_states,ihpsi_current)
do j = 1, N_states
e_corr(i,j) = psi_non_ref_coef(i,j) * ihpsi_current(j) * inv_norm_psi_ref(j)
accu(j) += e_corr(i,j)
enddo
enddo
print *, 'accu = ',accu
double precision :: hjj,diag_h_mat_elem
do i = 1, N_det_non_ref
H_jj(i) = 0.d0
n_h = holes_part(i,1)
n_p = holes_part(i,2)
integer :: degree
! do j = 1, N_det_non_ref
! call get_excitation_degree(psi_non_ref(1,1,i),psi_non_ref(1,1,j),degree,N_int)
! if(degree .gt. 2)then
! if(n_h + holes_part(j,1) .gt. 2 .or. n_p + holes_part(j,2) .gt. 2 ) then
! H_jj(i) += e_corr(j,1)
! endif
! endif
! enddo
call filter_not_connected(psi_non_ref,psi_non_ref(1,1,i),N_int,N_det_non_ref,idx)
do j = 1, idx(0)
if(n_h + holes_part(idx(j),1) .gt. 2 .or. n_p + holes_part(idx(j),2) .gt. 2 ) then
H_jj(i) += e_corr(idx(j),1)
endif
enddo
enddo
do i=1,N_Det
H_jj_total(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
call get_s2(psi_det(1,1,i),psi_det(1,1,i),N_int,S2_jj(i))
enddo
do i = 1, N_det_non_ref
H_jj_total(idx_non_ref(i)) += H_jj(i)
enddo
print *, 'coef'
call davidson_diag_hjj_sjj(psi_det,CI_eigenvectors_sc2_no_amp,H_jj_total,S2_jj,CI_electronic_energy_sc2_no_amp,size(CI_eigenvectors_sc2_no_amp,1),N_Det,N_states,N_states_diag,N_int,6)
do i = 1, N_det
hjj = diag_h_mat_elem(psi_det(1,1,i),N_int)
! if(hjj<-210.d0)then
! call debug_det(psi_det(1,1,i),N_int)
! print *, CI_eigenvectors_sc2_no_amp((i),1),hjj, H_jj_total(i)
! endif
enddo
print *, 'ref',N_det_ref
do i =1, N_det_ref
call debug_det(psi_det(1,1,idx_ref(i)),N_int)
print *, CI_eigenvectors_sc2_no_amp(idx_ref(i),1), H_jj_total(idx_ref(i))
enddo
print *, 'non ref',N_det_non_ref
do i=1, N_det_non_ref
hjj = diag_h_mat_elem(psi_non_ref(1,1,i),N_int)
! print *, CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1),H_jj_total(idx_non_ref(i)), H_jj(i)
! if(dabs(CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1)).gt.1.d-1)then
! if(hjj<-210.d0)then
! call debug_det(psi_det(1,1,idx_non_ref(i)),N_int)
! write(*,'(10(F16.10,X))') CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1),hjj, H_jj(i),H_jj_total(idx_non_ref(i))
! endif
enddo
! do i = 1, N_det
! print *, CI_eigenvectors_sc2_no_amp(i,1)
! enddo
do i=1,N_states_diag
CI_eigenvectors_s2_sc2_no_amp(i) = S2_jj(i)
enddo
deallocate(e_corr,ihpsi_current,accu,H_jj,idx,H_jj_total,s2_jj,holes_part)
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_energy_sc2_no_amp, (N_states_diag) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
END_DOC
integer :: j
character*(8) :: st
call write_time(output_determinants)
do j=1,min(N_det,N_states_diag)
CI_energy_sc2_no_amp(j) = CI_electronic_energy_sc2_no_amp(j) + nuclear_repulsion
enddo
do j=1,min(N_det,N_states)
write(st,'(I4)') j
call write_double(output_determinants,CI_energy_sc2_no_amp(j),'Energy of state '//trim(st))
call write_double(output_determinants,CI_eigenvectors_s2_sc2_no_amp(j),'S^2 of state '//trim(st))
enddo
END_PROVIDER
subroutine diagonalize_CI_sc2_no_amp
implicit none
integer :: i,j
do j=1,N_states
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors_sc2_no_amp(i,j)
enddo
enddo
SOFT_TOUCH ci_eigenvectors_s2_sc2_no_amp ci_eigenvectors_sc2_no_amp ci_electronic_energy_sc2_no_amp ci_energy_sc2_no_amp psi_coef
end

14
plugins/mrsc2_no_amp/sc2_no_amp.irp.f Normal file
View File

@ -0,0 +1,14 @@
program pouet
provide ao_bielec_integrals_in_map
call bla
end
subroutine bla
implicit none
integer :: i
do i = 1, 10
call diagonalize_CI_sc2_no_amp
TOUCH psi_coef
enddo
print *, "E+PT2 = ", ci_energy_sc2_no_amp(:)
end

272
promela/integrals.pml Normal file
View File

@ -0,0 +1,272 @@
#define NPROC 1
#define BUFSIZE 2
#define NTASKS 3
mtype = { NONE, OK, WRONG_STATE, TERMINATE, GETPSI, PUTPSI, NEWJOB, ENDJOB, SETRUNNING,
SETWAITING, SETSTOPPED, CONNECT, DISCONNECT, ADDTASK, DELTASK, TASKDONE, GETTASK,
PSI, TASK, PUTPSI_REPLY, WAITING, RUNNING, STOPPED
}
typedef rep_message {
mtype m = NONE;
byte value = 0;
}
typedef req_message {
mtype m = NONE;
byte state = 0;
byte value = 0;
chan reply = [BUFSIZE] of { rep_message };
}
#define send_req( MESSAGE, VALUE ) msg.m=MESSAGE ; msg.value=VALUE ; msg.state=state; rep_socket ! msg; msg.reply ? reply
chan rep_socket = [NPROC] of { req_message };
chan pull_socket = [NPROC] of { byte };
chan pair_socket = [NPROC] of { req_message };
chan task_queue = [NTASKS+2] of { byte };
chan pub_socket = [NTASKS+2] of { mtype };
bit socket_up = 0;
mtype global_state; /* Sent by pub */
active proctype qp_run() {
bit psi = 0;
bit address_tcp = 0;
bit address_inproc = 0;
bit running = 0;
byte status = 0;
byte state = 0;
byte ntasks = 0;
req_message msg;
rep_message reply;
byte nclients = 0;
byte task;
socket_up = 1;
running = 1;
do
// :: ( (running == 0) && (nclients == 0) && (ntasks == 0) ) -> break
:: ( running == 0 ) -> break
:: else ->
rep_socket ? msg;
printf("req: "); printm(msg.m); printf("\t%d\n",msg.value);
if
:: ( msg.m == TERMINATE ) ->
assert (state != 0);
assert (msg.state == state);
running = 0;
reply.m = OK;
:: ( msg.m == PUTPSI ) ->
assert (state != 0);
assert (msg.state == state);
assert (psi == 0);
psi = 1;
reply.m = PUTPSI_REPLY;
:: ( msg.m == GETPSI ) ->
assert (state != 0);
assert (msg.state == state);
assert (psi == 1);
reply.m = PSI;
:: ( msg.m == NEWJOB ) ->
assert (state == 0);
state = msg.value;
pair_socket ! WAITING;
reply.m = OK;
reply.value = state;
:: ( msg.m == ENDJOB ) ->
assert (state != 0);
assert (msg.state == state);
state = 0;
pair_socket ! WAITING;
reply.m = OK;
:: ( msg.m == ADDTASK ) ->
assert (state != 0);
assert (msg.state == state);
task_queue ! msg.value;
ntasks++;
reply.m = OK;
:: ( msg.m == GETTASK ) ->
assert (nclients > 0);
assert (state != 0);
assert (msg.state == state);
if
:: ( task_queue ?[task] ) ->
pair_socket ! WAITING;
reply.m = TASK;
task_queue ? reply.value
:: else ->
pair_socket ! RUNNING;
reply.m = NONE;
reply.value = 255;
fi;
:: ( msg.m == TASKDONE) ->
assert (state != 0);
assert (msg.state == state);
assert (nclients > 0);
assert (ntasks > 0);
reply.m = OK;
:: ( msg.m == DELTASK ) ->
assert (state != 0);
assert (msg.state == state);
ntasks--;
if
:: (ntasks > 0) -> reply.value = 1;
:: else -> reply.value = 0;
fi;
reply.m = OK;
:: ( msg.m == CONNECT ) ->
assert ( state != 0 )
nclients++;
reply.m = OK;
reply.value = state;
:: ( msg.m == DISCONNECT ) ->
assert ( msg.state == state )
nclients--;
reply.m = OK;
:: ( msg.m == STOPPED ) ->
pair_socket ! STOPPED;
reply.m = OK;
:: ( msg.m == WAITING ) ->
pair_socket ! WAITING;
reply.m = OK;
:: ( msg.m == RUNNING ) ->
assert ( state != 0 );
pair_socket ! RUNNING;
reply.m = OK;
fi
msg.reply ! reply
od
pair_socket ! STOPPED;
socket_up = 0;
}
active proctype master() {
req_message msg;
rep_message reply;
byte state = 0;
byte count;
run pub_thread();
/* New parallel job */
state=1;
send_req( NEWJOB, state );
assert (reply.m == OK);
/* Add tasks */
count = 0;
do
:: (count == NTASKS) -> break;
:: else ->
count++;
send_req( ADDTASK, count );
assert (reply.m == OK);
od
/* Run collector */
run collector(state);
/* Run slaves */
count = 0;
do
:: (count == NPROC) -> break;
:: else -> count++; run slave();
od
}
proctype slave() {
req_message msg;
rep_message reply;
byte task;
byte state;
msg.m=CONNECT;
msg.state = 0;
if
:: (!socket_up) -> goto exit;
:: else -> skip;
fi
rep_socket ! msg;
if
:: (!socket_up) -> goto exit;
:: else -> skip;
fi
msg.reply ? reply;
state = reply.value;
task = 1;
do
:: (task == 255) -> break;
:: else ->
send_req( GETTASK, 0);
if
:: (reply.m == NONE) ->
task = 255;
:: (reply.m == TASK) ->
/* Compute task */
task = reply.value;
send_req( TASKDONE, task);
assert (reply.m == OK);
pull_socket ! task;
fi
od
send_req( DISCONNECT, 0);
assert (reply.m == OK);
exit: skip;
}
proctype collector(byte state) {
byte task;
req_message msg;
rep_message reply;
bit loop = 1;
do
:: (loop == 0) -> break
:: else ->
pull_socket ? task;
/* Handle result */
send_req(DELTASK, task);
assert (reply.m == OK);
loop = reply.value;
od
send_req( TERMINATE, 0);
assert (reply.m == OK);
}
proctype pub_thread() {
mtype state = WAITING;
do
:: (state == STOPPED) -> break;
:: (pair_socket ? [state]) ->
pair_socket ? state;
global_state = state;
od
}

5
src/AO_Basis/aos.irp.f
View File

@ -182,7 +182,7 @@ integer function ao_power_index(nx,ny,nz)
end
BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)]
BEGIN_PROVIDER [ character*(128), l_to_charater, (0:7)]
BEGIN_DOC
! character corresponding to the "L" value of an AO orbital
END_DOC
@ -192,6 +192,9 @@ BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)]
l_to_charater(2)='D'
l_to_charater(3)='F'
l_to_charater(4)='G'
l_to_charater(5)='H'
l_to_charater(6)='I'
l_to_charater(7)='J'
END_PROVIDER

9
src/Bitmask/bitmasks.irp.f
View File

@ -6,6 +6,7 @@ BEGIN_PROVIDER [ integer, N_int ]
! Number of 64-bit integers needed to represent determinants as binary strings
END_DOC
N_int = (mo_tot_num-1)/bit_kind_size + 1
call write_int(6,N_int, 'N_int')
END_PROVIDER
@ -386,6 +387,8 @@ END_PROVIDER
n_virt_orb += popcnt(virt_bitmask(i,1))
enddo
endif
call write_int(6,n_inact_orb, 'Number of inactive MOs')
call write_int(6,n_virt_orb, 'Number of virtual MOs')
END_PROVIDER
@ -554,7 +557,7 @@ END_PROVIDER
&BEGIN_PROVIDER [ integer, n_core_orb]
implicit none
BEGIN_DOC
! Core orbitals bitmask
! Core + deleted orbitals bitmask
END_DOC
integer :: i,j
n_core_orb = 0
@ -563,7 +566,7 @@ END_PROVIDER
core_bitmask(i,2) = xor(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,2),virt_bitmask(i,1)))
n_core_orb += popcnt(core_bitmask(i,1))
enddo
print*,'n_core_orb = ',n_core_orb
call write_int(6,n_core_orb,'Number of core MOs')
END_PROVIDER
@ -598,7 +601,7 @@ BEGIN_PROVIDER [ integer, n_act_orb]
do i = 1, N_int
n_act_orb += popcnt(cas_bitmask(i,1,1))
enddo
print*,'n_act_orb = ',n_act_orb
call write_int(6,n_act_orb, 'Number of active MOs')
END_PROVIDER
BEGIN_PROVIDER [integer, list_act, (n_act_orb)]

6
src/Davidson/EZFIO.cfg
View File

@ -28,3 +28,9 @@ doc: If true, disk space is used to store the vectors
default: False
interface: ezfio,provider,ocaml
[distributed_davidson]
type: logical
doc: If true, use the distributed algorithm
default: False
interface: ezfio,provider,ocaml

237
src/Davidson/davidson_parallel.irp.f
View File

@ -20,15 +20,16 @@ subroutine davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep)
double precision :: s2, hij
logical, allocatable :: wrotten(:)
PROVIDE dav_det ref_bitmask_energy
allocate(wrotten(bs))
wrotten = .false.
PROVIDE dav_det
ii=0
sh = blockb
do sh2=1,shortcut_(0,1)
exa = 0
do ni=1,N_int
exa = popcnt(xor(version_(1,sh,1), version_(1,sh2,1)))
do ni=2,N_int
exa = exa + popcnt(xor(version_(ni,sh,1), version_(ni,sh2,1)))
end do
if(exa > 2) cycle
@ -43,14 +44,18 @@ subroutine davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep)
do j=shortcut_(sh2,1), shortcut_(sh2+1,1)-1
if(i == j) cycle
org_j = sort_idx_(j,1)
ext = exa
do ni=1,N_int
ext = exa + popcnt(xor(sorted_i(1), sorted_(1,j,1)))
if(ext > 4) cycle
do ni=2,N_int
ext = ext + popcnt(xor(sorted_i(ni), sorted_(ni,j,1)))
if(ext > 4) exit
end do
if(ext <= 4) then
call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2)
org_j = sort_idx_(j,1)
call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2)
! call i_h_j (sorted_(1,j,1),sorted_(1,i,1),n_int,hij)
! call get_s2(sorted_(1,j,1),sorted_(1,i,1),n_int,s2)
if(.not. wrotten(ii)) then
wrotten(ii) = .true.
idx(ii) = org_i
@ -58,8 +63,8 @@ subroutine davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep)
st (:,ii) = 0d0
end if
do istate=1,N_states_diag
vt (istate,ii) += hij*dav_ut(istate,org_j)
st (istate,ii) += s2*dav_ut(istate,org_j)
vt (istate,ii) = vt (istate,ii) +hij*dav_ut(istate,org_j)
st (istate,ii) = st (istate,ii) +s2*dav_ut(istate,org_j)
enddo
endif
enddo
@ -67,32 +72,40 @@ subroutine davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep)
enddo
if (blockb <= shortcut_(0,2)) then
if ( blockb <= shortcut_(0,2) ) then
sh=blockb
do sh2=sh, shortcut_(0,2), shortcut_(0,1)
do i=blockb2+shortcut_(sh2,2),shortcut_(sh2+1,2)-1, istep
ii += 1
if (ii>bs) then
print *, irp_here
stop 'ii>bs'
endif
org_i = sort_idx_(i,2)
do j=shortcut_(sh2,2),shortcut_(sh2+1,2)-1
if(i == j) cycle
org_j = sort_idx_(j,2)
ext = 0
do ni=1,N_int
ext = popcnt(xor(sorted_(1,i,2), sorted_(1,j,2)))
if (ext > 4) cycle
do ni=2,N_int
ext = ext + popcnt(xor(sorted_(ni,i,2), sorted_(ni,j,2)))
if (ext > 4) exit
end do
if(ext == 4) then
call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2)
if(.not. wrotten(ii)) then
wrotten(ii) = .true.
idx(ii) = org_i
vt (:,ii) = 0d0
st (:,ii) = 0d0
end if
do istate=1,N_states_diag
vt (istate,ii) += hij*dav_ut(istate,org_j)
st (istate,ii) += s2*dav_ut(istate,org_j)
enddo
call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2)
! call i_h_j (sorted_(1,j,2),sorted_(1,i,2),n_int,hij)
! call get_s2(sorted_(1,j,2),sorted_(1,i,2),n_int,s2)
if(.not. wrotten(ii)) then
wrotten(ii) = .true.
idx(ii) = org_i
vt (:,ii) = 0d0
st (:,ii) = 0d0
end if
do istate=1,N_states_diag
vt (istate,ii) = vt (istate,ii) +hij*dav_ut(istate,org_j)
st (istate,ii) = st (istate,ii) +s2*dav_ut(istate,org_j)
enddo
end if
end do
end do
@ -128,10 +141,8 @@ subroutine davidson_collect(N, idx, vt, st , v0t, s0t)
integer :: i, j, k
!DIR$ IVDEP
do i=1,N
k = idx(i)
!DIR$ IVDEP
do j=1,N_states_diag
v0t(j,k) = v0t(j,k) + vt(j,i)
s0t(j,k) = s0t(j,k) + st(j,i)
@ -140,53 +151,42 @@ subroutine davidson_collect(N, idx, vt, st , v0t, s0t)
end subroutine
subroutine davidson_init(zmq_to_qp_run_socket,n,n_st_8,ut)
subroutine davidson_init(zmq_to_qp_run_socket,dets_in,u,n0,n,n_st,update_dets)
use f77_zmq
implicit none
integer(ZMQ_PTR), intent(out) :: zmq_to_qp_run_socket
integer, intent(in) :: n, n_st_8
double precision, intent(in) :: ut(n_st_8,n)
integer, intent(in) :: n0,n, n_st, update_dets
double precision, intent(in) :: u(n0,n_st)
integer(bit_kind), intent(in) :: dets_in(N_int,2,n)
integer :: i,k
dav_size = n
touch dav_size
if (update_dets == 1) then
dav_size = n
touch dav_size
do i=1,dav_size
do k=1,N_int
dav_det(k,1,i) = dets_in(k,1,i)
dav_det(k,2,i) = dets_in(k,2,i)
enddo
enddo
touch dav_det
endif
do i=1,n
do k=1,N_int
dav_det(k,1,i) = psi_det(k,1,i)
dav_det(k,2,i) = psi_det(k,2,i)
enddo
enddo
do i=1,n
do k=1,N_states_diag
dav_ut(k,i) = ut(k,i)
do k=1,n_st
dav_ut(k,i) = u(i,k)
enddo
enddo
touch dav_det dav_ut
soft_touch dav_ut
call new_parallel_job(zmq_to_qp_run_socket,"davidson")
end subroutine
subroutine davidson_add_task(zmq_to_qp_run_socket, blockb, blockb2, istep)
use f77_zmq
implicit none
integer(ZMQ_PTR) ,intent(in) :: zmq_to_qp_run_socket
integer ,intent(in) :: blockb, blockb2, istep
character*(512) :: task
write(task,*) blockb, blockb2, istep
call add_task_to_taskserver(zmq_to_qp_run_socket, task)
end subroutine
subroutine davidson_slave_inproc(i)
implicit none
integer, intent(in) :: i
@ -281,6 +281,7 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id)
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
call davidson_push_results(zmq_socket_push, blockb, blockb2, N, idx, vt, st, task_id)
end do
deallocate(idx, vt, st)
end subroutine
@ -320,6 +321,15 @@ subroutine davidson_push_results(zmq_socket_push, blockb, blocke, N, idx, vt, st
rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0)
if(rc /= 4) stop "davidson_push_results failed to push task_id"
! Activate is zmq_socket_push is a REQ
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
end subroutine
@ -358,6 +368,14 @@ subroutine davidson_pull_results(zmq_socket_pull, blockb, blocke, N, idx, vt, st
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
if(rc /= 4) stop "davidson_pull_results failed to pull task_id"
! Activate if zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
if (rc /= 4) then
print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
stop 'error'
endif
end subroutine
@ -390,8 +408,8 @@ subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0, LD
allocate(v0t(N_states_diag, dav_size))
allocate(s0t(N_states_diag, dav_size))
v0t = 00.d0
s0t = 00.d0
v0t = 0.d0
s0t = 0.d0
more = 1
@ -404,9 +422,7 @@ subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0, LD
deallocate(idx,vt,st)
integer :: i,j
!DIR$ IVDEP
do j=1,N_states_diag
!DIR$ IVDEP
do i=1,dav_size
v0(i,j) = v0t(j,i)
s0(i,j) = s0t(j,i)
@ -434,37 +450,22 @@ subroutine davidson_run(zmq_to_qp_run_socket , v0, s0, LDA)
double precision , intent(inout) :: v0(LDA, N_states_diag)
double precision , intent(inout) :: s0(LDA, N_states_diag)
call zmq_set_running(zmq_to_qp_run_socket)
zmq_collector = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket()
i = omp_get_thread_num()
PROVIDE nproc
!$OMP PARALLEL NUM_THREADS(nproc+2) PRIVATE(i)
i = omp_get_thread_num()
if (i == 0 ) then
call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0, LDA)
call end_zmq_to_qp_run_socket(zmq_collector)
call end_zmq_pull_socket(zmq_socket_pull)
call davidson_miniserver_end()
else if (i == 1 ) then
call davidson_miniserver_run ()
else
call davidson_slave_inproc(i)
endif
!$OMP END PARALLEL
call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0, LDA)
call end_zmq_to_qp_run_socket(zmq_collector)
call end_zmq_pull_socket(zmq_socket_pull)
call davidson_miniserver_end()
call end_parallel_job(zmq_to_qp_run_socket, 'davidson')
end subroutine
subroutine davidson_miniserver_run()
subroutine davidson_miniserver_run(update_dets)
use f77_zmq
implicit none
integer update_dets
integer(ZMQ_PTR) responder
character*(64) address
character(len=:), allocatable :: buffer
@ -473,18 +474,23 @@ subroutine davidson_miniserver_run()
allocate (character(len=20) :: buffer)
address = 'tcp://*:11223'
PROVIDE dav_det dav_ut dav_size
responder = f77_zmq_socket(zmq_context, ZMQ_REP)
rc = f77_zmq_bind(responder,address)
do
rc = f77_zmq_recv(responder, buffer, 5, 0)
if (buffer(1:rc) /= 'end') then
rc = f77_zmq_send (responder, dav_size, 4, ZMQ_SNDMORE)
rc = f77_zmq_send (responder, dav_det, 16*N_int*dav_size, ZMQ_SNDMORE)
rc = f77_zmq_send (responder, dav_ut, 8*dav_size*N_states_diag, 0)
else
if (buffer(1:rc) == 'end') then
rc = f77_zmq_send (responder, "end", 3, 0)
exit
else if (buffer(1:rc) == 'det') then
rc = f77_zmq_send (responder, dav_size, 4, ZMQ_SNDMORE)
rc = f77_zmq_send (responder, dav_det, 16*N_int*dav_size, 0)
else if (buffer(1:rc) == 'ut') then
rc = f77_zmq_send (responder, update_dets, 4, ZMQ_SNDMORE)
rc = f77_zmq_send (responder, dav_size, 4, ZMQ_SNDMORE)
rc = f77_zmq_send (responder, dav_ut, 8*dav_size*N_states_diag, 0)
endif
enddo
@ -511,34 +517,63 @@ subroutine davidson_miniserver_end()
end subroutine
subroutine davidson_miniserver_get()
subroutine davidson_miniserver_get(force_update)
implicit none
use f77_zmq
logical, intent(in) :: force_update
integer(ZMQ_PTR) requester
character*(64) address
character*(20) buffer
integer rc
integer rc, update_dets
address = trim(qp_run_address)//':11223'
requester = f77_zmq_socket(zmq_context, ZMQ_REQ)
rc = f77_zmq_connect(requester,address)
rc = f77_zmq_send(requester, "Hello", 5, 0)
rc = f77_zmq_recv(requester, dav_size, 4, 0)
TOUCH dav_size
rc = f77_zmq_recv(requester, dav_det, 16*N_int*dav_size, 0)
rc = f77_zmq_recv(requester, dav_ut, 8*dav_size*N_states_diag, 0)
TOUCH dav_det dav_ut
rc = f77_zmq_send(requester, 'ut', 2, 0)
rc = f77_zmq_recv(requester, update_dets, 4, 0)
if (rc /= 4) then
print *, irp_here, ': f77_zmq_recv(requester, update_dets, 4, 0)'
print *, irp_here, ': rc = ', rc
endif
rc = f77_zmq_recv(requester, dav_size, 4, 0)
if (rc /= 4) then
print *, irp_here, ': f77_zmq_recv(requester, dav_size, 4, 0)'
print *, irp_here, ': rc = ', rc
endif
if (update_dets == 1 .or. force_update) then
TOUCH dav_size
endif
rc = f77_zmq_recv(requester, dav_ut, 8*dav_size*N_states_diag, 0)
if (rc /= 8*dav_size*N_states_diag) then
print *, irp_here, ': f77_zmq_recv(requester, dav_ut, 8*dav_size*N_states_diag, 0)'
print *, irp_here, ': rc = ', rc
endif
SOFT_TOUCH dav_ut
if (update_dets == 1 .or. force_update) then
rc = f77_zmq_send(requester, 'det', 3, 0)
rc = f77_zmq_recv(requester, dav_size, 4, 0)
if (rc /= 4) then
print *, irp_here, ': f77_zmq_recv(requester, dav_size, 4, 0)'
print *, irp_here, ': rc = ', rc
endif
rc = f77_zmq_recv(requester, dav_det, 16*N_int*dav_size, 0)
if (rc /= 16*N_int*dav_size) then
print *, irp_here, ': f77_zmq_recv(requester, dav_det, 16*N_int*dav_size, 0)'
print *, irp_here, ': rc = ', rc
endif
SOFT_TOUCH dav_det
endif
end subroutine
BEGIN_PROVIDER [ integer(bit_kind), dav_det, (N_int, 2, dav_size) ]
&BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ]
use bitmasks
implicit none
BEGIN_DOC
@ -546,7 +581,19 @@ end subroutine
!
! Touched in davidson_miniserver_get
END_DOC
integer :: i,k
dav_det = 0_bit_kind
END_PROVIDER
BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ]
use bitmasks
implicit none
BEGIN_DOC
! Temporary arrays for parallel davidson
!
! Touched in davidson_miniserver_get
END_DOC
dav_ut = -huge(1.d0)
END_PROVIDER

6
src/Davidson/davidson_slave.irp.f
View File

@ -7,6 +7,7 @@ program davidson_slave
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states_diag)
character*(64) :: state
logical :: force_update
call provide_everything
call switch_qp_run_to_master
@ -16,11 +17,12 @@ program davidson_slave
state = 'Waiting'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
force_update = .True.
do
call wait_for_state(zmq_state,state)
if(trim(state) /= "davidson") exit
call davidson_miniserver_get()
call davidson_miniserver_get(force_update)
force_update = .False.
integer :: rc, i

27
src/Davidson/diagonalization_hs2.irp.f
View File

@ -110,7 +110,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
character*(16384) :: write_buffer
double precision :: to_print(3,N_st)
double precision :: cpu, wall
integer :: shift, shift2, itermax
integer :: shift, shift2, itermax, update_dets
double precision :: r1, r2
logical :: state_ok(N_st_diag*davidson_sze_max)
include 'constants.include.F'
@ -122,6 +122,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
stop -1
endif
integer, external :: align_double
sze_8 = align_double(sze)
itermax = max(3,min(davidson_sze_max, sze/N_st_diag))
PROVIDE nuclear_repulsion expected_s2
call write_time(iunit)
@ -134,6 +138,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
call write_int(iunit,N_st,'Number of states')
call write_int(iunit,N_st_diag,'Number of states in diagonalization')
call write_int(iunit,sze,'Number of determinants')
r1 = 8.d0*(3.d0*dble(sze_8*N_st_diag*itermax+5.d0*(N_st_diag*itermax)**2 &
+ 4.d0*(N_st_diag*itermax))/(1024.d0**3))
call write_double(iunit, r1, 'Memory(Gb)')
write(iunit,'(A)') ''
write_buffer = '===== '
do i=1,N_st
@ -151,14 +158,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
enddo
write(iunit,'(A)') trim(write_buffer)
integer, external :: align_double
sze_8 = align_double(sze)
itermax = max(3,min(davidson_sze_max, sze/N_st_diag))
allocate( &
! Large
W(sze_8,N_st_diag*itermax), &
U(sze_8,N_st_diag*itermax), &
S(sze_8,N_st_diag*itermax), &
! Small
h(N_st_diag*itermax,N_st_diag*itermax), &
y(N_st_diag*itermax,N_st_diag*itermax), &
s_(N_st_diag*itermax,N_st_diag*itermax), &
@ -204,6 +211,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
enddo
update_dets = 1
do while (.not.converged)
do k=1,N_st_diag
@ -223,8 +232,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! -----------------------------------------
! call H_S2_u_0_nstates_zmq(W(1,shift+1),S(1,shift+1),U(1,shift+1),H_jj,S2_jj,sze,dets_in,Nint,N_st_diag,sze_8)
call H_S2_u_0_nstates(W(1,shift+1),S(1,shift+1),U(1,shift+1),H_jj,S2_jj,sze,dets_in,Nint,N_st_diag,sze_8)
if (distributed_davidson) then
call H_S2_u_0_nstates_zmq(W(1,shift+1),S(1,shift+1),U(1,shift+1),H_jj,S2_jj,sze,dets_in,Nint,N_st_diag,sze_8,update_dets)
else
call H_S2_u_0_nstates(W(1,shift+1),S(1,shift+1),U(1,shift+1),H_jj,S2_jj,sze,dets_in,Nint,N_st_diag,sze_8)
endif
update_dets = 0
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>

41
src/Davidson/find_reference.irp.f Normal file
View File

@ -0,0 +1,41 @@
subroutine find_reference(thresh,n_ref,result)
implicit none
double precision, intent(in) :: thresh
integer, intent(out) :: result(N_det),n_ref
integer :: i,j,istate
double precision :: i_H_psi_array(1), E0, hii, norm
double precision :: de
integer(bit_kind), allocatable :: psi_ref_(:,:,:)
double precision, allocatable :: psi_ref_coef_(:,:)
allocate(psi_ref_coef_(N_det,1), psi_ref_(N_int,2,N_det))
n_ref = 1
result(1) = 1
istate = 1
psi_ref_coef_(1,1) = psi_coef(1,istate)
psi_ref_(:,:,1) = psi_det(:,:,1)
norm = psi_ref_coef_(1,1) * psi_ref_coef_(1,1)
call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1))
print *, ''
print *, 'Reference determinants'
print *, '======================'
print *, ''
print *, n_ref, ': E0 = ', E0 + nuclear_repulsion
call debug_det(psi_ref_(1,1,n_ref),N_int)
do i=2,N_det
call i_h_psi(psi_det(1,1,i),psi_ref_(1,1,1),psi_ref_coef_(1,istate),N_int, &
n_ref,size(psi_ref_coef_,1),1,i_H_psi_array)
call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hii)
de = i_H_psi_array(istate)**2 / (E0 - hii)
if (dabs(de) > thresh) then
n_ref += 1
result(n_ref) = i
psi_ref_(:,:,n_ref) = psi_det(:,:,i)
psi_ref_coef_(n_ref,1) = psi_coef(i,istate)
call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1))
print *, n_ref, ': E0 = ', E0 + nuclear_repulsion
call debug_det(psi_ref_(1,1,n_ref),N_int)
endif
enddo
end

5
src/Davidson/parameters.irp.f
View File

@ -18,6 +18,11 @@ subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged
double precision :: E(N_st), time
double precision, allocatable, save :: energy_old(:)
if (iterations < 2) then
converged = .False.
return
endif
if (.not.allocated(energy_old)) then
allocate(energy_old(N_st))
energy_old = 0.d0

626
src/Davidson/u0Hu0.irp.f
View File

@ -32,272 +32,18 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
! Computes v_0 = H|u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
!
END_DOC
integer, intent(in) :: N_st,n,Nint, sze_8
double precision, intent(out) :: v_0(sze_8,N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
double precision, intent(in) :: H_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: hij
double precision, allocatable :: vt(:,:)
double precision, allocatable :: ut(:,:)
integer :: i,j,k,l, jj,ii
integer :: i0, j0
integer, allocatable :: shortcut(:,:), sort_idx(:,:)
integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate
integer :: N_st_8
integer, external :: align_double
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut
N_st_8 = align_double(N_st)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
allocate(ut(N_st_8,n))
v_0 = 0.d0
do i=1,n
do istate=1,N_st
ut(istate,i) = u_0(i,istate)
enddo
enddo
call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)&
!$OMP SHARED(n,H_jj,keys_tmp,ut,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,N_st_8)
allocate(vt(N_st_8,n))
Vt = 0.d0
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,1)
do sh2=1,shortcut(0,1)
exa = popcnt(xor(version(1,sh,1), version(1,sh2,1)))
if(exa > 2) then
cycle
end if
do ni=2,Nint
exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1)))
end do
if(exa > 2) then
cycle
end if
do i=shortcut(sh,1),shortcut(sh+1,1)-1
org_i = sort_idx(i,1)
do ni=1,Nint
sorted_i(ni) = sorted(ni,i,1)
enddo
jloop: do j=shortcut(sh2,1),shortcut(sh2+1,1)-1
org_j = sort_idx(j,1)
ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1)))
if(ext > 4) then
cycle jloop
endif
do ni=2,Nint
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
if(ext > 4) then
cycle jloop
endif
end do
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
do istate=1,N_st
vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j)
enddo
enddo jloop
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,2)
do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2)
do j=shortcut(sh,2),shortcut(sh+1,2)-1
org_j = sort_idx(j,2)
ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2)))
do ni=2,Nint
ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2)))
end do
if(ext /= 4) then
cycle
endif
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
do istate=1,N_st
vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j)
enddo
end do
end do
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
do istate=1,N_st
do i=n,1,-1
v_0(i,istate) = v_0(i,istate) + vt(istate,i)
enddo
enddo
!$OMP END CRITICAL
deallocate(vt)
!$OMP END PARALLEL
do istate=1,N_st
do i=1,n
v_0(i,istate) += H_jj(i) * u_0(i,istate)
enddo
enddo
deallocate (shortcut, sort_idx, sorted, version, ut)
end
BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]
implicit none
BEGIN_DOC
! Energy of the current wave function
END_DOC
call u_0_H_u_0(psi_energy,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size)
END_PROVIDER
subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
use f77_zmq
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
!
! S2_jj : array of <j|S^2|j>
END_DOC
integer, intent(in) :: N_st,n,Nint, sze_8
double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
double precision, intent(in) :: H_jj(n), S2_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: hij,s2
double precision, allocatable :: ut(:,:)
integer :: i,j,k,l, jj,ii
integer :: i0, j0
integer, allocatable :: shortcut(:,:), sort_idx(:)
integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate
integer :: N_st_8
integer, external :: align_double
integer :: blockb, blockb2, istep
double precision :: ave_workload, workload, target_workload_inv
integer(ZMQ_PTR) :: handler
if(N_st /= N_states_diag .or. sze_8 < N_det) stop "assert fail in H_S2_u_0_nstates"
N_st_8 = N_st ! align_double(N_st)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n), sorted(Nint,n), version(Nint,n))
allocate(ut(N_st_8,n))
v_0 = 0.d0
s_0 = 0.d0
do i=1,n
do istate=1,N_st
ut(istate,i) = u_0(i,istate)
enddo
enddo
call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut(0,1), version, n, Nint)
call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut(0,2), version, n, Nint)
blockb = shortcut(0,1)
call davidson_init(handler,n,N_st_8,ut)
ave_workload = 0.d0
do sh=1,shortcut(0,1)
ave_workload += shortcut(0,1)
ave_workload += (shortcut(sh+1,1) - shortcut(sh,1))**2
do i=sh, shortcut(0,2), shortcut(0,1)
do j=i, min(i, shortcut(0,2))
ave_workload += (shortcut(j+1,2) - shortcut(j, 2))**2
end do
end do
enddo
ave_workload = ave_workload/dble(shortcut(0,1))
target_workload_inv = 0.001d0/ave_workload
do sh=1,shortcut(0,1),1
workload = shortcut(0,1)+dble(shortcut(sh+1,1) - shortcut(sh,1))**2
do i=sh, shortcut(0,2), shortcut(0,1)
do j=i, min(i, shortcut(0,2))
workload += (shortcut(j+1,2) - shortcut(j, 2))**2
end do
end do
istep = 1+ int(workload*target_workload_inv)
do blockb2=0, istep-1
call davidson_add_task(handler, sh, blockb2, istep)
enddo
enddo
call davidson_run(handler, v_0, s_0, size(v_0,1))
do istate=1,N_st
do i=1,n
v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate)
s_0(i,istate) = s_0(i,istate) + s2_jj(i)* u_0(i,istate)
enddo
enddo
deallocate(shortcut, sort_idx, sorted, version)
deallocate(ut)
end
subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
!
! S2_jj : array of <j|S^2|j>
END_DOC
integer, intent(in) :: N_st,n,Nint, sze_8
double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
double precision, intent(in) :: H_jj(n), S2_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: hij,s2
double precision, allocatable :: vt(:,:), ut(:,:), st(:,:)
integer :: i,j,k,l, jj,ii
@ -311,8 +57,6 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
integer :: N_st_8
integer, external :: align_double
integer :: blockb, blockb2, istep
double precision :: ave_workload, workload, target_workload_inv
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut, st
@ -324,17 +68,16 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
allocate(ut(N_st_8,n))
allocate( ut(N_st_8,n))
v_0 = 0.d0
s_0 = 0.d0
call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)&
!$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,s_0,sorted,shortcut,sort_idx,version,N_st,N_st_8)
!$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,sorted,shortcut,sort_idx,version,N_st,N_st_8)
allocate(vt(N_st_8,n),st(N_st_8,n))
Vt = 0.d0
St = 0.d0
@ -347,7 +90,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic)
!$OMP DO SCHEDULE(guided)
do sh=1,shortcut(0,2)
do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2)
@ -380,7 +123,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic)
!$OMP DO SCHEDULE(guided)
do sh=1,shortcut(0,1)
do sh2=1,shortcut(0,1)
if (sh==sh2) cycle
@ -492,14 +235,367 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
enddo
!$OMP END DO
!$OMP CRITICAL (u0Hu0)
do istate=1,N_st
do i=1,n
!$OMP ATOMIC
v_0(i,istate) = v_0(i,istate) + vt(istate,i)
enddo
enddo
deallocate(vt,st)
!$OMP END PARALLEL
do istate=1,N_st
do i=1,n
v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate)
enddo
enddo
deallocate (shortcut, sort_idx, sorted, version, ut)
end
BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]
implicit none
BEGIN_DOC
! Energy of the current wave function
END_DOC
call u_0_H_u_0(psi_energy,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size)
END_PROVIDER
subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8,update_dets)
use omp_lib
use bitmasks
use f77_zmq
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
!
! S2_jj : array of <j|S^2|j>
END_DOC
integer, intent(in) :: N_st,n,Nint, sze_8, update_dets
double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
double precision, intent(in) :: H_jj(n), S2_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: hij,s2
integer :: i,j,k,l, jj,ii
integer :: i0, j0, ithread
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate
integer :: N_st_8
integer, external :: align_double
integer :: blockb2, istep
double precision :: ave_workload, workload, target_workload_inv
integer(ZMQ_PTR) :: handler
if(N_st /= N_states_diag .or. sze_8 < N_det) stop "assert fail in H_S2_u_0_nstates"
N_st_8 = N_st ! align_double(N_st)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy
v_0 = 0.d0
s_0 = 0.d0
call davidson_init(handler,keys_tmp,u_0,size(u_0,1),n,N_st,update_dets)
ave_workload = 0.d0
do sh=1,shortcut_(0,1)
ave_workload += shortcut_(0,1)
ave_workload += (shortcut_(sh+1,1) - shortcut_(sh,1))**2
do i=sh, shortcut_(0,2), shortcut_(0,1)
do j=i, min(i, shortcut_(0,2))
ave_workload += (shortcut_(j+1,2) - shortcut_(j, 2))**2
end do
end do
enddo
ave_workload = ave_workload/dble(shortcut_(0,1))
target_workload_inv = 0.01d0/ave_workload
PROVIDE nproc
character(len=:), allocatable :: task
task = repeat(' ', iposmax)
character(32) :: tmp_task
integer :: ipos, iposmax
iposmax = shortcut_(0,1)+32
ipos = 1
do sh=1,shortcut_(0,1),1
workload = shortcut_(0,1)+dble(shortcut_(sh+1,1) - shortcut_(sh,1))**2
do i=sh, shortcut_(0,2), shortcut_(0,1)
do j=i, min(i, shortcut_(0,2))
workload += (shortcut_(j+1,2) - shortcut_(j, 2))**2
end do
end do
! istep = 1+ int(workload*target_workload_inv)
istep = 1
do blockb2=0, istep-1
write(tmp_task,'(3(I9,X),''|'',X)') sh, blockb2, istep
task = task//tmp_task
ipos += 32
if (ipos+32 > iposmax) then
call add_task_to_taskserver(handler, trim(task))
ipos=1
task = ''
endif
enddo
enddo
if (ipos>1) then
call add_task_to_taskserver(handler, trim(task))
endif
!$OMP PARALLEL NUM_THREADS(nproc+2) PRIVATE(ithread)
ithread = omp_get_thread_num()
if (ithread == 0 ) then
call zmq_set_running(handler)
call davidson_run(handler, v_0, s_0, size(v_0,1))
else if (ithread == 1 ) then
call davidson_miniserver_run (update_dets)
else
call davidson_slave_inproc(ithread)
endif
!$OMP END PARALLEL
call end_parallel_job(handler, 'davidson')
do istate=1,N_st
do i=1,n
v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate)
s_0(i,istate) = s_0(i,istate) + s2_jj(i)* u_0(i,istate)
enddo
enddo
end
subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
!
! S2_jj : array of <j|S^2|j>
END_DOC
integer, intent(in) :: N_st,n,Nint, sze_8
double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
double precision, intent(in) :: H_jj(n), S2_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: hij,s2
double precision, allocatable :: vt(:,:), ut(:,:), st(:,:)
integer :: i,j,k,l, jj,ii
integer :: i0, j0
integer, allocatable :: shortcut(:,:), sort_idx(:,:)
integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate
integer :: N_st_8
integer, external :: align_double
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut, st
N_st_8 = align_double(N_st)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
allocate( ut(N_st_8,n))
v_0 = 0.d0
s_0 = 0.d0
call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)&
!$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,s_0,sorted,shortcut,sort_idx,version,N_st,N_st_8)
allocate(vt(N_st_8,n),st(N_st_8,n))
Vt = 0.d0
St = 0.d0
!$OMP DO
do i=1,n
do istate=1,N_st
ut(istate,i) = u_0(sort_idx(i,2),istate)
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(guided)
do sh=1,shortcut(0,2)
do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2)
do j=shortcut(sh,2),shortcut(sh+1,2)-1
org_j = sort_idx(j,2)
ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2)))
if (ext > 4) cycle
do ni=2,Nint
ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2)))
if (ext > 4) exit
end do
if(ext == 4) then
call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij)
call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2)
do istate=1,n_st
vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j)
st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j)
enddo
end if
end do
end do
enddo
!$OMP END DO
!$OMP DO
do i=1,n
do istate=1,N_st
ut(istate,i) = u_0(sort_idx(i,1),istate)
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(guided)
do sh=1,shortcut(0,1)
do sh2=1,shortcut(0,1)
if (sh==sh2) cycle
exa = 0
do ni=1,Nint
exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1)))
end do
if(exa > 2) then
cycle
end if
do i=shortcut(sh,1),shortcut(sh+1,1)-1
org_i = sort_idx(i,1)
do ni=1,Nint
sorted_i(ni) = sorted(ni,i,1)
enddo
do j=shortcut(sh2,1),shortcut(sh2+1,1)-1
ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1)))
if (ext > 4) cycle
do ni=2,Nint
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
if (ext > 4) exit
end do
if(ext <= 4) then
org_j = sort_idx(j,1)
call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij)
if (hij /= 0.d0) then
do istate=1,n_st
vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j)
enddo
endif
if (ext /= 2) then
call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2)
if (s2 /= 0.d0) then
do istate=1,n_st
st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j)
enddo
endif
endif
endif
enddo
enddo
enddo
exa = 0
do i=shortcut(sh,1),shortcut(sh+1,1)-1
org_i = sort_idx(i,1)
do ni=1,Nint
sorted_i(ni) = sorted(ni,i,1)
enddo
do j=shortcut(sh,1),i-1
ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1)))
if (ext > 4) cycle
do ni=2,Nint
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
if (ext > 4) exit
end do
if(ext <= 4) then
org_j = sort_idx(j,1)
call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij)
if (hij /= 0.d0) then
do istate=1,n_st
vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j)
enddo
endif
if (ext /= 2) then
call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2)
if (s2 /= 0.d0) then
do istate=1,n_st
st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j)
enddo
endif
endif
endif
enddo
do j=i+1,shortcut(sh+1,1)-1
ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1)))
if (ext > 4) cycle
do ni=2,Nint
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
if (ext > 4) exit
end do
if(ext <= 4) then
org_j = sort_idx(j,1)
call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij)
if (hij /= 0.d0) then
do istate=1,n_st
vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j)
enddo
endif
if (ext /= 2) then
call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2)
if (s2 /= 0.d0) then
do istate=1,n_st
st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j)
enddo
endif
endif
endif
enddo
enddo
enddo
!$OMP END DO
do istate=1,N_st
do i=1,n
!$OMP ATOMIC
v_0(i,istate) = v_0(i,istate) + vt(istate,i)
!$OMP ATOMIC
s_0(i,istate) = s_0(i,istate) + st(istate,i)
enddo
enddo
!$OMP END CRITICAL (u0Hu0)
deallocate(vt,st)
!$OMP END PARALLEL

2
src/Determinants/EZFIO.cfg
View File

@ -38,7 +38,7 @@ default: False
type: logical
doc: Force the wave function to be an eigenfunction of S^2
interface: ezfio,provider,ocaml
default: False
default: True
[threshold_generators]
type: Threshold

9
src/Determinants/Fock_diag.irp.f
View File

@ -19,6 +19,15 @@ subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint)
fock_diag_tmp = 0.d0
E0 = 0.d0
if (Ne(1) /= elec_alpha_num) then
print *, 'Error in build_fock_tmp (alpha)', Ne(1), Ne(2)
stop -1
endif
if (Ne(2) /= elec_beta_num) then
print *, 'Error in build_fock_tmp (beta)', Ne(1), Ne(2)
stop -1
endif
! Occupied MOs
do ii=1,elec_alpha_num
i = occ(ii,1)

22
src/Determinants/H_apply.irp.f
View File

@ -438,12 +438,12 @@ subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,i_generator,N_st,t
endif
! Activate if zmq_socket_push is a REQ
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
end
subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,i_generator,N_st,n,task_id)
@ -509,11 +509,11 @@ subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,i_generator,N_st,n
endif
! Activate if zmq_socket_pull is a REP
! rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_pull, 0, 4, 0)'
! stop 'error'
! endif
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_pull, 0, 4, 0)'
stop 'error'
endif
end

2
src/Determinants/H_apply_zmq.template.f
View File

@ -38,7 +38,7 @@ subroutine $subroutine($params_main)
do i_generator=1,N_det_generators
$skip
write(task,*) i_generator
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task))
enddo
allocate ( pt2_generators(N_states,N_det_generators), &

88
src/Determinants/density_matrix.irp.f
View File

@ -78,25 +78,33 @@ END_PROVIDER
double precision :: ck, cl, ckl
double precision :: phase
integer :: h1,h2,p1,p2,s1,s2, degree
integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int,2)
integer :: exc(0:2,2,2),n_occ(2)
double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:)
integer :: krow, kcol, lrow, lcol
one_body_dm_mo_alpha = 0.d0
one_body_dm_mo_beta = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
!$OMP tmp_a, tmp_b, n_occ)&
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
!$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
!$OMP mo_tot_num)
!$OMP mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, &
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
tmp_a = 0.d0
tmp_b = 0.d0
!$OMP DO SCHEDULE(dynamic)
!$OMP DO SCHEDULE(guided)
do k=1,N_det
call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int)
krow = psi_bilinear_matrix_rows(k)
kcol = psi_bilinear_matrix_columns(k)
tmp_det(:,1) = psi_det(:,1, krow)
tmp_det(:,2) = psi_det(:,2, kcol)
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
do m=1,N_states
ck = psi_coef(k,m)*psi_coef(k,m)
ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
do l=1,elec_alpha_num
j = occ(l,1)
tmp_a(j,j,m) += ck
@ -106,24 +114,61 @@ END_PROVIDER
tmp_b(j,j,m) += ck
enddo
enddo
do l=1,k-1
call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
if (degree /= 1) then
cycle
l = k+1
lrow = psi_bilinear_matrix_rows(l)
lcol = psi_bilinear_matrix_columns(l)
do while ( lcol == kcol )
tmp_det2(:,1) = psi_det(:,1, lrow)
tmp_det2(:,2) = psi_det(:,2, lcol)
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
if (degree == 1) then
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do m=1,N_states
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
if (s1==1) then
tmp_a(h1,p1,m) += ckl
tmp_a(p1,h1,m) += ckl
else
tmp_b(h1,p1,m) += ckl
tmp_b(p1,h1,m) += ckl
endif
enddo
endif
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do m=1,N_states
ckl = psi_coef(k,m) * psi_coef(l,m) * phase
if (s1==1) then
tmp_a(h1,p1,m) += ckl
tmp_a(p1,h1,m) += ckl
else
tmp_b(h1,p1,m) += ckl
tmp_b(p1,h1,m) += ckl
endif
enddo
l = l+1
if (l>N_det) exit
lrow = psi_bilinear_matrix_rows(l)
lcol = psi_bilinear_matrix_columns(l)
enddo
l = k+1
lrow = psi_bilinear_matrix_transp_rows(l)
lcol = psi_bilinear_matrix_transp_columns(l)
do while ( lrow == krow )
tmp_det2(:,1) = psi_det(:,1, lrow)
tmp_det2(:,2) = psi_det(:,2, lcol)
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
if (degree == 1) then
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do m=1,N_states
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_transp_values(l,m) * phase
if (s1==1) then
tmp_a(h1,p1,m) += ckl
tmp_a(p1,h1,m) += ckl
else
tmp_b(h1,p1,m) += ckl
tmp_b(p1,h1,m) += ckl
endif
enddo
endif
l = l+1
if (l>N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l)
lcol = psi_bilinear_matrix_transp_columns(l)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
@ -244,7 +289,6 @@ subroutine set_natural_mos
double precision, allocatable :: tmp(:,:)
label = "Natural"
! call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,label,-1)
call mo_as_svd_vectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,label)
end

98
src/Determinants/filter_connected.irp.f
View File

@ -1,4 +1,102 @@
subroutine filter_not_connected(key1,key2,Nint,sze,idx)
use bitmasks
implicit none
BEGIN_DOC
! Returns the array idx which contains the index of the
!
! determinants in the array key1 that DO NOT interact
!
! via the H operator with key2.
!
! idx(0) is the number of determinants that DO NOT interact with key1
END_DOC
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,j,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze >= 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) &
+ popcnt( xor( key1(1,2,i), key2(1,2)))
if (degree_x2 > 4) then
idx(l) = i
l = l+1
else
cycle
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 > 4) then
idx(l) = i
l = l+1
else
cycle
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 > 4) then
idx(l) = i
l = l+1
else
cycle
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do j=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +&
popcnt(xor( key1(j,2,i), key2(j,2)))
if (degree_x2 > 4) then
idx(l) = i
l = l+1
endif
enddo
if (degree_x2 <= 5) then
exit
endif
enddo
endif
idx(0) = l-1
end
subroutine filter_connected(key1,key2,Nint,sze,idx)
use bitmasks
implicit none

44
src/Determinants/occ_pattern.irp.f
View File

@ -36,7 +36,7 @@ subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint)
amax -= popcnt( o(k,2) )
enddo
sze = int( min(binom_func(bmax, amax), 1.d8) )
sze = sze*sze
sze = 2*sze*sze + 16
end
@ -246,14 +246,22 @@ subroutine make_s2_eigenfunction
integer :: i,j,k
integer :: smax, s
integer(bit_kind), allocatable :: d(:,:,:), det_buffer(:,:,:)
integer :: N_det_new
integer :: N_det_new, ithread, omp_get_thread_num
integer, parameter :: bufsze = 1000
logical, external :: is_in_wavefunction
allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) )
smax = 1
N_det_new = 0
call write_int(6,N_occ_pattern,'Number of occupation patterns')
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(N_occ_pattern, psi_occ_pattern, elec_alpha_num,N_int) &
!$OMP PRIVATE(s,ithread, d, det_buffer, smax, N_det_new,i,j,k)
N_det_new = 0
call occ_pattern_to_dets_size(psi_occ_pattern(1,1,1),s,elec_alpha_num,N_int)
allocate (d(N_int,2,s), det_buffer(N_int,2,bufsze) )
smax = s
ithread=0
!$ ithread = omp_get_thread_num()
!$OMP DO
do i=1,N_occ_pattern
call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int)
s += 1
@ -270,40 +278,26 @@ subroutine make_s2_eigenfunction
det_buffer(k,1,N_det_new) = d(k,1,j)
det_buffer(k,2,N_det_new) = d(k,2,j)
enddo
! integer :: ne(2)
! ne(:) = 0
! do k=1,N_int
! ne(1) += popcnt(d(k,1,j))
! ne(2) += popcnt(d(k,2,j))
! enddo
! if (ne(1) /= elec_alpha_num) then
! call debug_det(d(1,1,j),N_int)
! stop "ALPHA"
! endif
! if (ne(2) /= elec_beta_num) then
! call debug_det(d(1,1,j),N_int)
! stop "BETA"
! endif
if (N_det_new == bufsze) then
call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,0)
call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread)
N_det_new = 0
endif
endif
enddo
enddo
!$OMP END DO NOWAIT
if (N_det_new > 0) then
call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,0)
! call fill_H_apply_buffer_no_selection_first_order_coef(N_det_new,det_buffer,N_int,0)
call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,ithread)
endif
!$OMP BARRIER
deallocate(d,det_buffer)
!$OMP END PARALLEL
call copy_H_apply_buffer_to_wf
SOFT_TOUCH N_det psi_coef psi_det
print *, 'Added determinants for S^2'
! logical :: found
! call remove_duplicates_in_psi_det(found)
call write_time(6)
end

46
src/Determinants/print_wf.irp.f
View File

@ -28,32 +28,32 @@ subroutine routine
if(degree == 0)then
print*,'Reference determinant '
else
call i_H_j(psi_det(1,1,i),psi_det(1,1,1),N_int,hij)
call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij)
call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*,'phase = ',phase
if(degree == 1)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
if(s1 == 1)then
norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
else
norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
endif
print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map)
double precision :: hmono,hdouble
call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
print*,'hmono = ',hmono
print*,'hdouble = ',hdouble
print*,'hmono+hdouble = ',hmono+hdouble
print*,'hij = ',hij
else
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
print*,'s2',s2
print*,'h2,p2 = ',h2,p2
print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
endif
! if(degree == 1)then
! print*,'s1',s1
! print*,'h1,p1 = ',h1,p1
! if(s1 == 1)then
! norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
! else
! norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
! endif
! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map)
! double precision :: hmono,hdouble
! call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
! print*,'hmono = ',hmono
! print*,'hdouble = ',hdouble
! print*,'hmono+hdouble = ',hmono+hdouble
! print*,'hij = ',hij
! else
! print*,'s1',s1
! print*,'h1,p1 = ',h1,p1
! print*,'s2',s2
! print*,'h2,p2 = ',h2,p2
! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
! endif
print*,'<Ref| H |D_I> = ',hij
endif

3
src/Determinants/s2.irp.f
View File

@ -223,13 +223,12 @@ subroutine S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
do istate=1,N_st
do i=n,1,-1
!$OMP ATOMIC
v_0(i,istate) = v_0(i,istate) + vt(i,istate)
enddo
enddo
!$OMP END CRITICAL
deallocate(vt)
!$OMP END PARALLEL

48
src/Determinants/slater_rules.irp.f
View File

@ -925,22 +925,29 @@ subroutine create_minilist(key_mask, fullList, miniList, idx_miniList, N_fullLis
N_miniList = 0
integer :: e_ab
e_ab = n_a+n_b
do i=1,N_fullList
e_a = n_a - popcnt(iand(fullList(1, 1, i), key_mask(1, 1)))
e_b = n_b - popcnt(iand(fullList(1, 2, i), key_mask(1, 2)))
e_a = e_ab - popcnt(iand(fullList(1, 1, i), key_mask(1, 1))) &
- popcnt(iand(fullList(1, 2, i), key_mask(1, 2)))
do ni=2,nint
e_a -= popcnt(iand(fullList(ni, 1, i), key_mask(ni, 1)))
e_b -= popcnt(iand(fullList(ni, 2, i), key_mask(ni, 2)))
e_a = e_a - popcnt(iand(fullList(ni, 1, i), key_mask(ni, 1))) &
- popcnt(iand(fullList(ni, 2, i), key_mask(ni, 2)))
end do
if(e_a + e_b <= 2) then
N_miniList = N_miniList + 1
do ni=1,Nint
miniList(ni,1,N_miniList) = fullList(ni,1,i)
miniList(ni,2,N_miniList) = fullList(ni,2,i)
enddo
idx_miniList(N_miniList) = i
end if
if(e_a > 2) then
cycle
endif
N_miniList = N_miniList + 1
miniList(1,1,N_miniList) = fullList(1,1,i)
miniList(1,2,N_miniList) = fullList(1,2,i)
do ni=2,Nint
miniList(ni,1,N_miniList) = fullList(ni,1,i)
miniList(ni,2,N_miniList) = fullList(ni,2,i)
enddo
idx_miniList(N_miniList) = i
end do
end subroutine
@ -1041,13 +1048,15 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
integer, allocatable :: idx(:)
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
allocate(idx(0:Ndet))
i_H_psi_array = 0.d0
call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
@ -1089,7 +1098,7 @@ subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
integer, allocatable :: idx(:)
BEGIN_DOC
! Computes <i|H|Psi> = \sum_J c_J <i|H|J>.
!
@ -1102,6 +1111,7 @@ subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
allocate(idx(0:Ndet))
i_H_psi_array = 0.d0
call filter_connected_i_H_psi0(keys,key,Nint,N_minilist,idx)
@ -1148,7 +1158,8 @@ subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet),n_interact
integer,allocatable :: idx(:)
integer :: n_interact
BEGIN_DOC
! <key|H|psi> for the various Nstates
END_DOC
@ -1158,6 +1169,7 @@ subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
allocate(idx(0:Ndet))
i_H_psi_array = 0.d0
call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
n_interact = 0
@ -1207,7 +1219,7 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
integer,allocatable :: idx(:)
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
@ -1215,6 +1227,7 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
i_H_psi_array = 0.d0
allocate(idx(0:Ndet))
call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat)
do ii=1,idx(0)
i = idx(ii)
@ -1254,7 +1267,7 @@ subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_a
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
integer,allocatable :: idx(:)
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
@ -1262,6 +1275,7 @@ subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_a
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
i_H_psi_array = 0.d0
allocate(idx(0:Ndet))
call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat)
print*,'--------'
do ii=1,idx(0)

48
src/Determinants/spindeterminants.irp.f
View File

@ -393,6 +393,8 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
BEGIN_DOC
! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
! D_a^t C D_b
!
! Rows are alpha determinants and columns are beta.
END_DOC
integer :: i,j,k, l
integer(bit_kind) :: tmp_det(N_int,2)
@ -421,10 +423,54 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
call isort(to_sort, iorder, N_det)
call iset_order(psi_bilinear_matrix_rows,iorder,N_det)
call iset_order(psi_bilinear_matrix_columns,iorder,N_det)
call dset_order(psi_bilinear_matrix_values,iorder,N_det)
do l=1,N_states
call dset_order(psi_bilinear_matrix_values(1,l),iorder,N_det)
enddo
deallocate(iorder,to_sort)
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ]
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows, (N_det) ]
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ]
use bitmasks
implicit none
BEGIN_DOC
! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
! D_a^t C D_b
!
! Rows are Beta determinants and columns are alpha
END_DOC
integer :: i,j,k,l
PROVIDE psi_coef_sorted_bit
integer, allocatable :: iorder(:), to_sort(:)
allocate(iorder(N_det), to_sort(N_det))
do l=1,N_states
do k=1,N_det
psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l)
enddo
enddo
do k=1,N_det
psi_bilinear_matrix_transp_columns(k) = psi_bilinear_matrix_columns(k)
psi_bilinear_matrix_transp_rows (k) = psi_bilinear_matrix_rows (k)
i = psi_bilinear_matrix_transp_columns(k)
j = psi_bilinear_matrix_transp_rows (k)
to_sort(k) = N_det_beta_unique * (j-1) + i
iorder(k) = k
enddo
call isort(to_sort, iorder, N_det)
call iset_order(psi_bilinear_matrix_transp_rows,iorder,N_det)
call iset_order(psi_bilinear_matrix_transp_columns,iorder,N_det)
do l=1,N_states
call dset_order(psi_bilinear_matrix_transp_values(1,l),iorder,N_det)
enddo
deallocate(iorder,to_sort)
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_det_beta_unique,N_states) ]
implicit none
BEGIN_DOC

13
src/Integrals_Bielec/ao_bi_integrals.irp.f
View File

@ -346,6 +346,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
integer :: n_integrals, rc
integer :: kk, m, j1, i1, lmax
character*(64) :: fmt
integral = ao_bielec_integral(1,1,1,1)
@ -365,14 +366,16 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
call cpu_time(cpu_1)
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
character*(32) :: task
call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals')
do l=ao_num,1,-1
write(task,*) "triangle ", l
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
character(len=:), allocatable :: task
allocate(character(len=ao_num*12) :: task)
write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
do l=1,ao_num
write(task,fmt) (i,l, i=1,l)
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task))
enddo
deallocate(task)
call zmq_set_running(zmq_to_qp_run_socket)

22
src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
View File

@ -57,12 +57,12 @@ subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value,
endif
! Activate is zmq_socket_push is a REQ
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
end
@ -187,11 +187,11 @@ subroutine ao_bielec_integrals_in_map_collector
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
! Activate if zmq_socket_pull is a REP
! rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
! if (rc /= 4) then
! print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
! stop 'error'
! endif
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
if (rc /= 4) then
print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
stop 'error'
endif
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)

2
src/Integrals_Bielec/mo_bi_integrals.irp.f
View File

@ -35,6 +35,8 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
print*, 'MO integrals provided'
return
else
PROVIDE ao_bielec_integrals_in_map
endif
if(no_vvvv_integrals)then

28
src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
View File

@ -51,6 +51,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
print*, 'Providing the nuclear electron pseudo integrals (local)'
call wall_time(wall_1)
wall_0 = wall_1
call cpu_time(cpu_1)
thread_num = 0
@ -102,7 +103,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
pseudo_n_k_transp (1,k), &
pseudo_dz_k_transp(1,k), &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
ao_pseudo_integral_local(i,j) = ao_pseudo_integral_local(i,j) +&
ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
@ -148,14 +148,9 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
print*, 'Providing the nuclear electron pseudo integrals (non-local)'
call wall_time(wall_1)
wall_0 = wall_1
call cpu_time(cpu_1)
thread_num = 0
!write(34,*) 'xxxNONLOCxxx'
!write(34,*) ' pseudo_lmax,pseudo_kmax', pseudo_lmax,pseudo_kmax
!write(34,*) ' pseudo_v_kl_transp(1,0,k)', pseudo_v_kl_transp
!write(34,*) ' pseudo_n_kl_transp(1,0,k)', pseudo_n_kl_transp
!write(34,*) ' pseudo_dz_kl_transp(1,0,k)', pseudo_dz_kl_transp
!write(34,*) 'xxxNONLOCxxx'
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -170,7 +165,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
!$ thread_num = omp_get_thread_num()
!$OMP DO SCHEDULE (guided)
!
do j = 1, ao_num
num_A = ao_nucl(j)
@ -207,15 +202,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
pseudo_n_kl_transp(1,0,k), &
pseudo_dz_kl_transp(1,0,k), &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
! write(34,*) i,j,k
! write(34,*) &
! A_center,power_A,alpha,B_center,power_B,beta,C_center, &
! Vpseudo(pseudo_lmax,pseudo_kmax, &
! pseudo_v_kl_transp(1,0,k), &
! pseudo_n_kl_transp(1,0,k), &
! pseudo_dz_kl_transp(1,0,k), &
! A_center,power_A,alpha,B_center,power_B,beta,C_center)
! write(34,*) ''
enddo
ao_pseudo_integral_non_local(i,j) = ao_pseudo_integral_non_local(i,j) +&
ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
@ -232,12 +218,12 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
endif
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_v_k_transp, (pseudo_klocmax,nucl_num) ]

46
src/Integrals_Monoelec/pseudopot.f90
View File

@ -15,14 +15,10 @@ double precision function Vps &
implicit none
integer n_a(3),n_b(3)
double precision g_a,g_b,a(3),b(3),c(3)
integer kmax_max,lmax_max
parameter (kmax_max=2,lmax_max=2)
integer lmax,kmax,n_kl(kmax_max,0:lmax_max)
double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max)
integer klocmax_max
parameter (klocmax_max=10)
integer klocmax,n_k(klocmax_max)
double precision v_k(klocmax_max),dz_k(klocmax_max)
integer lmax,kmax,n_kl(kmax,0:lmax)
double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
integer klocmax,n_k(klocmax)
double precision v_k(klocmax),dz_k(klocmax)
double precision Vloc,Vpseudo
Vps=Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) &
@ -36,13 +32,10 @@ double precision function Vps_num &
implicit none
integer n_a(3),n_b(3)
double precision g_a,g_b,a(3),b(3),c(3),rmax
integer kmax_max,lmax_max
parameter (kmax_max=2,lmax_max=2)
integer lmax,kmax,n_kl(kmax_max,0:lmax_max)
double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max)
integer klocmax_max;parameter (klocmax_max=10)
integer klocmax,n_k(klocmax_max)
double precision v_k(klocmax_max),dz_k(klocmax_max)
integer lmax,kmax,n_kl(kmax,0:lmax)
double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
integer klocmax,n_k(klocmax)
double precision v_k(klocmax),dz_k(klocmax)
double precision Vloc_num,Vpseudo_num,v1,v2
integer npts,nptsgrid
nptsgrid=50
@ -54,11 +47,9 @@ end
double precision function Vloc_num(npts_over,xmax,klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c)
implicit none
integer klocmax_max
parameter (klocmax_max=10)
integer klocmax
double precision v_k(klocmax_max),dz_k(klocmax_max)
integer n_k(klocmax_max)
double precision v_k(klocmax),dz_k(klocmax)
integer n_k(klocmax)
integer npts_over,ix,iy,iz
double precision xmax,dx,x,y,z
double precision a(3),b(3),c(3),term,r,orb_phi,g_a,g_b,ac(3),bc(3)
@ -705,12 +696,9 @@ end
double precision function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c)
implicit none
integer klocmax_max,lmax_max,ntot_max
parameter (klocmax_max=10,lmax_max=2)
parameter (ntot_max=10)
integer klocmax
double precision v_k(klocmax_max),dz_k(klocmax_max),crochet,bigA
integer n_k(klocmax_max)
double precision v_k(klocmax),dz_k(klocmax),crochet,bigA
integer n_k(klocmax)
double precision a(3),g_a,b(3),g_b,c(3),d(3)
integer n_a(3),n_b(3),ntotA,ntotB,ntot,m
integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p
@ -719,6 +707,7 @@ double precision,allocatable :: array_R_loc(:,:,:)
double precision,allocatable :: array_coefs(:,:,:,:,:,:)
double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
fourpi=4.d0*dacos(-1.d0)
f=fourpi**1.5d0
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
@ -755,8 +744,8 @@ double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
dreal=2.d0*d2
allocate (array_R_loc(-2:ntot_max+klocmax_max,klocmax_max,0:ntot_max))
allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max))
allocate (array_R_loc(-2:ntot+klocmax,klocmax,0:ntot))
allocate (array_coefs(0:ntot,0:ntot,0:ntot,0:ntot,0:ntot,0:ntot))
do ktot=-2,ntotA+ntotB+klocmax
do l=0,ntot
@ -2111,9 +2100,7 @@ end
! r : Distance between the Atomic Orbital center and the considered point
double precision function ylm_orb(l,m,c,a,n_a,g_a,r)
implicit none
integer lmax_max,ntot_max
parameter (lmax_max=2)
parameter (ntot_max=14)
integer lmax_max
integer l,m
double precision a(3),g_a,c(3)
double precision prod,binom_func,accu,bigI,ylm,bessel_mod
@ -2131,7 +2118,6 @@ factor=fourpi*dexp(-arg)
areal=2.d0*g_a*ac
ntotA=n_a(1)+n_a(2)+n_a(3)
if(ntotA.gt.ntot_max)stop 'increase ntot_max'
if(ac.eq.0.d0)then
ylm_orb=dsqrt(fourpi)*r**ntotA*dexp(-g_a*r**2)*bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))

8
src/MO_Basis/EZFIO.cfg
View File

@ -20,7 +20,13 @@ doc: MO occupation numbers
interface: ezfio
size: (mo_basis.mo_tot_num)
[mo_class]
type: character*(32)
doc: c: core, i: inactive, a: active, v: virtual, d: deleted
interface: ezfio, provider
size: (mo_basis.mo_tot_num)
[ao_md5]
type: character*(32)
doc: Ao_md5
interface: ezfio
interface: ezfio

80
src/MO_Basis/cholesky_mo.irp.f
View File

@ -1,8 +1,20 @@
subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
implicit none
BEGIN_DOC
! Cholesky decomposition of AO Density matrix to
! generate MOs
! Cholesky decomposition of AO Density matrix
!
! n : Number of AOs
! m : Number of MOs
!
! P(LDP,n) : Density matrix in AO basis
!
! C(LDC,m) : MOs
!
! tol_in : tolerance
!
! rank : Nomber of local MOs (output)
!
END_DOC
integer, intent(in) :: n,m, LDC, LDP
double precision, intent(in) :: P(LDP,n)
@ -15,9 +27,6 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
integer :: ipiv(n)
double precision:: tol
double precision, allocatable :: W(:,:), work(:)
!DEC$ ATTRIBUTES ALIGN: 32 :: W
!DEC$ ATTRIBUTES ALIGN: 32 :: work
!DEC$ ATTRIBUTES ALIGN: 32 :: ipiv
allocate(W(LDC,n),work(2*n))
tol=tol_in
@ -41,40 +50,37 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
deallocate(W,work)
end
BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ]
subroutine svd_mo(n,m,P,LDP,C,LDC)
implicit none
BEGIN_DOC
! Density matrix in MO basis
END_DOC
integer :: i,j,k
mo_density_matrix = 0.d0
do k=1,mo_tot_num
if (mo_occ(k) == 0.d0) then
cycle
endif
do j=1,ao_num
do i=1,ao_num
mo_density_matrix(i,j) = mo_density_matrix(i,j) + &
mo_occ(k) * mo_coef(i,k) * mo_coef(j,k)
enddo
enddo
enddo
END_PROVIDER
! Singular value decomposition of the AO Density matrix
!
! n : Number of AOs
BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ]
implicit none
BEGIN_DOC
! Density matrix in MO basis (virtual MOs)
! m : Number of MOs
!
! P(LDP,n) : Density matrix in AO basis
!
! C(LDC,m) : MOs
!
! tol_in : tolerance
!
! rank : Nomber of local MOs (output)
!
END_DOC
integer :: i,j,k
mo_density_matrix_virtual = 0.d0
do k=1,mo_tot_num
do j=1,ao_num
do i=1,ao_num
mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + &
(2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k)
enddo
enddo
enddo
END_PROVIDER
integer, intent(in) :: n,m, LDC, LDP
double precision, intent(in) :: P(LDP,n)
double precision, intent(out) :: C(LDC,m)
integer :: info
integer :: i,k
integer :: ipiv(n)
double precision:: tol
double precision, allocatable :: W(:,:), work(:)
allocate(W(LDC,n),work(2*n))
call svd(P,LDP,C,LDC,W,size(W,1),m,n)
deallocate(W,work)
end

1
src/MO_Basis/mos.irp.f
View File

@ -258,3 +258,4 @@ subroutine mix_mo_jk(j,k)
enddo
end

45
src/Nuclei/nuclei.irp.f
View File

@ -169,7 +169,7 @@ END_PROVIDER
'Nuclear repulsion energy')
END_PROVIDER
BEGIN_PROVIDER [ character*(128), element_name, (36)]
BEGIN_PROVIDER [ character*(128), element_name, (78)]
BEGIN_DOC
! Array of the name of element, sorted by nuclear charge (integer)
END_DOC
@ -209,4 +209,47 @@ BEGIN_PROVIDER [ character*(128), element_name, (36)]
element_name(34) = 'Se'
element_name(35) = 'Br'
element_name(36) = 'Kr'
element_name(37) = 'Rb'
element_name(38) = 'Sr'
element_name(39) = 'Y'
element_name(40) = 'Zr'
element_name(41) = 'Nb'
element_name(42) = 'Mo'
element_name(43) = 'Tc'
element_name(44) = 'Ru'
element_name(45) = 'Rh'
element_name(46) = 'Pd'
element_name(47) = 'Ag'
element_name(48) = 'Cd'
element_name(49) = 'In'
element_name(50) = 'Sn'
element_name(51) = 'Sb'
element_name(52) = 'Te'
element_name(53) = 'I'
element_name(54) = 'Xe'
element_name(55) = 'Cs'
element_name(56) = 'Ba'
element_name(57) = 'La'
element_name(58) = 'Ce'
element_name(59) = 'Pr'
element_name(60) = 'Nd'
element_name(61) = 'Pm'
element_name(62) = 'Sm'
element_name(63) = 'Eu'
element_name(64) = 'Gd'
element_name(65) = 'Tb'
element_name(66) = 'Dy'
element_name(67) = 'Ho'
element_name(68) = 'Er'
element_name(69) = 'Tm'
element_name(70) = 'Yb'
element_name(71) = 'Lu'
element_name(72) = 'Hf'
element_name(73) = 'Ta'
element_name(74) = 'W'
element_name(75) = 'Re'
element_name(76) = 'Os'
element_name(77) = 'Ir'
element_name(78) = 'Pt'
END_PROVIDER

4
src/Utils/LinearAlgebra.irp.f
View File

@ -152,8 +152,8 @@ subroutine ortho_qr(A,LDA,m,n)
LWORK=2*WORK(1)
deallocate(WORK)
allocate(WORK(LWORK))
call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
call dorgqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO)
call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
call dorgqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO)
deallocate(WORK,jpvt,tau)
end

21
src/Utils/map_functions.irp.f
View File

@ -73,10 +73,11 @@ subroutine map_load_from_disk(filename,map)
implicit none
character*(*), intent(in) :: filename
type(map_type), intent(inout) :: map
double precision :: x
type(c_ptr) :: c_pointer(3)
integer :: fd(3)
integer*8 :: i,k
integer :: n_elements
integer*8 :: i,k, l
integer :: n_elements, j
@ -95,7 +96,9 @@ subroutine map_load_from_disk(filename,map)
call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .True., c_pointer(3))
call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /))
l = 0_8
k = 1_8
x = 0.d0
do i=0_8, map % map_size
deallocate(map % map(i) % value)
deallocate(map % map(i) % key)
@ -106,9 +109,21 @@ subroutine map_load_from_disk(filename,map)
k = map % consolidated_idx (i+2)
map % map(i) % map_size = n_elements
map % map(i) % n_elements = n_elements
! Load memory from disk
do j=1,n_elements
x = x + map % map(i) % value(j)
l = iand(l,map % map(i) % key(j))
if (map % map(i) % value(j) > 1.e30) then
stop 'Error in integrals file'
endif
if (map % map(i) % key(j) < 0) then
stop 'Error in integrals file'
endif
enddo
enddo
map % sorted = x>0 .or. l == 0_8
map % n_elements = k-1
map % sorted = .True.
map % sorted = map % sorted .or. .True.
map % consolidated = .True.
end

74
src/ZMQ/utils.irp.f
View File

@ -94,7 +94,7 @@ subroutine switch_qp_run_to_master
print *, 'This run should be started with the qp_run command'
stop -1
endif
qp_run_address = trim(buffer)
qp_run_address = adjustl(buffer)
print *, 'Switched to qp_run master : ', trim(qp_run_address)
integer :: i
@ -235,8 +235,8 @@ function new_zmq_pull_socket()
if (zmq_context == 0_ZMQ_PTR) then
stop 'zmq_context is uninitialized'
endif
new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL)
! new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP)
! new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL)
new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP)
call omp_unset_lock(zmq_lock)
if (new_zmq_pull_socket == 0_ZMQ_PTR) then
stop 'Unable to create zmq pull socket'
@ -312,8 +312,8 @@ function new_zmq_push_socket(thread)
if (zmq_context == 0_ZMQ_PTR) then
stop 'zmq_context is uninitialized'
endif
new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_PUSH)
! new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
! new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_PUSH)
new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
call omp_unset_lock(zmq_lock)
if (new_zmq_push_socket == 0_ZMQ_PTR) then
stop 'Unable to create zmq push socket'
@ -684,10 +684,43 @@ subroutine add_task_to_taskserver(zmq_to_qp_run_socket,task)
character*(*), intent(in) :: task
integer :: rc, sze
character*(512) :: message
character(len=:), allocatable :: message
message='add_task '//trim(zmq_state)//' '//trim(task)
sze = len(message)
rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
if (rc /= sze) then
print *, rc, sze
print *, irp_here,': f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)'
stop 'error'
endif
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, sze-1, 0)
if (message(1:rc) /= 'ok') then
print *, trim(task)
print *, 'Unable to add the next task'
stop -1
endif
end
subroutine add_task_to_taskserver_send(zmq_to_qp_run_socket,task)
use f77_zmq
implicit none
BEGIN_DOC
! Get a task from the task server
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
character*(*), intent(in) :: task
integer :: rc, sze
character(len=:), allocatable :: message
sze = len(trim(task))+12+len(trim(zmq_state))
message = repeat(' ',sze)
write(message,*) 'add_task '//trim(zmq_state)//' '//trim(task)
sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
if (rc /= sze) then
print *, rc, sze
@ -695,10 +728,20 @@ subroutine add_task_to_taskserver(zmq_to_qp_run_socket,task)
stop 'error'
endif
end
subroutine add_task_to_taskserver_recv(zmq_to_qp_run_socket)
use f77_zmq
implicit none
BEGIN_DOC
! Get a task from the task server
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer :: rc, sze
character*(512) :: message
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
message = trim(message(1:rc))
if (trim(message) /= 'ok') then
print *, trim(task)
if (message(1:rc) /= 'ok') then
print *, 'Unable to add the next task'
stop -1
endif
@ -726,8 +769,7 @@ subroutine task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_id)
endif
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
message = trim(message(1:rc))
if (trim(message) /= 'ok') then
if (trim(message(1:rc)) /= 'ok') then
print *, 'Unable to send task_done message'
stop -1
endif
@ -752,17 +794,17 @@ subroutine get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id,task)
write(message,*) 'get_task '//trim(zmq_state), worker_id
sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
if (rc /= sze) then
print *, irp_here, ':f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)'
stop 'error'
endif
message = repeat(' ',512)
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
message = trim(message(1:rc))
read(message,*) reply
read(message(1:rc),*) reply
if (trim(reply) == 'get_task_reply') then
read(message,*) reply, task_id
read(message(1:rc),*) reply, task_id
rc = 15
do while (message(rc:rc) == ' ')
rc += 1

4
tests/bats/cassd.bats
View File

@ -15,12 +15,12 @@ source $QP_ROOT/tests/bats/common.bats.sh
energy="$(ezfio get cas_sd_zmq energy_pt2)"
eq $energy -76.231084536315 5.E-5
ezfio set determinants n_det_max 2048
ezfio set determinants n_det_max 1024
ezfio set determinants read_wf True
ezfio set perturbation do_pt2_end True
qp_run cassd_zmq $INPUT
ezfio set determinants read_wf False
energy="$(ezfio get cas_sd_zmq energy)"
eq $energy -76.2300887947446 2.E-5
eq $energy -76.2225863580749 2.E-5
}

8
tests/bats/mrcepa0.bats
View File

@ -16,7 +16,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy_pt2)"
eq $energy -76.23752746236 1.e-4
eq $energy -76.2382106224545 1.e-4
}
@test "MRCC H2O cc-pVDZ" {
@ -32,7 +32,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy_pt2)"
eq $energy -76.237469267705 2.e-4
eq $energy -76.2381673136696 2.e-4
}
@test "MRSC2 H2O cc-pVDZ" {
@ -48,7 +48,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy_pt2)"
eq $energy -76.2347764009137 2.e-4
eq $energy -76.235786994991 2.e-4
}
@test "MRCEPA0 H2O cc-pVDZ" {
@ -64,6 +64,6 @@ source $QP_ROOT/tests/bats/common.bats.sh
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy_pt2)"
eq $energy -76.2406942855164 2.e-4
eq $energy -76.2411829210128 2.e-4
}