Merge pull request #187 from scemama/master

Accelerated Davidson

README.md

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

config/gfortran_debug.cfg

@@ -51,7 +51,7 @@ FCFLAGS : -Ofast
 # -g           : Extra debugging information
 #
 [DEBUG]
-FCFLAGS : -g -msse4.2 
+FCFLAGS : -g -msse4.2  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant 
 
 # OpenMP flags
 #################

ocaml/Basis.ml

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

ocaml/Basis.mli

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

ocaml/Input_determinants_by_hand.ml

@@ -7,6 +7,7 @@ module Determinants_by_hand : sig
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
+      n_states               : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -18,11 +19,14 @@ module Determinants_by_hand : sig
   val to_rst : t -> Rst_string.t
   val of_rst : Rst_string.t -> t option
   val read_n_int : unit -> N_int_number.t
+  val update_ndet : Det_number.t -> unit
+  val extract_state : States_number.t -> unit
 end = struct
   type t = 
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
+      n_states               : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -129,12 +133,12 @@ end = struct
     |> Array.map ~f:Det_coef.of_float
   ;;
 
-  let write_psi_coef ~n_det c =
+  let write_psi_coef ~n_det ~n_states c =
     let n_det = Det_number.to_int n_det
     and c = Array.to_list c
             |> List.map ~f:Det_coef.to_float
     and n_states = 
-      read_n_states () |> States_number.to_int
+      States_number.to_int n_states
     in
     Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
     |> Ezfio.set_determinants_psi_coef 
@@ -200,6 +204,7 @@ end = struct
           expected_s2            = read_expected_s2 ()          ;
           psi_coef               = read_psi_coef ()             ;
           psi_det                = read_psi_det ()              ;
+          n_states               = read_n_states ()             ;
         }
     else
       failwith "No molecular orbitals, so no determinants"
@@ -222,12 +227,14 @@ end = struct
               expected_s2          ;
               psi_coef             ;
               psi_det              ;
+              n_states             ;
             } =
      write_n_int n_int ;
      write_bit_kind bit_kind;
      write_n_det n_det;
+     write_n_states n_states;
      write_expected_s2 expected_s2;
-     write_psi_coef ~n_det:n_det psi_coef ;
+     write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
      write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
   ;;
 
@@ -298,6 +305,7 @@ Determinants ::
 n_int                  = %s
 bit_kind               = %s
 n_det                  = %s
+n_states               = %s
 expected_s2            = %s
 psi_coef               = %s
 psi_det                = %s
@@ -305,6 +313,7 @@ psi_det                = %s
      (b.n_int         |> N_int_number.to_string)
      (b.bit_kind      |> Bit_kind.to_string)
      (b.n_det         |> Det_number.to_string)
+     (b.n_states      |> States_number.to_string)
      (b.expected_s2   |> Positive_float.to_string)
      (b.psi_coef  |> Array.to_list |> List.map ~f:Det_coef.to_string
       |> String.concat ~sep:", ")
@@ -433,14 +442,83 @@ psi_det                = %s
       |> Bit_kind.to_int)
     and n_int =
       Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
+    and n_states =
+      Printf.sprintf "(n_states %d)" (States_number.to_int @@ read_n_states ())
     in
     let s = 
-       String.concat [ header ; bitkind ; n_int ; psi_coef ; psi_det]
+       String.concat [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det]
     in
 
+
+
+
     Generic_input_of_rst.evaluate_sexp t_of_sexp s
   ;;
 
+  let update_ndet n_det_new =
+    Printf.printf "Reducing n_det to %d\n" (Det_number.to_int n_det_new);
+    let n_det_new =
+      Det_number.to_int n_det_new
+    in
+    let det =
+      read ()
+    in
+    let n_det_old, n_states =
+      Det_number.to_int det.n_det,
+      States_number.to_int det.n_states
+    in
+    if n_det_new = n_det_old then
+      ()
+    ;
+    if n_det_new > n_det_new then
+      failwith @@ Printf.sprintf "Requested n_det should be less than %d" n_det_old
+    ;
+    for j=0 to (n_states-1) do
+      let ishift_old, ishift_new =
+        j*n_det_old,
+        j*n_det_new
+      in
+      for i=0 to (n_det_new-1) do
+        det.psi_coef.(i+ishift_new) <- det.psi_coef.(i+ishift_old)
+      done
+    done
+    ;
+    let new_det =
+      { det with n_det = (Det_number.of_int n_det_new) }
+    in
+    write new_det
+  ;;
+
+  let extract_state istate =
+    Printf.printf "Extracting state %d\n" (States_number.to_int istate);
+    let det =
+      read ()
+    in
+    let n_det, n_states =
+      Det_number.to_int det.n_det,
+      States_number.to_int det.n_states
+    in
+    if (States_number.to_int istate) > n_states then
+      failwith "State to extract should not be greater than n_states"
+    ;
+    let j =
+      (States_number.to_int istate) - 1
+    in
+    begin
+      if (j>0) then
+        let ishift =
+          j*n_det
+        in
+        for i=0 to (n_det-1) do
+          det.psi_coef.(i) <- det.psi_coef.(i+ishift)
+        done
+    end;
+    let new_det =
+      { det with n_states = (States_number.of_int 1) }
+    in
+    write new_det
+  ;;
+
 end
 
 

ocaml/Makefile

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

ocaml/Message.ml

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

ocaml/Message_lexer.mll

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

ocaml/Symmetry.ml

@@ -85,7 +85,7 @@ module Xyz = struct
   let of_string s =
     let flush state accu number =
       let n = 
-        if (number = "") then 0
+        if (number = "") then 1
         else (Int.of_string number) 
       in
       match state with

ocaml/TaskServer.ml

@@ -47,6 +47,14 @@ let debug str =
 let zmq_context =
   ZMQ.Context.create ()
 
+let () =
+  let nproc =
+    match Sys.getenv "OMP_NUM_THREADS" with
+    | Some m -> int_of_string m
+    | None -> 2
+  in
+  ZMQ.Context.set_io_threads zmq_context nproc
+
 
 let bind_socket ~socket_type ~socket ~port =
   let rec loop = function
@@ -62,7 +70,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 +115,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 +137,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 +305,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 +318,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 +346,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 +356,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 +423,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 +439,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;

ocaml/create_git_sha1.sh

@@ -2,7 +2,7 @@
 
 SHA1=$(git log -1  | head -1 | cut -d ' ' -f 2)
 DATE=$(git log -1  | grep Date | cut -d ':' -f 2-)
-MESSAGE=$(git log -1  | tail -1)
+MESSAGE=$(git log -1  | tail -1 | sed 's/"/\\"/g')
 cat << EOF > Git.ml
 open Core.Std
 let sha1 = "$SHA1" |> String.strip

ocaml/qp_create_ezfio_from_xyz.ml

@@ -21,6 +21,9 @@ let spec =
      ~doc:" Compute AOs in the Cartesian basis set (6d, 10f, ...)"
   +> anon ("(xyz_file|zmt_file)" %: file )
 
+type element =
+| Element of Element.t
+| Int_elem of (Nucl_number.t * Element.t)
 
 (** Handle dummy atoms placed on bonds *)
 let dummy_centers ~threshold ~molecule ~nuclei =
@@ -115,17 +118,14 @@ let run ?o b c d m p cart xyz_file =
   (* Open basis set channels *)
   let basis_channel element =
     let key =
-      Element.to_string element
+      match element with
+      | Element e -> Element.to_string e
+      | Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i)  (Element.to_string e)
     in
     match Hashtbl.find basis_table key with
     | Some in_channel -> 
         in_channel
-    | None -> 
-         let msg = 
-           Printf.sprintf "%s is not defined in basis %s.%!"
-           (Element.to_long_string element) b ;
-         in
-         failwith msg
+    | None -> raise Not_found
   in
   
   let temp_filename =
@@ -189,12 +189,21 @@ let run ?o b c d m p cart xyz_file =
       | Some (key, basis) -> (*Aux basis *)
         begin
           let elem  = 
-            Element.of_string key
+            try
+              Element (Element.of_string key)
+            with Element.ElementError _ ->
+              let result = 
+                match (String.split ~on:',' key) with
+                | i :: k :: [] -> (Nucl_number.of_int @@ int_of_string i, Element.of_string k)
+                | _ -> failwith "Expected format is int,Element:basis"
+              in Int_elem result
           and basis =
             String.lowercase basis
           in
           let key = 
-             Element.to_string elem
+             match elem with
+             | Element e -> Element.to_string e
+             | Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i) (Element.to_string e)
           in
           let new_channel =
             fetch_channel basis
@@ -202,7 +211,13 @@ let run ?o b c d m p cart xyz_file =
           begin
             match Hashtbl.add basis_table ~key:key ~data:new_channel with
             | `Ok -> ()
-            | `Duplicate -> failwith ("Duplicate definition of basis for "^(Element.to_long_string elem))
+            | `Duplicate ->
+              let e = 
+                match elem with
+                | Element e -> e
+                | Int_elem (_,e) -> e
+              in
+              failwith ("Duplicate definition of basis for "^(Element.to_long_string e))
           end
        end
     end;
@@ -537,7 +552,20 @@ let run ?o b c d m p cart xyz_file =
               | Element.X -> Element.H
               | e -> e
             in
-            Basis.read_element (basis_channel x.Atom.element) i e
+            let key =
+              Int_elem (i,x.Atom.element)
+            in
+            try
+              Basis.read_element (basis_channel key) i e
+            with Not_found ->
+              let key =
+                Element x.Atom.element
+              in
+              try
+                Basis.read_element (basis_channel key) i e
+              with Not_found ->
+                failwith (Printf.sprintf "Basis not found for atom %d (%s)" (Nucl_number.to_int i)
+                 (Element.to_string x.Atom.element) )
           with
           | End_of_file -> failwith
                 ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
@@ -647,6 +675,7 @@ atoms are taken from the same basis set, otherwise specific elements can be
 defined as follows:
 
  -b \"cc-pcvdz | H:cc-pvdz | C:6-31g\"
+ -b \"cc-pvtz | 1,H:sto-3g | 3,H:6-31g\"
 
 If a file with the same name as the basis set exists, this file will be read.
 Otherwise, the basis set is obtained from the database.

ocaml/qptypes_generator.ml

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

plugins/Alavi/.gitignore

@@ -1,23 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Utils
-alavi_graph
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/Alavi/README.rst

@@ -1,23 +0,0 @@
-=====
-alavi
-=====
-
-Documentation
-=============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-`alavi_graph <http://github.com/LCPQ/quantum_package/tree/master/src/Alavi/alavi_graph.irp.f#L1>`_
-  Undocumented
-
-Needed Modules
-==============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-.. image:: tree_dependency.png
-
-* `Determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants>`_
-

plugins/Alavi/alavi_graph.irp.f

@@ -1,28 +0,0 @@
-program alavi_graph
-  implicit none
-  integer :: exc(0:2,2,2),h1,p1,h2,p2,s1,s2
-  double precision :: phase
-
-  read_wf = .True.
-  touch read_wf
-
-  integer :: k,degree
-  double precision :: hii
-
-  do k=1,N_det
-      call get_excitation_degree(psi_det(1,1,1),psi_det(1,1,k),degree,N_int)
-      call  i_H_j(psi_det(1,1,k),psi_det(1,1,k),N_int,hii) 
-      print*, k,abs(psi_coef(k,1)), hii,degree
-
-!      if (degree == 2) then
-!      	call get_excitation(psi_det(1,1,1),psi_det(1,1,k),exc,degree,phase,N_int)
-!      	call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-!		print*, h1,h2,hii, abs(psi_coef(k,1))
-!      endif
-!      
-
-
-  enddo
-end
-
-!plot "test.dat" u (abs($2)):(abs($3)):4  w p palette

plugins/All_singles/.gitignore

@@ -1,5 +0,0 @@
-IRPF90_temp/
-IRPF90_man/
-irpf90.make
-irpf90_entities
-tags

plugins/CAS_SD/.gitignore

@@ -1,34 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Davidson
-Determinants
-Electrons
-Ezfio_files
-Generators_CAS
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-Utils
-ZMQ
-cas_s
-cas_s_selected
-cas_sd
-cas_sd_selected
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/CAS_SD_ZMQ/EZFIO.cfg

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

plugins/CAS_SD_ZMQ/cassd_zmq.irp.f

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

plugins/CAS_SD_ZMQ/ezfio_interface.irp.f

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

plugins/CAS_SD_ZMQ/run_selection_slave.irp.f

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

plugins/CAS_SD_ZMQ/target_pt2_ratio_cassd.irp.f

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

plugins/CID/.gitignore

@@ -1,28 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-cid
-cid_lapack
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/CID_SC2_selected/.gitignore

@@ -1,31 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-CID
-CID_selected
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-cid_sc2_selection
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/CID_selected/.gitignore

@@ -1,30 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-CID
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-cid_selection
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/CIS/.gitignore

@@ -1,28 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-cis
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-super_ci
-tags

plugins/CISD/.gitignore

@@ -1,29 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-ZMQ
-cisd
-cisd_lapack
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/CISD_SC2_selected/.gitignore

@@ -1,31 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-CISD
-CISD_selected
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-cisd_sc2_selection
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/CISD_selected/.gitignore

@@ -1,31 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-CISD
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-ZMQ
-cisd_selection
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/Casino/.gitignore

@@ -1,23 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Utils
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-save_for_casino
-tags

plugins/Casino/save_for_casino.irp.f

@@ -5,7 +5,7 @@ subroutine save_casino
  integer                        :: getUnitAndOpen, iunit
  integer, allocatable           :: itmp(:)
  integer                        :: n_ao_new
- real, allocatable              :: rtmp(:)
+ double precision, allocatable              :: rtmp(:)
  PROVIDE ezfio_filename
 
  iunit = getUnitAndOpen('gwfn.data','w')

plugins/DDCI_selected/.gitignore

@@ -1,29 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Generators_CAS
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-Utils
-ddci
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/DensityMatrix/.gitignore

@@ -1,13 +0,0 @@
-#
-# Do not modify this file. Add your ignored files to the gitignore
-# (without the dot at the beginning) file.
-#
-IRPF90_temp
-IRPF90_man
-irpf90.make
-tags
-Makefile.depend
-irpf90_entities
-build.ninja
-.ninja_log
-.ninja_deps

plugins/FCIdump/.gitignore

@@ -1,24 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Utils
-ZMQ
-ezfio_interface.irp.f
-fcidump
-irpf90.make
-irpf90_entities
-tags

plugins/Full_CI/.gitignore

@@ -1,34 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Davidson
-Determinants
-Electrons
-Ezfio_files
-Generators_full
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-Utils
-ZMQ
-ezfio_interface.irp.f
-full_ci
-full_ci_no_skip
-irpf90.make
-irpf90_entities
-tags
-target_pt2
-var_pt2_ratio

plugins/Full_CI/full_ci.irp.f

@@ -3,6 +3,11 @@ program full_ci
   integer                        :: i,k
 
   
+  print *,  '===================================================================='
+  print *,  'This program is slow. Consider using the Full_CI_ZMQ module instead.'
+  print *,  '===================================================================='
+  call sleep(2)
+
   double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
   integer                        :: N_st, degree
   N_st = N_states

plugins/Full_CI_ZMQ/.gitignore

@@ -1,5 +0,0 @@
-IRPF90_temp/
-IRPF90_man/
-irpf90.make
-irpf90_entities
-tags

plugins/Full_CI_ZMQ/energy.irp.f

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

plugins/Full_CI_ZMQ/fci_zmq.irp.f

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

plugins/Full_CI_ZMQ/pt2_slave.irp.f

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

plugins/Full_CI_ZMQ/pt2_stoch.irp.f

@@ -0,0 +1,38 @@
+program pt2_stoch
+  implicit none
+  read_wf = .True.
+  SOFT_TOUCH 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, relative_error
+
+  allocate (pt2(N_states))
+  pt2 = 0.d0
+  
+  E_CI_before = pt2_E0_denominator(1) + nuclear_repulsion
+  threshold_selectors = 1.d0
+  threshold_generators = 1d0 
+  relative_error = 1.d-3
+  call ZMQ_pt2(pt2, relative_error)
+  print *,  'Final step'
+  print *,  'N_det    = ', N_det
+  print *,  'PT2      = ', pt2
+  print *,  'E        = ', E_CI_before
+  print *,  'E+PT2    = ', E_CI_before+pt2
+  print *,  '-----'
+  call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before+pt2(1))
+end
+
+

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -0,0 +1,580 @@
+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, k, 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 pt2_weight
+
+  !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,1X,I9,''|'')') 0, tbc(i)
+        ipos += 20
+        if (ipos > 63980) 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,1X,I9,''|'')') j, tbc(i)
+          ipos += 20
+          if (ipos > 63980) 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 = 0.d0
+       do i=1,N_det_generators
+         do k=1,N_states
+            pt2(k) = pt2(k) + pt2_detail(k,i)
+         enddo
+       enddo
+    endif
+
+    tbc(0) = 0
+    if (pt2(1) /= 0.d0) then
+      exit
+    endif
+  end do
+
+  deallocate(pt2_detail, comb, computed, tbc)
+
+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) )
+  do i=1,N_det_generators
+    actually_computed(i) = computed(i)
+  enddo
+  
+  parts_to_get(:) = 1
+  if(fragment_first > 0) then
+    do i=1,fragment_first
+      parts_to_get(i) = fragment_count
+    enddo
+  endif
+
+  do i=1,tbc(0)
+    actually_computed(tbc(i)) = .false.
+  end do
+  
+  orgTBDcomb = int(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(1:N_states, index(i)) += pt2_mwen(1:N_states,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 = int(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
+      else
+        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)
+      endif
+    end if
+  end do pullLoop
+
+  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
+  integer, intent(inout) :: Ncomb
+  double precision, intent(out) :: comb(Ncomb)
+  integer, intent(inout) :: tbc(0:size_tbc)
+  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 ((i>1000).and.(icount > n)) then
+        call get_filling_teeth(computed, tbc)
+        icount = 0
+        n = ishft(tbc_save,-4)
+      endif
+  enddo
+  call get_filling_teeth(computed, tbc)
+
+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
+  
+  do i=1,N_det_generators
+    pt2_weight(i)  = pt2_weight(i) / pt2_cweight(N_det_generators)
+    pt2_cweight(i) = pt2_cweight(i) / pt2_cweight(N_det_generators)
+  enddo
+  
+  norm_left = 1d0
+  
+  comb_step = 1d0/dfloat(comb_teeth)
+  first_det_of_comb = 1
+  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
+
+
+
+
+
+

plugins/Full_CI_ZMQ/run_pt2_slave.irp.f

@@ -0,0 +1,172 @@
+
+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_generators))
+  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_generators)
+  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_generators)
+  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, 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)
+
+  do i=N+1,N_det_generators
+    pt2_detail(1:N_states,i) = 0.d0
+  enddo
+end subroutine
+ 
+ 
+BEGIN_PROVIDER [ double precision, pt2_workload, (N_det_generators) ]
+  integer :: i
+  do i=1,N_det_generators
+    pt2_workload(i) = dfloat(N_det_generators - i + 1)**2
+  end do
+  pt2_workload = pt2_workload / sum(pt2_workload)
+END_PROVIDER
+            

plugins/Full_CI_ZMQ/run_selection_slave.irp.f

@@ -26,7 +26,6 @@ subroutine run_selection_slave(thread,iproc,energy)
   call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
   if(worker_id == -1) then
     print *, "WORKER -1"
-    !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)
     return
@@ -41,20 +40,16 @@ 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)
-        call create_selection_buffer(N, N*3, buf2)
+        call create_selection_buffer(N, N*2, buf2)
       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
@@ -66,8 +61,10 @@ subroutine run_selection_slave(thread,iproc,energy)
         call push_selection_results(zmq_socket_push, pt2, buf, task_id(1), ctask)
         do i=1,buf%cur
           call add_to_selection_buffer(buf2, buf%det(1,1,i), buf%val(i))
+          if (buf2%cur == buf2%N) then
+            call sort_selection_buffer(buf2)
+          endif
         enddo
-        call sort_selection_buffer(buf2)
         buf%mini = buf2%mini
         pt2 = 0d0
         buf%cur = 0
@@ -115,7 +112,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 +146,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
  
  

plugins/Full_CI_ZMQ/selection.irp.f

@@ -1,5 +1,13 @@
 use bitmasks
 
+BEGIN_PROVIDER [ integer, fragment_count ]
+  implicit none
+  BEGIN_DOC
+  ! Number of fragments for the deterministic part
+  END_DOC
+  fragment_count = (elec_alpha_num-n_core_orb)**2
+END_PROVIDER
+
 
 double precision function integral8(i,j,k,l)
   implicit none
@@ -39,7 +47,7 @@ subroutine assert(cond, msg)
   logical, intent(in) :: cond
   
   if(.not. cond) then
-    print *, "assert fail: "//msg
+    print *, "assert failed: "//msg
     stop
   end if
 end subroutine
@@ -50,7 +58,7 @@ subroutine get_mask_phase(det, phasemask)
   implicit none
 
   integer(bit_kind), intent(in) :: det(N_int, 2)
-  integer(1), intent(out) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1), intent(out) :: phasemask(2,N_int*bit_kind_size)
   integer :: s, ni, i
   logical :: change
 
@@ -60,18 +68,18 @@ subroutine get_mask_phase(det, phasemask)
     do ni=1,N_int
       do i=0,bit_kind_size-1
         if(BTEST(det(ni, s), i)) change = .not. change
-          if(change) phasemask((ni-1)*bit_kind_size + i + 1, s) = 1_1
+          if(change) phasemask(s, (ni-1)*bit_kind_size + i + 1) = 1_1
       end do
     end do
   end do
 end subroutine
 
 
-subroutine select_connected(i_generator,E0,pt2,b)
+subroutine select_connected(i_generator,E0,pt2,b,subset)
   use bitmasks
   use selection_types
   implicit none
-  integer, intent(in)            :: i_generator
+  integer, intent(in)            :: i_generator, subset
   type(selection_buffer), intent(inout) :: b
   double precision, intent(inout)  :: pt2(N_states)
   integer :: k,l
@@ -90,8 +98,7 @@ subroutine select_connected(i_generator,E0,pt2,b)
       particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) )
 
     enddo
-    call select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b)
-    call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b)
+    call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b,subset)
   enddo
 end subroutine
 
@@ -100,186 +107,30 @@ double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2)
   use bitmasks
   implicit none
 
-  integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1), intent(in) :: phasemask(2,*)
   integer, intent(in) :: s1, s2, h1, h2, p1, p2
   logical :: change
-  integer(1) :: np
-  double precision, parameter :: res(0:1) = (/1d0, -1d0/)
+  integer(1) :: np1
+  integer :: np
+  double precision, save :: res(0:1) = (/1d0, -1d0/)
 
-  np = phasemask(h1,s1) + phasemask(p1,s1) + phasemask(h2,s2) + phasemask(p2,s2)
-  if(p1 < h1) np = np + 1_1
-  if(p2 < h2) np = np + 1_1
+  np1 = phasemask(s1,h1) + phasemask(s1,p1) + phasemask(s2,h2) + phasemask(s2,p2)
+  np = np1
+  if(p1 < h1) np = np + 1
+  if(p2 < h2) np = np + 1
   
-  if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1_1
-  get_phase_bi = res(iand(np,1_1))
-end function
+  if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1
+  get_phase_bi = res(iand(np,1))
+end 
 
 
 
-! Selection single
-! ----------------
-
-subroutine select_singles(i_gen,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf)
-  use bitmasks
-  use selection_types
-  implicit none
-  BEGIN_DOC
-! Select determinants connected to i_det by H
-  END_DOC
-  integer, intent(in)             :: i_gen
-  integer(bit_kind), intent(in)   :: hole_mask(N_int,2), particle_mask(N_int,2)
-  double precision, intent(in)    :: fock_diag_tmp(mo_tot_num)
-  double precision, intent(in)    :: E0(N_states)
-  double precision, intent(inout) :: pt2(N_states)
-  type(selection_buffer), intent(inout) :: buf
-  
-  double precision                :: vect(N_states, mo_tot_num)
-  logical                         :: bannedOrb(mo_tot_num)
-  integer                         :: i, j, k
-  integer                         :: h1,h2,s1,s2,i1,i2,ib,sp
-  integer(bit_kind)               :: hole(N_int,2), particle(N_int,2), mask(N_int, 2)
-  logical                         :: fullMatch, ok
-  
-  
-  do k=1,N_int
-    hole    (k,1) = iand(psi_det_generators(k,1,i_gen), hole_mask(k,1))
-    hole    (k,2) = iand(psi_det_generators(k,2,i_gen), hole_mask(k,2))
-    particle(k,1) = iand(not(psi_det_generators(k,1,i_gen)), particle_mask(k,1))
-    particle(k,2) = iand(not(psi_det_generators(k,2,i_gen)), particle_mask(k,2))
-  enddo
-
-  ! Create lists of holes and particles
-  ! -----------------------------------
-
-  integer                        :: N_holes(2), N_particles(2)
-  integer                        :: hole_list(N_int*bit_kind_size,2)
-  integer                        :: particle_list(N_int*bit_kind_size,2)
-
-  call bitstring_to_list_ab(hole    , hole_list    , N_holes    , N_int)
-  call bitstring_to_list_ab(particle, particle_list, N_particles, N_int)
-
-  do sp=1,2
-    do i=1, N_holes(sp)
-      h1 = hole_list(i,sp)
-      call apply_hole(psi_det_generators(1,1,i_gen), sp, h1, mask, ok, N_int)
-      bannedOrb = .true.
-      do j=1,N_particles(sp)
-        bannedOrb(particle_list(j, sp)) = .false.
-      end do
-      call spot_hasBeen(mask, sp, psi_det_sorted, i_gen, N_det, bannedOrb, fullMatch)
-      if(fullMatch) cycle
-      vect = 0d0
-      call splash_p(mask, sp, psi_selectors(1,1,i_gen), psi_phasemask(1,1,i_gen), psi_selectors_coef_transp(1,i_gen), N_det_selectors - i_gen + 1, bannedOrb, vect)
-      call fill_buffer_single(i_gen, sp, h1, bannedOrb, fock_diag_tmp, E0, pt2, vect, buf)
-    end do
-  enddo
-end subroutine
-
-
-subroutine fill_buffer_single(i_generator, sp, h1, bannedOrb, fock_diag_tmp, E0, pt2, vect, buf)
-  use bitmasks
-  use selection_types
-  implicit none
-  
-  integer, intent(in) :: i_generator, sp, h1
-  double precision, intent(in) :: vect(N_states, mo_tot_num)
-  logical, intent(in) :: bannedOrb(mo_tot_num)
-  double precision, intent(in)           :: fock_diag_tmp(mo_tot_num)
-  double precision, intent(in)    :: E0(N_states)
-  double precision, intent(inout) :: pt2(N_states)
-  type(selection_buffer), intent(inout) :: buf
-  logical :: ok
-  integer :: s1, s2, p1, p2, ib, istate
-  integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
-  double precision :: e_pert, delta_E, val, Hii, max_e_pert, tmp
-  double precision, external :: diag_H_mat_elem_fock
-  
-  
-  call apply_hole(psi_det_generators(1,1,i_generator), sp, h1, mask, ok, N_int)
-  
-  do p1=1,mo_tot_num
-    if(bannedOrb(p1)) cycle
-    if(vect(1, p1) == 0d0) cycle
-    call apply_particle(mask, sp, p1, det, ok, N_int)
-    
-    
-    Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
-    max_e_pert = 0d0
-    
-    do istate=1,N_states
-      val = vect(istate, p1) + vect(istate, p1)
-      delta_E = E0(istate) - Hii
-      tmp = dsqrt(delta_E * delta_E + val * val)
-      if (delta_E < 0.d0) then
-        tmp = -tmp
-      endif
-      e_pert = 0.5d0 * ( tmp - delta_E)
-      pt2(istate) += e_pert
-      if(dabs(e_pert) > dabs(max_e_pert)) max_e_pert = e_pert
-    end do
-    
-    if(dabs(max_e_pert) > buf%mini) call add_to_selection_buffer(buf, det, max_e_pert)
-  end do
-end subroutine
-
-
-subroutine splash_p(mask, sp, det, phasemask, coefs, N_sel, bannedOrb, vect)
-  use bitmasks
-  implicit none
-
-  integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int,2,N_sel)
-  integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2, N_sel)
-  double precision, intent(in) :: coefs(N_states, N_sel)
-  integer, intent(in) :: sp, N_sel
-  logical, intent(inout) :: bannedOrb(mo_tot_num)
-  double precision, intent(inout)     :: vect(N_states, mo_tot_num)
-
-  integer :: i, j, h(0:2,2), p(0:3,2), nt
-  integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
-  
-  do i=1,N_int
-    negMask(i,1) = not(mask(i,1))
-    negMask(i,2) = not(mask(i,2))
-  end do
-
-  do i=1, N_sel
-    nt = 0
-    do j=1,N_int
-      mobMask(j,1) = iand(negMask(j,1), det(j,1,i))
-      mobMask(j,2) = iand(negMask(j,2), det(j,2,i))
-      nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
-    end do
-
-    if(nt > 3) cycle
-      
-    do j=1,N_int
-      perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
-      perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
-    end do
-
-    call bitstring_to_list(perMask(1,1), h(1,1), h(0,1), N_int)
-    call bitstring_to_list(perMask(1,2), h(1,2), h(0,2), N_int)
-
-    call bitstring_to_list(mobMask(1,1), p(1,1), p(0,1), N_int)
-    call bitstring_to_list(mobMask(1,2), p(1,2), p(0,2), N_int)
-    
-    if(nt == 3) then
-      call get_m2(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i))
-    else if(nt == 2) then
-      call get_m1(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i))
-    else
-      call get_m0(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i))
-    end if
-  end do
-end subroutine
-
-
 subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   use bitmasks
   implicit none
   
   integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
-  integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
   logical, intent(in) :: bannedOrb(mo_tot_num)
   double precision, intent(in) :: coefs(N_states)
   double precision, intent(inout) :: vect(N_states, mo_tot_num)
@@ -329,7 +180,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
       vect(:, puti) += hij * coefs
     end if
   end if
-end subroutine
+end 
 
 
 
@@ -338,7 +189,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   implicit none
   
   integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
-  integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
   logical, intent(in) :: bannedOrb(mo_tot_num)
   double precision, intent(in) :: coefs(N_states)
   double precision, intent(inout) :: vect(N_states, mo_tot_num)
@@ -392,7 +243,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   call apply_particle(mask, sp, p1, det, ok,  N_int)
   call i_h_j(gen, det, N_int, hij)
   vect(:, p1) += hij * coefs
-end subroutine
+end 
 
 
 subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
@@ -400,7 +251,7 @@ subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   implicit none
   
   integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
-  integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
   logical, intent(in) :: bannedOrb(mo_tot_num)
   double precision, intent(in) :: coefs(N_states)
   double precision, intent(inout) :: vect(N_states, mo_tot_num)
@@ -418,69 +269,14 @@ subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
     call i_h_j(gen, det, N_int, hij)
     vect(:, i) += hij * coefs
   end do
-end subroutine
+end 
 
-
-subroutine spot_hasBeen(mask, sp, det, i_gen, N, banned, fullMatch)
-  use bitmasks
-  implicit none
-  
-  integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N)
-  integer, intent(in) :: i_gen, N, sp
-  logical, intent(inout) :: banned(mo_tot_num)
-  logical, intent(out) :: fullMatch
-
-
-  integer :: i, j, na, nb, list(3), nt
-  integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2)
-
-  fullMatch = .false.
-
-  do i=1,N_int
-    negMask(i,1) = not(mask(i,1))
-    negMask(i,2) = not(mask(i,2))
-  end do
-
-  genl : do i=1, N
-    nt = 0
-    
-    do j=1, N_int
-      myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1))
-      myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2))
-      nt += popcnt(myMask(j, 1)) + popcnt(myMask(j, 2))
-    end do
-    
-    if(nt > 3) cycle
-    
-    if(nt <= 2 .and. i < i_gen) then
-      fullMatch = .true.
-      return
-    end if
-    
-    call bitstring_to_list(myMask(1,sp), list(1), na, N_int)
-    
-    if(nt == 3 .and. i < i_gen) then
-      do j=1,na
-        banned(list(j)) = .true.
-      end do
-    else if(nt == 1 .and. na == 1) then
-      banned(list(1)) = .true.
-    end if
-  end do genl
-end subroutine
-
-
-
-
-! Selection double
-! ----------------
- 
-subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf)
+subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf,subset)
   use bitmasks
   use selection_types
   implicit none
   
-  integer, intent(in)            :: i_generator
+  integer, intent(in)            :: i_generator, subset
   integer(bit_kind), intent(in)  :: hole_mask(N_int,2), particle_mask(N_int,2)
   double precision, intent(in)   :: fock_diag_tmp(mo_tot_num)
   double precision, intent(in)   :: E0(N_states)
@@ -496,6 +292,14 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
   integer,allocatable               :: preinteresting(:), prefullinteresting(:), interesting(:), fullinteresting(:)
   integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
   
+  logical :: monoAdo, monoBdo;
+  integer :: maskInd
+
+  PROVIDE fragment_count
+
+  monoAdo = .true.
+  monoBdo = .true.
+  
   allocate(minilist(N_int, 2, N_det_selectors), fullminilist(N_int, 2, N_det))
   allocate(preinteresting(0:N_det_selectors), prefullinteresting(0:N_det), interesting(0:N_det_selectors), fullinteresting(0:N_det))
   
@@ -513,7 +317,24 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
   call bitstring_to_list_ab(hole    , hole_list    , N_holes    , N_int)
   call bitstring_to_list_ab(particle, particle_list, N_particles, N_int)
 
+!  ! ======
+!  ! If the subset doesn't exist, return
+!  logical :: will_compute
+!  will_compute = subset == 0
+!
+!  if (.not.will_compute) then
+!    maskInd = N_holes(1)*N_holes(2) + N_holes(2)*((N_holes(2)-1)/2) + N_holes(1)*((N_holes(1)-1)/2)
+!    will_compute = (maskInd >= subset)
+!    if (.not.will_compute) then
+!      return
+!    endif
+!  endif
+!  ! ======
+
   
+  integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
+  allocate (preinteresting_det(N_int,2,N_det))
+
   preinteresting(0) = 0
   prefullinteresting(0) = 0
   
@@ -523,17 +344,23 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
   end do
   
   do i=1,N_det
-    nt = 0
-    do j=1,N_int
+    mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
+    mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i))
+    nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) 
+    do j=2,N_int
       mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i))
       mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i))
-      nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+      nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) 
     end do
 
     if(nt <= 4) then
       if(i <= N_det_selectors) then
         preinteresting(0) += 1
         preinteresting(preinteresting(0)) = i
+        do j=1,N_int
+          preinteresting_det(j,1,preinteresting(0)) = psi_det_sorted(j,1,i)
+          preinteresting_det(j,2,preinteresting(0)) = psi_det_sorted(j,2,i)
+        enddo
       else if(nt <= 2) then
         prefullinteresting(0) += 1
         prefullinteresting(prefullinteresting(0)) = i
@@ -541,37 +368,48 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
     end if
   end do
   
-        
+
+  maskInd = -1
+  integer :: nb_count
   do s1=1,2
     do i1=N_holes(s1),1,-1   ! Generate low excitations first
       h1 = hole_list(i1,s1)
       call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int)
       
-      do i=1,N_int
-        negMask(i,1) = not(pmask(i,1))
-        negMask(i,2) = not(pmask(i,2))
-      end do
+      negMask = not(pmask)
       
       interesting(0) = 0
       fullinteresting(0) = 0
       
       do ii=1,preinteresting(0)
         i = preinteresting(ii)
-        nt = 0
-        do j=1,N_int
-          mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i))
-          mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i))
-          nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+        mobMask(1,1) = iand(negMask(1,1), preinteresting_det(1,1,ii))
+        mobMask(1,2) = iand(negMask(1,2), preinteresting_det(1,2,ii))
+        nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+        do j=2,N_int
+          mobMask(j,1) = iand(negMask(j,1), preinteresting_det(j,1,ii))
+          mobMask(j,2) = iand(negMask(j,2), preinteresting_det(j,2,ii))
+          nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
         end do
         
         if(nt <= 4) then
           interesting(0) += 1
           interesting(interesting(0)) = i
-          minilist(:,:,interesting(0)) = psi_det_sorted(:,:,i)
+          minilist(1,1,interesting(0)) = preinteresting_det(1,1,ii)
+          minilist(1,2,interesting(0)) = preinteresting_det(1,2,ii)
+          do j=2,N_int
+            minilist(j,1,interesting(0)) = preinteresting_det(j,1,ii)
+            minilist(j,2,interesting(0)) = preinteresting_det(j,2,ii)
+          enddo
           if(nt <= 2) then
             fullinteresting(0) += 1
             fullinteresting(fullinteresting(0)) = i
-            fullminilist(:,:,fullinteresting(0)) = psi_det_sorted(:,:,i)
+            fullminilist(1,1,fullinteresting(0)) = preinteresting_det(1,1,ii)
+            fullminilist(1,2,fullinteresting(0)) = preinteresting_det(1,2,ii)
+            do j=2,N_int
+              fullminilist(j,1,fullinteresting(0)) = preinteresting_det(j,1,ii)
+              fullminilist(j,2,fullinteresting(0)) = preinteresting_det(j,2,ii)
+            enddo
           end if
         end if
       end do
@@ -579,56 +417,83 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
       do ii=1,prefullinteresting(0)
         i = prefullinteresting(ii)
         nt = 0
-        do j=1,N_int
+        mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
+        mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i))
+        nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+        do j=2,N_int
           mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i))
           mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i))
-          nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+          nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
         end do
 
         if(nt <= 2) then
           fullinteresting(0) += 1
           fullinteresting(fullinteresting(0)) = i
-          fullminilist(:,:,fullinteresting(0)) = psi_det_sorted(:,:,i)
+          fullminilist(1,1,fullinteresting(0)) = psi_det_sorted(1,1,i)
+          fullminilist(1,2,fullinteresting(0)) = psi_det_sorted(1,2,i)
+          do j=2,N_int
+            fullminilist(j,1,fullinteresting(0)) = psi_det_sorted(j,1,i)
+            fullminilist(j,2,fullinteresting(0)) = psi_det_sorted(j,2,i)
+          enddo
         end if
       end do
       
+
+
       do s2=s1,2
         sp = s1
+        
         if(s1 /= s2) sp = 3
         
         ib = 1
         if(s1 == s2) ib = i1+1
+        monoAdo = .true.
         do i2=N_holes(s2),ib,-1   ! Generate low excitations first
-          
-          h2 = hole_list(i2,s2)
-          call apply_hole(pmask, s2,h2, mask, ok, N_int)
-          
           logical                        :: banned(mo_tot_num, mo_tot_num,2)
           logical                        :: bannedOrb(mo_tot_num, 2)
           
+          h2 = hole_list(i2,s2)
+          call apply_hole(pmask, s2,h2, mask, ok, N_int)
           banned = .false.
-          
-          call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
-          
-          if(fullMatch) cycle
-          
-          bannedOrb(1:mo_tot_num, 1:2) = .true.
+          do j=1,mo_tot_num
+            bannedOrb(j, 1) = .true.
+            bannedOrb(j, 2) = .true.
+          enddo
           do s3=1,2
             do i=1,N_particles(s3)
               bannedOrb(particle_list(i,s3), s3) = .false.
             enddo
           enddo
-             
-          mat = 0d0
-          call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
-          call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf)
+          if(s1 /= s2) then
+            if(monoBdo) then
+              bannedOrb(h1,s1) = .false.
+            end if
+            if(monoAdo) then
+              bannedOrb(h2,s2) = .false.
+              monoAdo = .false.
+            end if
+          end if
+
+          maskInd += 1
+          if(subset == 0 .or. mod(maskInd, fragment_count) == (subset-1)) then  
+            
+            call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
+            if(fullMatch) cycle
+          
+            mat = 0d0
+            call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
+
+            call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf)
+          end if
         enddo
+        if(s1 /= s2) monoBdo = .false.
       enddo
     enddo
   enddo
 end subroutine
 
 
+
 subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf)
   use bitmasks
   use selection_types
@@ -670,7 +535,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       if(mat(1, p1, p2) == 0d0) cycle
       call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
       
-      
       Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
       max_e_pert = 0d0
       
@@ -679,11 +543,12 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         val = mat(istate, p1, p2) + mat(istate, p1, p2) 
         tmp = dsqrt(delta_E * delta_E + val * val)
         if (delta_E < 0.d0) then
-          tmp = -tmp
+            tmp = -tmp
         endif
         e_pert = 0.5d0 * ( tmp - delta_E)
         pt2(istate) = pt2(istate) + e_pert
         max_e_pert = min(e_pert,max_e_pert)
+!        ci(istate) = e_pert / mat(istate, p1, p2)
       end do
       
       if(dabs(max_e_pert) > buf%mini) then
@@ -691,18 +556,17 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       end if
     end do
   end do
-end subroutine
+end 
 
 
 subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
   use bitmasks
   implicit none
   
-  integer, intent(in) :: interesting(0:N_sel)
-  
-  integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel)
-  integer, intent(in) :: sp, i_gen, N_sel
-  logical, intent(inout) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2)
+  integer, intent(in)            :: sp, i_gen, N_sel
+  integer, intent(in)            :: interesting(0:N_sel)
+  integer(bit_kind),intent(in)   :: mask(N_int, 2), det(N_int, 2, N_sel)
+  logical, intent(inout)         :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2)
   double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
 
   integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt
@@ -710,6 +574,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
 !   logical :: bandon
 !   
 !   bandon = .false.
+  PROVIDE psi_phasemask psi_selectors_coef_transp
   mat = 0d0
   
   do i=1,N_int
@@ -719,35 +584,32 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
 
   do i=1, N_sel ! interesting(0)
     !i = interesting(ii)
+    if (interesting(i) < 0) then
+      stop 'prefetch interesting(i)'
+    endif
+
     
-    nt = 0
-    do j=1,N_int
+    mobMask(1,1) = iand(negMask(1,1), det(1,1,i))
+    mobMask(1,2) = iand(negMask(1,2), det(1,2,i))
+    nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+
+    if(nt > 4) cycle
+
+    do j=2,N_int
       mobMask(j,1) = iand(negMask(j,1), det(j,1,i))
       mobMask(j,2) = iand(negMask(j,2), det(j,2,i))
-      nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+      nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
     end do
 
     if(nt > 4) cycle
     
-    do j=1,N_int
-      perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
-      perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
-    end do
-
-    call bitstring_to_list(perMask(1,1), h(1,1), h(0,1), N_int)
-    call bitstring_to_list(perMask(1,2), h(1,2), h(0,2), N_int)
-
-    call bitstring_to_list(mobMask(1,1), p(1,1), p(0,1), N_int)
-    call bitstring_to_list(mobMask(1,2), p(1,2), p(0,2), N_int)
-    
-    if(interesting(i) < i_gen) then
-      if(nt == 4) call past_d2(banned, p, sp)
-      if(nt == 3) call past_d1(bannedOrb, p)
-    else
-      if(interesting(i) == i_gen) then
-!         bandon = .true.
+    if (interesting(i) == i_gen) then
         if(sp == 3) then
-          banned(:,:,2) = transpose(banned(:,:,1))
+          do j=1,mo_tot_num
+            do k=1,mo_tot_num
+              banned(j,k,2) = banned(k,j,1)
+            enddo
+          enddo
         else
           do k=1,mo_tot_num
           do l=k+1,mo_tot_num
@@ -755,17 +617,35 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
           end do
           end do
         end if
-      end if
-      if(nt == 4) then
-        call get_d2(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
-      else if(nt == 3) then
-        call get_d1(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
-      else
-        call get_d0(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
-      end if
+    end if
+
+    call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
+    call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
+    
+    perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
+    perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
+    do j=2,N_int
+      perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
+      perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
+    end do
+
+    call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
+    call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
+
+    if (interesting(i) >= i_gen) then
+        if(nt == 4) then
+          call get_d2(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
+        else if(nt == 3) then
+          call get_d1(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
+        else
+          call get_d0(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
+        end if
+    else 
+        if(nt == 4) call past_d2(banned, p, sp)
+        if(nt == 3) call past_d1(bannedOrb, p)
     end if
   end do
-end subroutine
+end 
 
 
 subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
@@ -773,7 +653,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   implicit none
 
   integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
-  integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
   logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
   double precision, intent(in) :: coefs(N_states)
   double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
@@ -822,20 +702,20 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         end if
       end do
     else
-      do i = 1,2
+      h1 = h(1,1)
+      h2 = h(1,2)
       do j = 1,2
-        puti = p(i, 1)
         putj = p(j, 2)
-        
-        if(banned(puti,putj,bant)) cycle
-        p1 = p(turn2(i), 1)
         p2 = p(turn2(j), 2)
-        h1 = h(1,1)
-        h2 = h(1,2)
-        
-        hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
-        mat(:, puti, putj) += coefs * hij
-      end do
+        do i = 1,2
+          puti = p(i, 1)
+          
+          if(banned(puti,putj,bant)) cycle
+          p1 = p(turn2(i), 1)
+          
+          hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
+          mat(:, puti, putj) += coefs * hij
+        end do
       end do
     end if
 
@@ -883,7 +763,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       end if
     end if
   end if
-end subroutine
+end 
 
 
 subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
@@ -891,7 +771,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   implicit none
 
   integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
-  integer(1),intent(in) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1),intent(in) :: phasemask(2,N_int*bit_kind_size)
   logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
   integer(bit_kind) :: det(N_int, 2)
   double precision, intent(in) :: coefs(N_states)
@@ -1050,7 +930,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         mat(:, p1, p2) += coefs * hij
       end do
     end do
-end subroutine
+end 
 
 
 
@@ -1060,7 +940,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   implicit none
 
   integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
-  integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
   logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
   integer(bit_kind) :: det(N_int, 2)
   double precision, intent(in) :: coefs(N_states)
@@ -1088,8 +968,8 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
           call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
           call i_h_j(gen, det, N_int, hij)
         else
-          hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
           phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
+          hij = integral8(p1, p2, h1, h2) * phase
         end if
         mat(:, p1, p2) += coefs(:) * hij
       end do
@@ -1112,7 +992,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       end do
     end do
   end if
-end subroutine
+end 
  
 
 subroutine past_d1(bannedOrb, p)
@@ -1128,7 +1008,7 @@ subroutine past_d1(bannedOrb, p)
       bannedOrb(p(i, s), s) = .true.
     end do
   end do
-end subroutine
+end 
 
 
 subroutine past_d2(banned, p, sp)
@@ -1153,7 +1033,7 @@ subroutine past_d2(banned, p, sp)
       end do
     end do
   end if
-end subroutine
+end 
 
 
 
@@ -1161,9 +1041,9 @@ subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
   use bitmasks
   implicit none
   
+  integer, intent(in) :: i_gen, N
   integer, intent(in) :: interesting(0:N)
   integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N)
-  integer, intent(in) :: i_gen, N
   logical, intent(inout) :: banned(mo_tot_num, mo_tot_num)
   logical, intent(out) :: fullMatch
 
@@ -1194,9 +1074,37 @@ subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
       myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2))
     end do
 
-    call bitstring_to_list(myMask(1,1), list(1), na, N_int)
-    call bitstring_to_list(myMask(1,2), list(na+1), nb, N_int)
+    call bitstring_to_list_in_selection(myMask(1,1), list(1), na, N_int)
+    call bitstring_to_list_in_selection(myMask(1,2), list(na+1), nb, N_int)
     banned(list(1), list(2)) = .true.
   end do genl
-end subroutine
+end 
 
+
+subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Gives the inidices(+1) of the bits set to 1 in the bit string
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: string(Nint)
+  integer, intent(out)           :: list(Nint*bit_kind_size)
+  integer, intent(out)           :: n_elements
+  
+  integer                        :: i, ishift
+  integer(bit_kind)              :: l
+  
+  n_elements = 0
+  ishift = 2
+  do i=1,Nint
+    l = string(i)
+    do while (l /= 0_bit_kind)
+      n_elements = n_elements+1
+      list(n_elements) = ishift+popcnt(l-1_bit_kind) - popcnt(l)
+      l = iand(l,l-1_bit_kind)
+    enddo
+    ishift = ishift + bit_kind_size
+  enddo
+  
+end

plugins/Full_CI_ZMQ/selection_buffer.irp.f

@@ -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)
@@ -41,29 +41,32 @@ subroutine sort_selection_buffer(b)
   implicit none
 
   type(selection_buffer), intent(inout) :: b
-  double precision, allocatable :: vals(:), absval(:)
+  double precision, allocatable::  absval(:)
   integer, allocatable :: iorder(:)
-  integer(bit_kind), allocatable :: detmp(:,:,:)
+  double precision, pointer :: vals(:)
+  integer(bit_kind), pointer :: detmp(:,:,:)
   integer :: i, nmwen
   logical, external :: detEq
   nmwen = min(b%N, b%cur)
 
 
-  allocate(iorder(b%cur), detmp(N_int, 2, nmwen), absval(b%cur), vals(nmwen))
+  allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), absval(b%cur), vals(size(b%val)))
   absval = -dabs(b%val(:b%cur))
   do i=1,b%cur
     iorder(i) = i
   end do
   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
+  do i=nmwen+1, size(vals)
+    vals(i) = 0.d0
+  enddo
+  deallocate(b%det, b%val)
+  b%det => detmp
+  b%val => vals
   b%mini = max(b%mini,dabs(b%val(b%N)))
   b%cur = nmwen
 end subroutine

plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f

@@ -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,7 +23,7 @@ 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
   
   call provide_everything
@@ -31,6 +31,7 @@ subroutine run_wf
   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()
 
@@ -52,7 +53,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 +63,29 @@ subroutine run_wf
       ! --------
 
       print *,  'Davidson'
-      call davidson_miniserver_get()
+      call davidson_slave_tcp(i)
+      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 davidson_slave_tcp(i)
+      call run_pt2_slave(0,i,energy,lstop)
       !$OMP END PARALLEL
-      print *,  'Davidson done'
+      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
 

plugins/Full_CI_ZMQ/selection_slave.irp.f

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

plugins/Full_CI_ZMQ/selection_types.f90

@@ -1,9 +1,9 @@
 module selection_types
   type selection_buffer
     integer :: N, cur
-    integer(8), allocatable :: det(:,:,:)
-    double precision, allocatable :: val(:)
-    double precision :: mini
+    integer(8)      , pointer :: det(:,:,:)
+    double precision, pointer :: val(:)
+    double precision          :: mini
   endtype
 end module
 

plugins/Full_CI_ZMQ/target_pt2_ratio_zmq.irp.f

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

plugins/Full_CI_ZMQ/target_pt2_zmq.irp.f

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

plugins/Full_CI_ZMQ/zmq_selection.irp.f

@@ -0,0 +1,126 @@
+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)
+  integer, parameter             :: maxtasks=10000
+  
+  
+  PROVIDE fragment_count
+
+  N = max(N_in,1)
+  if (.True.) then
+    PROVIDE pt2_e0_denominator
+    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 create_selection_buffer(N, N*2, b)
+  endif
+
+  character*(20*maxtasks) :: task
+  task = ' ' 
+
+  integer :: k
+  k=0
+  do i= 1, N_det_generators
+    k = k+1
+    write(task(20*(k-1)+1:20*k),'(I9,1X,I9,''|'')') i, N
+    if (k>=maxtasks) then
+       k=0
+       call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+    endif
+  end do
+  if (k > 0) then
+    call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+  endif
+  call zmq_set_running(zmq_to_qp_run_socket)
+
+  !$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_generators))
+  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
+

plugins/Generators_CAS/.gitignore

@@ -1,23 +0,0 @@
-# Automatically created by /home/razoa/quantum_package/scripts/module/module_handler.py 
-IRPF90_temp
-IRPF90_man
-irpf90_entities
-tags
-irpf90.make
-Makefile
-Makefile.depend
-build.ninja
-.ninja_log
-.ninja_deps
-ezfio_interface.irp.f
-Ezfio_files
-Determinants
-Integrals_Monoelec
-MO_Basis
-Utils
-Pseudo
-Bitmask
-AO_Basis
-Electrons
-Nuclei
-Integrals_Bielec

plugins/Generators_full/.gitignore

@@ -1,25 +0,0 @@
-# Automatically created by /home/razoa/quantum_package/scripts/module/module_handler.py 
-IRPF90_temp
-IRPF90_man
-irpf90_entities
-tags
-irpf90.make
-Makefile
-Makefile.depend
-build.ninja
-.ninja_log
-.ninja_deps
-ezfio_interface.irp.f
-Ezfio_files
-Determinants
-Integrals_Monoelec
-MO_Basis
-Utils
-Pseudo
-Bitmask
-AO_Basis
-Electrons
-MOGuess
-Nuclei
-Hartree_Fock
-Integrals_Bielec

plugins/Generators_restart/.gitignore

@@ -1,13 +0,0 @@
-#
-# Do not modify this file. Add your ignored files to the gitignore
-# (without the dot at the beginning) file.
-#
-IRPF90_temp
-IRPF90_man
-irpf90.make
-tags
-Makefile.depend
-irpf90_entities
-build.ninja
-.ninja_log
-.ninja_deps

plugins/Hartree_Fock/.gitignore

@@ -1,25 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Electrons
-Ezfio_files
-Huckel_guess
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-SCF
-Utils
-ZMQ
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-tags

plugins/Hartree_Fock/localize_mos.irp.f

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

plugins/MP2/.gitignore

@@ -1,31 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-SingleRefMethod
-Utils
-ZMQ
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-mp2
-mp2_wf
-tags

plugins/MRCC_Utils/.gitignore

@@ -1,33 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Davidson
-Determinants
-Electrons
-Ezfio_files
-Generators_full
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Psiref_CAS
-Psiref_Utils
-Selectors_full
-Utils
-ZMQ
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-mrcc_dummy
-tags

plugins/MRCC_Utils/H_apply.irp.f

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

plugins/MRCC_Utils/amplitudes.irp.f

@@ -23,33 +23,39 @@
  allocate(pathTo(N_det_non_ref))
 
  pathTo(:) = 0
- is_active_exc(:) = .false.
+ is_active_exc(:) = .True. 
  n_exc_active = 0
 
- do hh = 1, hh_shortcut(0)
-   do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
-     do II = 1, N_det_ref
+! do hh = 1, hh_shortcut(0)
+!   do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
+!     do II = 1, N_det_ref
+!
+!       call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
+!       if(.not. ok) cycle
+!
+!       call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
+!       if(.not. ok) cycle
+!
+!       ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
+!       if(ind == -1) cycle
+!
+!       logical, external :: is_a_two_holes_two_particles
+!       if (is_a_two_holes_two_particles(myDet)) then
+!         is_active_exc(pp) = .False.
+!       endif
 
-       call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
-       if(.not. ok) cycle
+!       ind = psi_non_ref_sorted_idx(ind)
+!       if(pathTo(ind) == 0) then
+!         pathTo(ind) = pp
+!       else
+!         is_active_exc(pp) = .true.
+!         is_active_exc(pathTo(ind)) = .true.
+!       end if
 
-       call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
-       if(.not. ok) cycle
+!     end do
+!   end do
+! end do
 
-       ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
-       if(ind == -1) cycle
-
-       ind = psi_non_ref_sorted_idx(ind)
-       if(pathTo(ind) == 0) then
-         pathTo(ind) = pp
-       else
-         is_active_exc(pp) = .true.
-         is_active_exc(pathTo(ind)) = .true.
-       end if
-     end do
-   end do
- end do
-!is_active_exc=.true.
  do hh = 1, hh_shortcut(0)
   do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
     if(is_active_exc(pp)) then
@@ -66,6 +72,32 @@
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ logical, has_a_unique_parent, (N_det_non_ref) ]
+ implicit none
+ BEGIN_DOC
+ ! True if the determinant in the non-reference has a unique parent
+ END_DOC
+ integer :: i,j,n
+ integer :: degree
+ do j=1,N_det_non_ref
+  has_a_unique_parent(j) = .True. 
+  n=0
+  do i=1,N_det_ref
+   call get_excitation_degree(psi_ref(1,1,i), psi_non_ref(1,1,j), degree, N_int)
+   if (degree < 2) then
+     n = n+1
+     if (n > 1) then
+       has_a_unique_parent(j) = .False.
+       exit
+     endif
+   endif
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+
 BEGIN_PROVIDER [ integer, n_exc_active_sze ]
  implicit none
  BEGIN_DOC
@@ -89,14 +121,13 @@ END_PROVIDER
  double precision               :: phase
  logical                        :: ok
  integer, external              :: searchDet
-
- PROVIDE psi_non_ref_sorted_idx psi_ref_coef 
+ 
 
  !$OMP PARALLEL default(none) shared(psi_non_ref, hh_exists, pp_exists, N_int,&
      !$OMP   active_excitation_to_determinants_val, active_excitation_to_determinants_idx)&
      !$OMP shared(hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, &
      !$OMP   psi_non_ref_sorted_idx, psi_ref, N_det_ref, N_states)&
-     !$OMP shared(is_active_exc, active_hh_idx, active_pp_idx, n_exc_active)&
+     !$OMP shared(active_hh_idx, active_pp_idx, n_exc_active)&
      !$OMP private(lref, pp, II, ok, myMask, myDet, ind, phase, wk, ppp, hh, s)
  allocate(lref(N_det_non_ref))
  !$OMP DO schedule(dynamic)
@@ -127,7 +158,6 @@ END_PROVIDER
        wk += 1
        do s=1,N_states
         active_excitation_to_determinants_val(s,wk, ppp) = psi_ref_coef(lref(i), s)
-
        enddo
        active_excitation_to_determinants_idx(wk, ppp) = i
      else if(lref(i) < 0) then
@@ -160,7 +190,7 @@ END_PROVIDER
  double precision, allocatable  :: t(:), A_val_mwen(:,:), As2_val_mwen(:,:)
  integer, allocatable           :: A_ind_mwen(:)
  double precision               :: sij
- PROVIDE psi_non_ref active_excitation_to_determinants_val
+ PROVIDE psi_non_ref 
 
  mrcc_AtA_ind(:) = 0
  mrcc_AtA_val(:,:) = 0.d0
@@ -168,6 +198,7 @@ END_PROVIDER
  mrcc_N_col(:) = 0
  AtA_size = 0
 
+
   !$OMP PARALLEL default(none) shared(k, active_excitation_to_determinants_idx,&
       !$OMP   active_excitation_to_determinants_val, hh_nex)         &
       !$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen,&

plugins/MRCC_Utils/davidson.irp.f

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

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -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) ]
@@ -291,11 +351,11 @@ logical function is_generable(det1, det2, Nint)
   integer, intent(in) :: Nint
   integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
   integer :: degree, f, exc(0:2, 2, 2), t
-  integer*2 :: h1, h2, p1, p2, s1, s2
+  integer :: h1, h2, p1, p2, s1, s2
   integer, external :: searchExc
   logical, external :: excEq
   double precision :: phase
-  integer*2 :: tmp_array(4)
+  integer :: tmp_array(4)
   
   is_generable = .false.
   call get_excitation(det1, det2, exc, degree, phase, Nint)
@@ -306,7 +366,7 @@ logical function is_generable(det1, det2, Nint)
   end if
   if(degree > 2) stop "?22??"
   
-  call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
+  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
   
   if(degree == 1) then
     h2 = h1
@@ -394,7 +454,7 @@ integer function searchExc(excs, exc, n)
   use bitmasks
   
   integer, intent(in) :: n
-  integer*2,intent(in) :: excs(4,n), exc(4)
+  integer,intent(in) :: excs(4,n), exc(4)
   integer :: l, h, c
   integer, external :: excCmp
   logical, external :: excEq
@@ -459,8 +519,8 @@ subroutine sort_exc(key, N_key)
   
 
   integer, intent(in)                   :: N_key
-  integer*2,intent(inout)       :: key(4,N_key)
-  integer*2                     :: tmp(4)
+  integer,intent(inout)       :: key(4,N_key)
+  integer                     :: tmp(4)
   integer                               :: i,ni
   
   
@@ -482,7 +542,7 @@ end subroutine
 
 logical function exc_inf(exc1, exc2)
   implicit none
-  integer*2,intent(in) :: exc1(4), exc2(4)
+  integer,intent(in) :: exc1(4), exc2(4)
   integer :: i
   exc_inf = .false.
   do i=1,4
@@ -504,9 +564,9 @@ subroutine tamise_exc(key, no, n, N_key)
 ! Uncodumented : TODO
   END_DOC
   integer,intent(in)            :: no, n, N_key
-  integer*2,intent(inout)       :: key(4, N_key)
+  integer,intent(inout)       :: key(4, N_key)
   integer                       :: k,j
-  integer*2                     :: tmp(4)
+  integer                     :: tmp(4)
   logical                       :: exc_inf
   integer                       :: ni
   
@@ -535,8 +595,9 @@ end subroutine
 
 subroutine dec_exc(exc, h1, h2, p1, p2)
   implicit none
-  integer :: exc(0:2,2,2), s1, s2, degree
-  integer*2, intent(out) :: h1, h2, p1, p2
+  integer, intent(in)            :: exc(0:2,2,2)
+  integer, intent(out)           :: h1, h2, p1, p2
+  integer                        :: degree, s1, s2
   
   degree = exc(0,1,1) + exc(0,1,2)
   
@@ -547,7 +608,7 @@ subroutine dec_exc(exc, h1, h2, p1, p2)
     
   if(degree == 0) return
   
-  call decode_exc_int2(exc, degree, h1, p1, h2, p2, s1, s2)
+  call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2)
   
   h1 += mo_tot_num * (s1-1)
   p1 += mo_tot_num * (s1-1)
@@ -579,7 +640,7 @@ end subroutine
 &BEGIN_PROVIDER [ integer, N_ex_exists ]
   implicit none
   integer :: exc(0:2, 2, 2), degree, n, on, s, l, i
-  integer*2 :: h1, h2, p1, p2
+  integer :: h1, h2, p1, p2
   double precision :: phase
   logical,allocatable :: hh(:,:) , pp(:,:)
   
@@ -632,12 +693,12 @@ END_PROVIDER
   double precision               :: phase
   
   
-  double precision, allocatable :: rho_mrcc_init(:)
+  double precision, allocatable :: rho_mrcc_inact(:)
   integer                        :: a_coll, at_roww
   
   print *, "TI", hh_nex, N_det_non_ref
 
-  allocate(rho_mrcc_init(N_det_non_ref))
+  allocate(rho_mrcc_inact(N_det_non_ref))
   allocate(x_new(hh_nex))
   allocate(x(hh_nex), AtB(hh_nex))
 
@@ -649,7 +710,7 @@ END_PROVIDER
         !$OMP private(at_row, a_col, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
         !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtB, mrcc_AtA_val, mrcc_AtA_ind, s, n_exc_active, active_pp_idx)
     
-    !$OMP DO schedule(dynamic, 100)
+    !$OMP DO schedule(static, 100)
     do at_roww = 1, n_exc_active ! hh_nex
       at_row = active_pp_idx(at_roww)
       do i=1,active_excitation_to_determinants_idx(0,at_roww)
@@ -668,7 +729,7 @@ END_PROVIDER
       X(a_col) = AtB(a_col)
     end do
     
-    rho_mrcc_init = 0d0
+    rho_mrcc_inact(:) = 0d0
     
     allocate(lref(N_det_ref))
     do hh = 1, hh_shortcut(0)
@@ -692,29 +753,23 @@ END_PROVIDER
         X(pp) =  AtB(pp) 
         do II=1,N_det_ref
           if(lref(II) > 0) then
-            rho_mrcc_init(lref(II)) = psi_ref_coef(II,s) * X(pp)
+            rho_mrcc_inact(lref(II)) = psi_ref_coef(II,s) * X(pp)
           else if(lref(II) < 0) then
-            rho_mrcc_init(-lref(II)) = -psi_ref_coef(II,s) * X(pp)
+            rho_mrcc_inact(-lref(II)) = -psi_ref_coef(II,s) * X(pp)
           end if
         end do
       end do
     end do
     deallocate(lref)
 
-    do i=1,N_det_non_ref
-      rho_mrcc(i,s) = rho_mrcc_init(i) 
-    enddo
-
     x_new = x
     
     double precision               :: factor, resold
     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
@@ -725,102 +780,108 @@ 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 ) then
+        exit
+      endif
+      if ( (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
 
-  do s=1,N_states
+    do i=1,N_det_non_ref
+      rho_mrcc(i,s) = 0.d0
+    enddo
 
     do a_coll=1,n_exc_active
       a_col = active_pp_idx(a_coll)
       do j=1,N_det_non_ref
         i = active_excitation_to_determinants_idx(j,a_coll)
         if (i==0) exit
+        if (rho_mrcc_inact(i) /= 0.d0) then
+          call debug_det(psi_non_ref(1,1,i),N_int)
+          stop
+        endif
         rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * dIj_unique(a_col,s)
       enddo
     end do
 
-    norm = 0.d0
-    do i=1,N_det_non_ref
-      norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
-    enddo
-    ! Norm now contains the norm of A.X
-    
+    double precision :: norm2_ref, norm2_inact, a, b, c, Delta
+    ! Psi = Psi_ref + Psi_inactive + f*Psi_active
+    ! Find f to normalize Psi
+
+    norm2_ref = 0.d0
     do i=1,N_det_ref
-      norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
+      norm2_ref = norm2_ref + psi_ref_coef(i,s)*psi_ref_coef(i,s)
     enddo
-    ! Norm now contains the norm of Psi + A.X
-    
+
+    a = 0.d0
+    do i=1,N_det_non_ref
+      a = a + rho_mrcc(i,s)*rho_mrcc(i,s)
+    enddo
+
+    norm = a + norm2_ref
     print *, "norm : ", sqrt(norm)
-   enddo
+
+    norm = sqrt((1.d0-norm2_ref)/a)
+
+    ! Renormalize Psi+A.X
+    do i=1,N_det_non_ref
+      rho_mrcc(i,s) = rho_mrcc(i,s) * norm 
+    enddo
+
+!norm = norm2_ref
+!do i=1,N_det_non_ref
+!  norm = norm + rho_mrcc(i,s)**2
+!enddo
+!print *,  'check', norm
+!stop
+
      
         
-   do s=1,N_states
      norm = 0.d0
-     double precision               :: f
+     double precision               :: f, g, gmax
+     gmax = 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 = 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)
+       norm = norm + (rho_mrcc(i,s)*f)**2
        rho_mrcc(i,s) = f
      enddo
-     ! norm now contains the norm of |T.Psi_0>
-     ! rho_mrcc now contains the f factors
+     ! rho_mrcc now contains the mu_i factors
 
-     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)
-     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)
-     if (dsqrt(norm) > 1.d0) then
+     if (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
-     ! rho_mrcc now contains the product of the scaling factors and the
-     ! normalization constant
-    
   end do
 
 END_PROVIDER
@@ -845,6 +906,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
@@ -872,11 +980,11 @@ double precision function get_dij(det1, det2, s, Nint)
   integer, intent(in) :: s, Nint
   integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
   integer :: degree, f, exc(0:2, 2, 2), t
-  integer*2 :: h1, h2, p1, p2, s1, s2
+  integer :: h1, h2, p1, p2, s1, s2
   integer, external :: searchExc
   logical, external :: excEq
   double precision :: phase
-  integer*2 :: tmp_array(4)
+  integer :: tmp_array(4)
   
   get_dij = 0d0
   call get_excitation(det1, det2, exc, degree, phase, Nint)
@@ -885,7 +993,7 @@ double precision function get_dij(det1, det2, s, Nint)
     stop "get_dij"
   end if
   
-  call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
+  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
   
   if(degree == 1) then
     h2 = h1
@@ -918,8 +1026,8 @@ double precision function get_dij(det1, det2, s, Nint)
 end function
 
 
- BEGIN_PROVIDER [ integer*2, hh_exists, (4, N_hh_exists) ]
-&BEGIN_PROVIDER [ integer*2, pp_exists, (4, N_pp_exists) ]
+ BEGIN_PROVIDER [ integer, hh_exists, (4, N_hh_exists) ]
+&BEGIN_PROVIDER [ integer, pp_exists, (4, N_pp_exists) ]
 &BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_hh_exists + 1) ]
 &BEGIN_PROVIDER [ integer, hh_nex ]
   implicit none
@@ -934,9 +1042,9 @@ end function
   ! hh_nex : Total number of excitation operators
   !
   END_DOC
-  integer*2,allocatable :: num(:,:)
+  integer,allocatable :: num(:,:)
   integer :: exc(0:2, 2, 2), degree, n, on, s, l, i
-  integer*2 :: h1, h2, p1, p2
+  integer :: h1, h2, p1, p2
   double precision :: phase
   logical, external :: excEq
   
@@ -962,24 +1070,40 @@ end function
   
   hh_shortcut(0) = 1
   hh_shortcut(1) = 1
-  hh_exists(:,1) = (/1_2, num(1,1), 1_2, num(2,1)/)
-  pp_exists(:,1) = (/1_2, num(3,1), 1_2, num(4,1)/)
+  hh_exists(:,1) = (/1, num(1,1), 1, num(2,1)/)
+  pp_exists(:,1) = (/1, num(3,1), 1, num(4,1)/)
   s = 1
   do i=2,n
     if(.not. excEq(num(1,i), num(1,s))) then
       s += 1
       num(:, s) = num(:, i)
-      pp_exists(:,s) = (/1_2, num(3,s), 1_2, num(4,s)/)
+      pp_exists(:,s) = (/1, num(3,s), 1, num(4,s)/)
       if(hh_exists(2, hh_shortcut(0)) /= num(1,s) .or. &
             hh_exists(4, hh_shortcut(0)) /= num(2,s)) then
         hh_shortcut(0) += 1
         hh_shortcut(hh_shortcut(0)) = s
-        hh_exists(:,hh_shortcut(0)) = (/1_2, num(1,s), 1_2, num(2,s)/)
+        hh_exists(:,hh_shortcut(0)) = (/1, num(1,s), 1, num(2,s)/)
       end if
     end if
   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
@@ -990,6 +1114,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
@@ -1005,7 +1130,7 @@ END_PROVIDER
 
 logical function excEq(exc1, exc2)
   implicit none
-  integer*2, intent(in) :: exc1(4), exc2(4)
+  integer, intent(in) :: exc1(4), exc2(4)
   integer :: i
   excEq = .false.
   do i=1, 4
@@ -1017,7 +1142,7 @@ end function
 
 integer function excCmp(exc1, exc2)
   implicit none
-  integer*2, intent(in) :: exc1(4), exc2(4)
+  integer, intent(in) :: exc1(4), exc2(4)
   integer :: i
   excCmp = 0
   do i=1, 4
@@ -1036,8 +1161,8 @@ subroutine apply_hole_local(det, exc, res, ok, Nint)
   use bitmasks
   implicit none
   integer, intent(in) :: Nint
-  integer*2, intent(in) :: exc(4)
-  integer*2 :: s1, s2, h1, h2
+  integer, intent(in) :: exc(4)
+  integer :: s1, s2, h1, h2
   integer(bit_kind),intent(in) :: det(Nint, 2)
   integer(bit_kind),intent(out) :: res(Nint, 2)
   logical, intent(out) :: ok
@@ -1073,8 +1198,8 @@ subroutine apply_particle_local(det, exc, res, ok, Nint)
   use bitmasks
   implicit none
   integer, intent(in) :: Nint
-  integer*2, intent(in) :: exc(4)
-  integer*2 :: s1, s2, p1, p2
+  integer, intent(in) :: exc(4)
+  integer :: s1, s2, p1, p2
   integer(bit_kind),intent(in) :: det(Nint, 2)
   integer(bit_kind),intent(out) :: res(Nint, 2)
   logical, intent(out) :: ok

plugins/MRPT_Utils/energies_cas.irp.f

@@ -898,7 +898,7 @@ END_PROVIDER
       enddo
       print*, '***'
       do i = 1, N_det+1
-       write(*,'(100(F16.10,X))')H_matrix(i,:)
+       write(*,'(100(F16.10,1X))')H_matrix(i,:)
       enddo
       call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1)
       corr_e_from_1h1p(state_target) += eigenvalues(1) - energy_cas_dyall(state_target)   
@@ -919,15 +919,15 @@ END_PROVIDER
        norm += psi_in_out_coef(i,state_target) * psi_in_out_coef(i,state_target)
       enddo
       print*, 'Coef '
-      write(*,'(100(X,F16.10))')psi_coef(1:N_det,state_target)
-      write(*,'(100(X,F16.10))')psi_in_out_coef(:,state_target)
+      write(*,'(100(1X,F16.10))')psi_coef(1:N_det,state_target)
+      write(*,'(100(1X,F16.10))')psi_in_out_coef(:,state_target)
       double precision :: coef_tmp(N_det)
       do i = 1, N_det
        coef_tmp(i) = psi_coef(i,1) * interact_psi0(i) / delta_e_alpha_beta(i,ispin)
       enddo
-      write(*,'(100(X,F16.10))')coef_tmp(:)
+      write(*,'(100(1X,F16.10))')coef_tmp(:)
       print*, 'naked interactions'
-      write(*,'(100(X,F16.10))')interact_psi0(:)
+      write(*,'(100(1X,F16.10))')interact_psi0(:)
       print*, ''
       
       print*, 'norm ',norm
@@ -953,10 +953,10 @@ END_PROVIDER
       enddo
     enddo
       print*, '***'
-      write(*,'(100(X,F16.10))')
-      write(*,'(100(X,F16.10))')delta_e_alpha_beta(:,2)
- !    write(*,'(100(X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,1,:)
- !    write(*,'(100(X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,2,:)
+      write(*,'(100(1X,F16.10))')
+      write(*,'(100(1X,F16.10))')delta_e_alpha_beta(:,2)
+ !    write(*,'(100(1X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,1,:)
+ !    write(*,'(100(1X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,2,:)
     print*, '---------------------------------------------------------------------------'
    enddo
   enddo
@@ -1089,11 +1089,11 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
        print*, 'e corr perturb EN',accu(state_target)
        print*, ''
        print*, 'coef diagonalized'
-       write(*,'(100(F16.10,X))')psi_in_out_coef(:,state_target)
+       write(*,'(100(F16.10,1X))')psi_in_out_coef(:,state_target)
        print*, 'coef_perturb'
-       write(*,'(100(F16.10,X))')coef_perturb(:)
+       write(*,'(100(F16.10,1X))')coef_perturb(:)
        print*, 'coef_perturb EN'
-       write(*,'(100(F16.10,X))')coef_perturb_bis(:)
+       write(*,'(100(F16.10,1X))')coef_perturb_bis(:)
       endif
       integer :: k
       do k = 1, N_det

plugins/MRPT_Utils/excitations_cas.irp.f

@@ -22,7 +22,7 @@ subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, &
 
   integer :: elec_num_tab_local(2)
   integer :: i,j,accu_elec,k
-  integer :: det_tmp(N_int), det_tmp_bis(N_int)
+  integer(bit_kind) :: det_tmp(N_int), det_tmp_bis(N_int)
   double precision :: phase
   double precision :: norm_factor
 

plugins/MRPT_Utils/ezfio_interface.irp.f

@@ -1,23 +0,0 @@
-! DO NOT MODIFY BY HAND
-! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
-! from file /home/scemama/quantum_package/src/MRPT_Utils/EZFIO.cfg
-
-
-BEGIN_PROVIDER [ logical, do_third_order_1h1p  ]
-  implicit none
-  BEGIN_DOC
-! If true, compute the third order contribution for the 1h1p
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  
-  call ezfio_has_mrpt_utils_do_third_order_1h1p(has)
-  if (has) then
-    call ezfio_get_mrpt_utils_do_third_order_1h1p(do_third_order_1h1p)
-  else
-    print *, 'mrpt_utils/do_third_order_1h1p not found in EZFIO file'
-    stop 1
-  endif
-
-END_PROVIDER

plugins/MRPT_Utils/new_way_second_order_coef.irp.f

@@ -210,7 +210,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
 
                  ! < det_tmp     | H | det_tmp_bis >  = F_{aorb,borb}
                  hab = (fock_operator_local(aorb,borb,kspin) ) * phase  
-            if(isnan(hab))then
+            if(hab /= hab)then ! check NaN
              print*, '1'
              stop
             endif
@@ -255,7 +255,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
                  call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
 !                ! < det_tmp | H | det_tmp_bis >  = F_{aorb,borb}
                  hab = fock_operator_local(aorb,borb,kspin) * phase  
-            if(isnan(hab))then
+            if(hab /= hab)then ! check NaN
              print*, '2'
              stop
             endif

plugins/Molden/.gitignore

@@ -1,18 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Utils
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-print_mo
-tags

plugins/Perturbation/.gitignore

@@ -1,26 +0,0 @@
-# Automatically created by /home/razoa/quantum_package/scripts/module/module_handler.py 
-IRPF90_temp
-IRPF90_man
-irpf90_entities
-tags
-irpf90.make
-Makefile
-Makefile.depend
-build.ninja
-.ninja_log
-.ninja_deps
-ezfio_interface.irp.f
-Ezfio_files
-Determinants
-Integrals_Monoelec
-MO_Basis
-Utils
-Pseudo
-Properties
-Bitmask
-AO_Basis
-Electrons
-MOGuess
-Nuclei
-Hartree_Fock
-Integrals_Bielec

plugins/Properties/.gitignore

@@ -1,25 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Utils
-ZMQ
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-print_hcc
-print_mulliken
-tags

plugins/Properties/delta_rho.irp.f

@@ -6,7 +6,7 @@
  z_min = 0.d0
  z_max = 10.d0
  delta_z = 0.005d0
- N_z_pts = (z_max - z_min)/delta_z
+ N_z_pts = int( (z_max - z_min)/delta_z )
  print*,'N_z_pts = ',N_z_pts
 
 END_PROVIDER

plugins/Properties/hyperfine_constants.irp.f

@@ -151,7 +151,7 @@ subroutine print_hcc
  integer :: i,j
  print*,'Z               AU           GAUSS              MHZ             cm^-1'
  do i = 1, nucl_num
-  write(*,'(I2,X,F4.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
+  write(*,'(I2,1X,F4.1,1X,4(F16.6,1X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
  enddo
 
 end

plugins/Properties/mulliken.irp.f

@@ -126,7 +126,7 @@ subroutine print_mulliken_sd
  accu = 0.d0
  do i = 1, ao_num
   accu += spin_gross_orbital_product(i)
-  write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
+  write(*,'(1X,I3,1X,A4,1X,I2,1X,A4,1X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
  enddo
  print*,'sum = ',accu
  accu = 0.d0
@@ -142,7 +142,7 @@ subroutine print_mulliken_sd
   accu = 0.d0
   do i = 0,  ao_l_max
    accu += spin_population_angular_momentum_per_atom(i,j)
-   write(*,'(XX,I3,XX,A4,X,A4,X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j)
+   write(*,'(1X,I3,1X,A4,1X,A4,1X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j)
    print*,'sum = ',accu
   enddo
  enddo

plugins/Psiref_CAS/.gitignore

@@ -1,26 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Davidson
-Determinants
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Psiref_Utils
-Utils
-ZMQ
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-overwrite_with_cas
-tags

plugins/Psiref_CAS/psi_ref.irp.f

@@ -67,3 +67,37 @@ END_PROVIDER
 
 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 

plugins/Psiref_Utils/.gitignore

@@ -1,29 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Generators_full
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-Utils
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-mrcc_general
-tags

plugins/Psiref_Utils/psi_ref_utils.irp.f

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

plugins/Psiref_threshold/.gitignore

@@ -1,29 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Generators_full
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Selectors_full
-Utils
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-mrcc_general
-tags

plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES

@@ -1 +0,0 @@
-Psiref_Utils

plugins/Psiref_threshold/README.rst

@@ -1,24 +0,0 @@
-=======================
-Psiref_threshold Module
-=======================
-
-
-Reference wave function is defined as all determinants with coefficients
-such that | c_i/c_max | > threshold.
-
-Documentation
-=============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-Needed Modules
-==============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-.. image:: tree_dependency.png
-
-* `Psiref_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils>`_
-

plugins/Psiref_threshold/psi_ref.irp.f

@@ -1,41 +0,0 @@
-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.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))
-    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
-

plugins/QmcChem/.gitignore

@@ -1,25 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Pseudo
-Utils
-ZMQ
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-save_for_qmcchem
-tags
-target_pt2_qmc

plugins/Selectors_full/.gitignore

@@ -1,25 +0,0 @@
-# Automatically created by /home/razoa/quantum_package/scripts/module/module_handler.py 
-IRPF90_temp
-IRPF90_man
-irpf90_entities
-tags
-irpf90.make
-Makefile
-Makefile.depend
-build.ninja
-.ninja_log
-.ninja_deps
-ezfio_interface.irp.f
-Ezfio_files
-Determinants
-Integrals_Monoelec
-MO_Basis
-Utils
-Pseudo
-Bitmask
-AO_Basis
-Electrons
-MOGuess
-Nuclei
-Hartree_Fock
-Integrals_Bielec

plugins/Selectors_full/zmq.irp.f

@@ -90,13 +90,13 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id, energy, size_energy)
   psi_det_size = psi_det_size_read
   TOUCH psi_det_size N_det N_states
 
-  rc = f77_zmq_recv(zmq_to_qp_run_socket,psi_det,N_int*2*N_det*bit_kind,ZMQ_SNDMORE)
+  rc = f77_zmq_recv(zmq_to_qp_run_socket,psi_det,N_int*2*N_det*bit_kind,0)
   if (rc /= N_int*2*N_det*bit_kind) then
     print *, 'f77_zmq_recv(zmq_to_qp_run_socket,psi_det,N_int*2*N_det*bit_kind,ZMQ_SNDMORE)'
     stop 'error'
   endif
 
-  rc = f77_zmq_recv(zmq_to_qp_run_socket,psi_coef,psi_det_size*N_states*8,ZMQ_SNDMORE)
+  rc = f77_zmq_recv(zmq_to_qp_run_socket,psi_coef,psi_det_size*N_states*8,0)
   if (rc /= psi_det_size*N_states*8) then
     print *, '77_zmq_recv(zmq_to_qp_run_socket,psi_coef,psi_det_size*N_states*8,ZMQ_SNDMORE)'
     stop 'error'

plugins/Selectors_no_sorted/.gitignore

@@ -1,13 +0,0 @@
-#
-# Do not modify this file. Add your ignored files to the gitignore
-# (without the dot at the beginning) file.
-#
-IRPF90_temp
-IRPF90_man
-irpf90.make
-tags
-Makefile.depend
-irpf90_entities
-build.ninja
-.ninja_log
-.ninja_deps

plugins/SingleRefMethod/.gitignore

@@ -1,19 +0,0 @@
-# Automatically created by /home/razoa/quantum_package/scripts/module/module_handler.py 
-IRPF90_temp
-IRPF90_man
-irpf90_entities
-tags
-irpf90.make
-Makefile
-Makefile.depend
-build.ninja
-.ninja_log
-.ninja_deps
-ezfio_interface.irp.f
-Ezfio_files
-MO_Basis
-Utils
-Bitmask
-AO_Basis
-Electrons
-Nuclei

plugins/analyze_wf/README.rst

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

plugins/analyze_wf/analyze_wf.irp.f

@@ -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,istate)
+    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.97d0) 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.03d0) 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

plugins/analyze_wf/occupation.irp.f

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

plugins/loc_cele/.gitignore

@@ -1,18 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Electrons
-Ezfio_files
-IRPF90_man
-IRPF90_temp
-MO_Basis
-Makefile
-Makefile.depend
-Nuclei
-Utils
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-loc_cele
-tags

plugins/mrcc_selected/dressing_slave.irp.f

@@ -1,601 +0,0 @@
-subroutine mrsc2_dressing_slave_tcp(i)
-  implicit none
-  integer, intent(in)            :: i
-  BEGIN_DOC
-! Task for parallel MR-SC2
-  END_DOC
-  call mrsc2_dressing_slave(0,i)
-end
-
-
-subroutine mrsc2_dressing_slave_inproc(i)
-  implicit none
-  integer, intent(in)            :: i
-  BEGIN_DOC
-! Task for parallel MR-SC2
-  END_DOC
-  call mrsc2_dressing_slave(1,i)
-end
-
-subroutine mrsc2_dressing_slave(thread,iproc)
-  use f77_zmq
-
-  implicit none
-  BEGIN_DOC
-! Task for parallel MR-SC2
-  END_DOC
-  integer,  intent(in)            :: thread, iproc
-!   integer                        :: j,l
-  integer                        :: rc 
-
-  integer                        :: worker_id, task_id 
-  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 
-
-  double precision, allocatable  :: delta(:,:,:), delta_s2(:,:,:)
-  
-
-
-  integer                         :: i_state, i, i_I, J, k, k2, k1, kk, ll, degree, degree2, m, l, deg, ni, m2
-  integer                         :: n(2)
-  integer                         :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, kn
-  logical                         :: ok
-  double precision                :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al
-  double precision                :: diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv(N_states), cj_inv(N_states)
-  double precision                :: contrib, contrib_s2, wall, iwall
-  double precision, allocatable   :: dleat(:,:,:), dleat_s2(:,:,:)
-  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
-  integer(bit_kind)               :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt
-  integer, external               :: get_index_in_psi_det_sorted_bit, searchDet, detCmp
-  logical, external               :: is_in_wavefunction, isInCassd, detEq
-  integer,allocatable :: komon(:)
-  logical :: komoned
-  !double precision, external :: get_dij
-     
-  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)
-
-  allocate (dleat(N_states, N_det_non_ref, 2), delta(N_states,0:N_det_non_ref, 2))
-  allocate (dleat_s2(N_states, N_det_non_ref, 2), delta_s2(N_states,0:N_det_non_ref, 2))
-  allocate(komon(0:N_det_non_ref))
-
-  do 
-    call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task)
-    if (task_id == 0) exit
-    read (task,*) i_I, J, k1, k2
-    do i_state=1, N_states
-      ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state)
-      cj_inv(i_state) = 1.d0 / psi_ref_coef(J,i_state)
-    end do
-    n = 0
-    delta(:,0,:) = 0d0
-    delta(:,:nlink(J),1) = 0d0
-    delta(:,:nlink(i_I),2) = 0d0
-    delta_s2(:,0,:) = 0d0
-    delta_s2(:,:nlink(J),1) = 0d0
-    delta_s2(:,:nlink(i_I),2) = 0d0
-    komon(0) = 0
-    komoned = .false.
-    
-    
-    
-    
-    do kk = k1, k2
-      k = det_cepa0_idx(linked(kk, i_I))
-      blok = blokMwen(kk, i_I)
-      
-      call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int)
-    
-      if(J /= i_I) then
-        call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp2,ok,N_int)
-        if(.not. ok) cycle
-          
-        l = searchDet(det_cepa0(1,1,cepa0_shortcut(blok)), det_tmp2, cepa0_shortcut(blok+1)-cepa0_shortcut(blok), N_int)
-        if(l == -1) cycle
-        ll = cepa0_shortcut(blok)-1+l
-        l = det_cepa0_idx(ll)
-        ll = child_num(ll, J)
-      else
-        l = k
-        ll = kk
-      end if
-      
-        
-      if(.not. komoned) then
-        m = 0
-        m2 = 0
-          
-        do while(m < nlink(i_I) .and. m2 < nlink(J))
-          m += 1
-          m2 += 1
-          if(linked(m, i_I) < linked(m2, J)) then
-            m2 -= 1
-            cycle
-          else if(linked(m, i_I) > linked(m2, J)) then
-            m -= 1
-            cycle
-          end if
-          i = det_cepa0_idx(linked(m, i_I))
-          
-          if(h_cache(J,i) == 0.d0) cycle
-          if(h_cache(i_I,i) == 0.d0) cycle
-          
-          komon(0) += 1
-          kn = komon(0)
-          komon(kn) = i
-          
-          do i_state = 1,N_states
-            dkI = h_cache(J,i) * dij(i_I, i, i_state)
-            dleat(i_state, kn, 1) = dkI
-            dleat(i_state, kn, 2) = dkI
-
-            dkI = s2_cache(J,i) * dij(i_I, i, i_state)
-            dleat_s2(i_state, kn, 1) = dkI
-            dleat_s2(i_state, kn, 2) = dkI
-          end do
-
-        end do
-          
-        komoned = .true.
-      end if
-      
-      integer :: hpmin(2)
-      hpmin(1) = 2 - HP(1,k)
-      hpmin(2) = 2 - HP(2,k)
-
-      do m = 1, komon(0)
-        
-        i = komon(m)
-        if(HP(1,i) <= hpmin(1) .and. HP(2,i) <= hpmin(2) ) then
-          cycle
-        end if
-        
-        call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
-        if(.not. ok) cycle
-          
-        do i_state = 1, N_states 
-          contrib =  dij(i_I, k, i_state) * dleat(i_state, m, 2)
-          contrib_s2 =  dij(i_I, k, i_state) * dleat_s2(i_state, m, 2)
-          delta(i_state,ll,1) += contrib
-          delta_s2(i_state,ll,1) += contrib_s2
-          if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
-            delta(i_state,0,1) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state)
-            delta_s2(i_state,0,1) -= contrib_s2 * ci_inv(i_state) * psi_non_ref_coef(l,i_state)
-          endif
-          
-          if(I_i == J) cycle
-          contrib =  dij(J, l, i_state) * dleat(i_state, m, 1)
-          contrib_s2 =  dij(J, l, i_state) * dleat_s2(i_state, m, 1)
-          delta(i_state,kk,2) += contrib
-          delta_s2(i_state,kk,2) += contrib_s2
-          if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
-            delta(i_state,0,2) -= contrib * cj_inv(i_state) * psi_non_ref_coef(k,i_state)
-            delta_s2(i_state,0,2) -= contrib_s2 * cj_inv(i_state) * psi_non_ref_coef(k,i_state)
-          end if
-        enddo !i_state
-      end do ! while
-    end do ! kk
-
-      
-      call push_mrsc2_results(zmq_socket_push, I_i, J, delta, delta_s2, task_id) 
-      call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
- 
-!     end if
-    
-  enddo
-
-  deallocate(delta)
-
-  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 push_mrsc2_results(zmq_socket_push, I_i, J, delta, delta_s2, task_id) 
-  use f77_zmq 
-  implicit none 
-  BEGIN_DOC 
-! Push integrals in the push socket 
-  END_DOC 
-
-  integer, intent(in)            :: i_I, J
-  integer(ZMQ_PTR), intent(in)   :: zmq_socket_push 
-  double precision,intent(inout) :: delta(N_states, 0:N_det_non_ref, 2)
-  double precision,intent(inout) :: delta_s2(N_states, 0:N_det_non_ref, 2)
-  integer, intent(in)            :: task_id 
-  integer                        :: rc , i_state,  i, kk, li
-  integer,allocatable            :: idx(:,:)
-  integer                        :: n(2)
-  logical :: ok
-  
-  allocate(idx(N_det_non_ref,2))
-  rc = f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE) 
-  if (rc /= 4) then 
-    print *, irp_here,  'f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE)' 
-    stop 'error' 
-  endif 
-
-  rc = f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE) 
-  if (rc /= 4) then 
-    print *, irp_here,  'f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE)' 
-    stop 'error' 
-  endif 
-  
-  
-  do kk=1,2
-    n(kk)=0
-    if(kk == 1) li = nlink(j)
-    if(kk == 2) li = nlink(i_I)
-    do i=1, li
-      ok = .false.
-      do i_state=1,N_states
-        if(delta(i_state, i, kk) /= 0d0) then
-          ok = .true.
-          exit
-        end if
-      end do
-      
-      if(ok) then
-        n(kk) += 1
-!         idx(n,kk) = i
-        if(kk == 1) then
-          idx(n(1),1) = det_cepa0_idx(linked(i, J))
-        else
-          idx(n(2),2) = det_cepa0_idx(linked(i, i_I))
-        end if
-        
-        do i_state=1, N_states
-            delta(i_state, n(kk), kk) = delta(i_state, i, kk)
-        end do
-      end if
-    end do
-    
-    rc = f77_zmq_send( zmq_socket_push, n(kk), 4, ZMQ_SNDMORE) 
-    if (rc /= 4) then 
-      print *, irp_here,  'f77_zmq_send( zmq_socket_push, n, 4, ZMQ_SNDMORE)' 
-      stop 'error' 
-    endif 
-    
-    if(n(kk) /= 0) then
-      rc = f77_zmq_send( zmq_socket_push, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE)  ! delta(1,0,1) = delta_I   delta(1,0,2) = delta_J 
-      if (rc /=  (n(kk)+1)*8*N_states) then 
-        print *, irp_here,  'f77_zmq_send( zmq_socket_push, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
-        stop 'error' 
-      endif 
-
-      rc = f77_zmq_send( zmq_socket_push, delta_s2(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE)  ! delta_s2(1,0,1) = delta_I   delta_s2(1,0,2) = delta_J 
-      if (rc /=  (n(kk)+1)*8*N_states) then 
-        print *, irp_here,  'f77_zmq_send( zmq_socket_push, delta_s2, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
-        stop 'error' 
-      endif 
-      
-      rc = f77_zmq_send( zmq_socket_push, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) 
-      if (rc /=  n(kk)*4) then 
-        print *, irp_here,  'f77_zmq_send( zmq_socket_push, delta, 8*n(kk), ZMQ_SNDMORE)'
-        stop 'error' 
-      endif 
-    end if
-  end do
-  
-  
-  rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0) 
-  if (rc /= 4) then 
-    print *, irp_here,  'f77_zmq_send( zmq_socket_push, task_id, 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
-
-
-
-subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, delta_s2, task_id) 
-  use f77_zmq 
-  implicit none 
-  BEGIN_DOC 
-! Push integrals in the push socket 
-  END_DOC 
-
-  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull 
-  integer, intent(out)           :: i_I, J, n(2)
-  double precision, intent(inout) :: delta(N_states, 0:N_det_non_ref, 2)
-  double precision, intent(inout) :: delta_s2(N_states, 0:N_det_non_ref, 2)
-  integer, intent(out)           :: task_id 
-  integer                        :: rc , i, kk
-  integer,intent(inout) :: idx(N_det_non_ref,2)
-  logical :: ok
-  
-  rc = f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE) 
-  if (rc /= 4) then 
-    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE)' 
-    stop 'error' 
-  endif 
-  
-  rc = f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE) 
-  if (rc /= 4) then 
-    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE)' 
-    stop 'error' 
-  endif 
-  
-  do kk = 1, 2
-    rc = f77_zmq_recv( zmq_socket_pull, n(kk), 4, ZMQ_SNDMORE) 
-    if (rc /= 4) then 
-      print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, n, 4, ZMQ_SNDMORE)' 
-      stop 'error' 
-    endif 
-    
-    if(n(kk) /= 0) then
-      rc = f77_zmq_recv( zmq_socket_pull, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) 
-      if (rc /= (n(kk)+1)*8*N_states) then 
-        print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
-        stop 'error' 
-      endif 
-      
-      rc = f77_zmq_recv( zmq_socket_pull, delta_s2(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) 
-      if (rc /= (n(kk)+1)*8*N_states) then 
-        print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, delta_s2, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
-        stop 'error' 
-      endif 
-      
-      rc = f77_zmq_recv( zmq_socket_pull, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) 
-      if (rc /= n(kk)*4) then 
-        print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, idx(1,kk), n(kk)*4, ZMQ_SNDMORE)'
-        stop 'error' 
-      endif 
-    end if
-  end do
-  
-  rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) 
-  if (rc /= 4) then 
-    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, task_id, 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
-
-
-
-subroutine mrsc2_dressing_collector(delta_ii_,delta_ij_,delta_ii_s2_,delta_ij_s2_)
-  use f77_zmq
-  implicit none
-  BEGIN_DOC 
-! Collects results from the AO integral calculation 
-  END_DOC 
- 
-  double precision,intent(inout)               :: delta_ij_(N_states,N_det_non_ref,N_det_ref) 
-  double precision,intent(inout)               :: delta_ii_(N_states,N_det_ref)
-  double precision,intent(inout)               :: delta_ij_s2_(N_states,N_det_non_ref,N_det_ref) 
-  double precision,intent(inout)               :: delta_ii_s2_(N_states,N_det_ref)
-
-!   integer                        :: j,l
-  integer                        :: rc 
-   
-  double precision, allocatable  :: delta(:,:,:), delta_s2(:,:,:) 
-   
-  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*8                      :: control, accu 
-  integer                        :: task_id, more 
- 
-  integer                        :: I_i, J, l, i_state, n(2), kk
-  integer,allocatable :: idx(:,:)
-  
-  delta_ii_(:,:) = 0d0
-  delta_ij_(:,:,:) = 0d0
-  delta_ii_s2_(:,:) = 0d0
-  delta_ij_s2_(:,:,:) = 0d0
-
-  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() 
-  zmq_socket_pull = new_zmq_pull_socket() 
- 
-  allocate ( delta(N_states,0:N_det_non_ref,2), delta_s2(N_states,0:N_det_non_ref,2) ) 
-  
-  allocate(idx(N_det_non_ref,2))
-  more = 1 
-  do while (more == 1) 
-         
-    call pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, delta_s2, task_id)
-    
-
-      do l=1, n(1)
-        do i_state=1,N_states
-          delta_ij_(i_state,idx(l,1),i_I) += delta(i_state,l,1)
-          delta_ij_s2_(i_state,idx(l,1),i_I) += delta_s2(i_state,l,1)
-        end do
-      end do
-      
-      do l=1, n(2)
-        do i_state=1,N_states
-          delta_ij_(i_state,idx(l,2),J) += delta(i_state,l,2)
-          delta_ij_s2_(i_state,idx(l,2),J) += delta_s2(i_state,l,2)
-        end do
-      end do
-
-    
-       if(n(1) /= 0) then 
-       do i_state=1,N_states
-         delta_ii_(i_state,i_I) += delta(i_state,0,1)
-         delta_ii_s2_(i_state,i_I) += delta_s2(i_state,0,1)
-       end do
-       end if
-
-      if(n(2) /= 0) then 
-       do i_state=1,N_states
-         delta_ii_(i_state,J) += delta(i_state,0,2)
-         delta_ii_s2_(i_state,J) += delta_s2(i_state,0,2)
-       end do
-       end if
-
-
-    if (task_id /= 0) then 
-      call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) 
-    endif 
-
- 
-  enddo 
-  deallocate( delta, delta_s2 )
- 
-  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
-  call end_zmq_pull_socket(zmq_socket_pull)
-
-end
-
-
-
-
- BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref,N_det_ref) ]
-&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_states,N_det_ref) ]
-&BEGIN_PROVIDER [ double precision, delta_ij_s2_old, (N_states,N_det_non_ref,N_det_ref) ]
-&BEGIN_PROVIDER [ double precision, delta_ii_s2_old, (N_states,N_det_ref) ]
-  implicit none
-  
-  integer                         :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2
-  integer                         :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, nex, nzer, ntot
-!   integer, allocatable            :: linked(:,:), blokMwen(:, :), nlink(:)
-  logical                         :: ok
-  double precision                :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI(N_states), HkI, ci_inv(N_states), dia_hla(N_states)
-  double precision                :: contrib, wall, iwall ! , searchance(N_det_ref)
-  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
-  integer(bit_kind)               :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt
-  integer, external               :: get_index_in_psi_det_sorted_bit, searchDet, detCmp
-  logical, external               :: is_in_wavefunction, isInCassd, detEq
-  character*(512)                :: task 
-  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket 
-  
-  integer :: KKsize = 1000000
-  
-  
-  call new_parallel_job(zmq_to_qp_run_socket,'mrsc2')
-
-
-  call wall_time(iwall)
-!   allocate(linked(N_det_non_ref, N_det_ref), blokMwen(N_det_non_ref, N_det_ref), nlink(N_det_ref))
-  
-
-!   searchance = 0d0
-!   do J = 1, N_det_ref
-!     nlink(J) = 0
-!     do blok=1,cepa0_shortcut(0)
-!     do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
-!       call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int)
-!       if(degree <= 2) then
-!         nlink(J) += 1
-!         linked(nlink(J),J) = k
-!         blokMwen(nlink(J),J) = blok
-!         searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok)))
-!       end if
-!     end do
-!     end do
-!   end do
-
-     
-     
-!   stop
-  nzer = 0
-  ntot = 0
-  do nex = 3, 0, -1
-    print *, "los ",nex
-    do I_s = N_det_ref, 1, -1
-!         if(mod(I_s,1) == 0) then
-!           call wall_time(wall)
-!           wall = wall-iwall
-!           print *, I_s, "/", N_det_ref, wall * (dfloat(N_det_ref) / dfloat(I_s)), wall, wall * (dfloat(N_det_ref) / dfloat(I_s))-wall
-!         end if
-
-
-      do J_s = 1, I_s
-      
-        call get_excitation_degree(psi_ref(1,1,J_s), psi_ref(1,1,I_s), degree, N_int)
-        if(degree /= nex) cycle
-        if(nex == 3) nzer = nzer + 1
-        ntot += 1
-!           if(degree > 3) then
-!             deg += 1
-!             cycle
-!           else if(degree == -10) then
-!             KKsize = 100000
-!           else
-!             KKsize = 1000000
-!           end if
-        
-        
-        
-        if(searchance(I_s) < searchance(J_s)) then
-          i_I = I_s
-          J = J_s
-        else
-          i_I = J_s
-          J = I_s
-        end if
-        
-        KKsize = nlink(1)
-        if(nex == 0) KKsize = int(float(nlink(1)) / float(nlink(i_I)) * (float(nlink(1)) / 64d0))
-        
-        !if(KKsize == 0) stop "ZZEO"
-        
-        do kk = 1 , nlink(i_I), KKsize
-          write(task,*) I_i, J, kk, int(min(kk+KKsize-1, nlink(i_I)))
-          call add_task_to_taskserver(zmq_to_qp_run_socket,task)
-        end do
-        
-  !       do kk = 1 , nlink(i_I)
-  !         k = linked(kk,i_I)
-  !         blok = blokMwen(kk,i_I)
-  !         write(task,*) I_i, J, k, blok
-  !         call add_task_to_taskserver(zmq_to_qp_run_socket,task)
-  !         
-  !       enddo !kk
-      enddo !J
-          
-    enddo !I
-  end do ! nex
-  print *, "tasked"
-!  integer(ZMQ_PTR)               ∷ collector_thread
-!  external                       ∷ ao_bielec_integrals_in_map_collector
-!  rc = pthread_create(collector_thread, mrsc2_dressing_collector)
-  print *, nzer, ntot, float(nzer) / float(ntot)
-  provide nproc
-  !$OMP PARALLEL DEFAULT(none)  SHARED(delta_ii_old,delta_ij_old,delta_ii_s2_old,delta_ij_s2_old)  PRIVATE(i) NUM_THREADS(nproc+1)
-      i = omp_get_thread_num()
-      if (i==0) then
-        call mrsc2_dressing_collector(delta_ii_old,delta_ij_old,delta_ii_s2_old,delta_ij_s2_old)
-      else
-        call mrsc2_dressing_slave_inproc(i)
-      endif
-  !$OMP END PARALLEL
-
-!  rc = pthread_join(collector_thread)
-  call end_parallel_job(zmq_to_qp_run_socket, 'mrsc2')
-    
-
-END_PROVIDER
-
-
-

plugins/mrcc_selected/ezfio_interface.irp.f

@@ -1,61 +0,0 @@
-! 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
-
-
-BEGIN_PROVIDER [ double precision, thresh_dressed_ci  ]
-  implicit none
-  BEGIN_DOC
-! Threshold on the convergence of the dressed CI energy
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  
-  call ezfio_has_mrcc_selected_thresh_dressed_ci(has)
-  if (has) then
-    call ezfio_get_mrcc_selected_thresh_dressed_ci(thresh_dressed_ci)
-  else
-    print *, 'mrcc_selected/thresh_dressed_ci not found in EZFIO file'
-    stop 1
-  endif
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, n_it_max_dressed_ci  ]
-  implicit none
-  BEGIN_DOC
-! Maximum number of dressed CI iterations
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  
-  call ezfio_has_mrcc_selected_n_it_max_dressed_ci(has)
-  if (has) then
-    call ezfio_get_mrcc_selected_n_it_max_dressed_ci(n_it_max_dressed_ci)
-  else
-    print *, 'mrcc_selected/n_it_max_dressed_ci not found in EZFIO file'
-    stop 1
-  endif
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, lambda_type  ]
-  implicit none
-  BEGIN_DOC
-! lambda type
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  
-  call ezfio_has_mrcc_selected_lambda_type(has)
-  if (has) then
-    call ezfio_get_mrcc_selected_lambda_type(lambda_type)
-  else
-    print *, 'mrcc_selected/lambda_type not found in EZFIO file'
-    stop 1
-  endif
-
-END_PROVIDER

plugins/mrcc_selected/mrcc_selected.irp.f

@@ -1,19 +0,0 @@
-program mrsc2sub
-  implicit none
-  double precision, allocatable  :: energy(:)
-  allocate (energy(N_states))
-  
-  !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
-  mrmode = 3
-  
-  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
-    call run_pt2(N_states,energy)
-  endif
-  deallocate(energy)
-end
-

plugins/mrcc_selected/mrcepa0_general.irp.f

@@ -1,245 +0,0 @@
-
-
-subroutine run(N_st,energy)
-  implicit none
-  
-  integer, intent(in) :: N_st 
-  double precision, intent(out) :: energy(N_st) 
-
-  integer :: i,j
-
-  double precision :: E_new, E_old, delta_e
-  integer :: iteration
-  double precision :: E_past(4)
-  
-  integer :: n_it_mrcc_max
-  double precision :: thresh_mrcc
-  double precision, allocatable :: lambda(:)
-  allocate (lambda(N_states))
-
-  
-  thresh_mrcc = thresh_dressed_ci
-  n_it_mrcc_max = n_it_max_dressed_ci
-
-  if(n_it_mrcc_max == 1) then
-    do j=1,N_states_diag
-      do i=1,N_det
-        psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
-      enddo
-    enddo
-    SOFT_TOUCH psi_coef ci_energy_dressed
-    call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
-    call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
-    call save_wavefunction
-    energy(:) = ci_energy_dressed(:)
-  else
-    E_new = 0.d0
-    delta_E = 1.d0
-    iteration = 0
-    lambda = 1.d0
-    do while (delta_E > thresh_mrcc)
-      iteration += 1
-      print *,  '===========================' 
-      print *,  'MRCEPA0 Iteration', iteration
-      print *,  '===========================' 
-      print *,  ''
-      E_old = sum(ci_energy_dressed)
-      call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy")
-      call diagonalize_ci_dressed(lambda)
-      E_new = sum(ci_energy_dressed)
-      delta_E = dabs(E_new - E_old)
-      call save_wavefunction
-      call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
-      if (iteration >= n_it_mrcc_max) then
-        exit
-      endif
-    enddo
-    call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy")
-    energy(:) = ci_energy_dressed(:)
-  endif
-end
-
-
-subroutine print_cas_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)
-  enddo
-  call write_double(6,ci_energy(1),"Initial CI energy")
-
-end
-
-
-
-
-subroutine run_pt2_old(N_st,energy) 
-  implicit none 
-  integer :: i,j,k 
-  integer, intent(in)          :: N_st 
-  double precision, intent(in) :: energy(N_st) 
-  double precision :: pt2_redundant(N_st), pt2(N_st)
-  double precision :: norm_pert(N_st),H_pert_diag(N_st)
-  
-  pt2_redundant = 0.d0 
-  pt2 = 0d0
-  !if(lambda_mrcc_pt2(0)  == 0) return
-  
-  print*,'Last iteration only to compute the PT2' 
-  
-  print * ,'Computing the redundant PT2 contribution'
-
-  if (mrmode == 1) then
-
-    N_det_generators = lambda_mrcc_kept(0) 
-    N_det_selectors = lambda_mrcc_kept(0) 
-
-    do i=1,N_det_generators
-      j = lambda_mrcc_kept(i)
-      do k=1,N_int
-        psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
-        psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
-        psi_selectors(k,1,i) = psi_non_ref(k,1,j)
-        psi_selectors(k,2,i) = psi_non_ref(k,2,j)
-      enddo
-      do k=1,N_st
-        psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
-        psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
-      enddo
-    enddo
-
-  else
-
-    N_det_generators = N_det_non_ref
-    N_det_selectors = N_det_non_ref
-
-    do i=1,N_det_generators
-      j = i
-      do k=1,N_int
-        psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
-        psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
-        psi_selectors(k,1,i) = psi_non_ref(k,1,j)
-        psi_selectors(k,2,i) = psi_non_ref(k,2,j)
-      enddo
-      do k=1,N_st
-        psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
-        psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
-      enddo
-    enddo
-
-  endif
-  
-  SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
-  SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
-
-  call H_apply_mrcepa_PT2(pt2_redundant, norm_pert, H_pert_diag,  N_st) 
-  
-  print * ,'Computing the remaining contribution'
-
-  threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
-  threshold_generators = max(threshold_generators,threshold_generators_pt2)
-
-  N_det_generators = N_det_non_ref + N_det_ref
-  N_det_selectors = N_det_non_ref + N_det_ref
-  
-  psi_det_generators(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
-  psi_selectors(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
-  psi_coef_generators(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
-  psi_selectors_coef(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
-  
-  do i=N_det_ref+1,N_det_generators
-    j = i-N_det_ref
-    do k=1,N_int
-      psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
-      psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
-      psi_selectors(k,1,i) = psi_non_ref(k,1,j)
-      psi_selectors(k,2,i) = psi_non_ref(k,2,j)
-    enddo
-    do k=1,N_st
-      psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
-      psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
-    enddo
-  enddo
-  
-  SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
-  SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
-
-  call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag,  N_st) 
- 
-
-  print *, "Redundant PT2 :",pt2_redundant
-  print *, "Full      PT2 :",pt2
-  print *, lambda_mrcc_kept(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1)
-  pt2 = pt2 - pt2_redundant
-  
-  print *,  'Final step' 
-  print *,  'N_det    = ', N_det 
-  print *,  'N_states = ', N_states 
-  print *,  'PT2      = ', pt2 
-  print *,  'E        = ', energy 
-  print *,  'E+PT2    = ', energy+pt2 
-  print *,  '-----' 
- 
-
-  call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1))
-
-end 
-
-subroutine run_pt2(N_st,energy) 
-  implicit none 
-  integer :: i,j,k 
-  integer, intent(in)          :: N_st 
-  double precision, intent(in) :: energy(N_st) 
-  double precision :: pt2(N_st)
-  double precision :: norm_pert(N_st),H_pert_diag(N_st)
-  
-  pt2 = 0d0
-  !if(lambda_mrcc_pt2(0)  == 0) return
-  
-  print*,'Last iteration only to compute the PT2' 
-  
-  N_det_generators = N_det_cas
-  N_det_selectors = N_det_non_ref
-
-  do i=1,N_det_generators
-    do k=1,N_int
-      psi_det_generators(k,1,i) = psi_ref(k,1,i)
-      psi_det_generators(k,2,i) = psi_ref(k,2,i)
-    enddo
-    do k=1,N_st
-      psi_coef_generators(i,k) = psi_ref_coef(i,k)
-    enddo
-  enddo
-  do i=1,N_det
-    do k=1,N_int
-      psi_selectors(k,1,i) = psi_det_sorted(k,1,i)
-      psi_selectors(k,2,i) = psi_det_sorted(k,2,i)
-    enddo
-    do k=1,N_st
-      psi_selectors_coef(i,k) = psi_coef_sorted(i,k)
-    enddo
-  enddo
-
-  SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
-  SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
-
-  call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag,  N_st) 
-  
-!  call ezfio_set_full_ci_energy_pt2(energy+pt2)
-
-  print *,  'Final step' 
-  print *,  'N_det    = ', N_det 
-  print *,  'N_states = ', N_states 
-  print *,  'PT2      = ', pt2 
-  print *,  'E        = ', energy 
-  print *,  'E+PT2    = ', energy+pt2 
-  print *,  '-----' 
-
-  call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1))
-
-end 
-