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

.gitignore

@@ -1,4 +1,5 @@
 quantum_package.rc
+config/ifort.cfg
 quantum_package_static.tar.gz
 build.ninja
 .ninja_log

.travis.yml

@@ -24,7 +24,7 @@ python:
 
 script: 
   - ./configure  --production ./config/gfortran.cfg
-  - source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles
+  - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD mrcepa0 All_singles
   - source ./quantum_package.rc ; ninja
   - source ./quantum_package.rc ; cd ocaml ; make ; cd -
   - source ./quantum_package.rc ; cd tests ; ./run_tests.sh #-v

README.md

@@ -7,11 +7,14 @@ Set of quantum chemistry programs and libraries.
 
 For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>), or below for the installation instructions.
 
+
+
 Demo
 ====
 
 [![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_295x166.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo")
 
+[![Frozen-core Full-CI energy of Ti](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png)](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png "Convergence of Ti in cc-pv{DTQ}Z")
 
 # Installation
 
@@ -155,7 +158,7 @@ Program exited with code 139.
 
 #### Why ?
 
-It's caused when we call the DGEM routine of LAPACK. 
+It's caused when we call the DGEMM routine of LAPACK. 
 
 ##### Fix
 

config/gfortran.cfg

@@ -35,7 +35,7 @@ OPENMP  : 1          ; Append OpenMP flags
 #                           -ffast-math and the Fortran-specific
 #                           -fno-protect-parens and -fstack-arrays.
 [OPT]
-FCFLAGS : -Ofast 
+FCFLAGS : -Ofast
 
 # Profiling flags
 #################

config/gfortran_avx.cfg

@@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -ffree-line-length-none -I . -mavx
+FC           : gfortran -ffree-line-length-none -I . -mavx -g
 LAPACK_LIB   : -llapack -lblas
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32

config/gfortran_debug.cfg

@@ -51,7 +51,7 @@ FCFLAGS : -Ofast
 # -g           : Extra debugging information
 #
 [DEBUG]
-FCFLAGS : -g -pedantic -msse4.2 
+FCFLAGS : -g -msse4.2 
 
 # OpenMP flags
 #################

config/ifort.cfg

@@ -32,14 +32,14 @@ OPENMP  : 1          ; Append OpenMP flags
 #
 [OPT]
 FC       : -traceback
-FCFLAGS  : -xSSE4.2 -O2 -ip -ftz -g -traceback
+FCFLAGS  : -xSSE4.2  -O2 -ip -ftz -g
 
 # Profiling flags
 #################
 #
 [PROFILE]
 FC       : -p -g -traceback
-FCFLAGS  : -xSSE4.2 -O2 -ip -ftz 
+FCFLAGS  : -xSSE4.2 -O2 -ip -ftz
 
 # Debugging flags
 #################
@@ -52,7 +52,7 @@ FCFLAGS  : -xSSE4.2 -O2 -ip -ftz
 #
 [DEBUG]
 FC      : -g -traceback
-FCFLAGS : -xSSE2 -C  
+FCFLAGS : -xSSE2 -C -fpe0 
 IRPF90_FLAGS : --openmp
 
 # OpenMP flags

config/sse4_avx2.cfg

@@ -0,0 +1,62 @@
+# Common flags
+##############
+#
+# -mkl=[parallel|sequential] : Use the MKL library
+# --ninja                 : Allow the utilisation of ninja. It is mandatory !
+# --align=32                 : Align all provided arrays on a 32-byte boundary
+#
+[COMMON]
+FC           : ifort
+LAPACK_LIB   : -mkl=parallel
+IRPF90       : irpf90
+IRPF90_FLAGS : --ninja --align=32
+
+# Global options
+################
+#
+# 1 : Activate
+# 0 : Deactivate
+# 
+[OPTION]
+MODE    : OPT        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
+CACHE   : 1          ; Enable cache_compile.py
+OPENMP  : 1          ; Append OpenMP flags
+
+# Optimization flags
+####################
+#
+# -xHost                     : Compile a binary optimized for the current architecture
+# -O2                        : O3 not better than O2.
+# -ip                        : Inter-procedural optimizations
+# -ftz                       : Flushes denormal results to zero
+#
+[OPT]
+FCFLAGS  : -axSSE4.2,AVX,CORE-AVX2 -O2 -ip -ftz -g -traceback
+
+# Profiling flags
+#################
+#
+[PROFILE]
+FC       : -p -g
+FCFLAGS  : -xSSE4.2 -O2 -ip -ftz
+
+# Debugging flags
+#################
+#
+# -traceback   : Activate backtrace on runtime
+# -fpe0        : All floating point exaceptions
+# -C           : Checks uninitialized variables,  array subscripts, etc...
+# -g           : Extra debugging information
+# -xSSE2       : Valgrind needs a very simple x86 executable
+#
+[DEBUG]
+FC      : -g -traceback
+FCFLAGS : -xSSE2 -C  -fpe0 
+
+# OpenMP flags
+#################
+#
+[OPENMP]
+FC           : -openmp
+IRPF90_FLAGS : --openmp 
+

configure

@@ -142,7 +142,7 @@ ezfio = Info(
     default_path=join(QP_ROOT_INSTALL, "EZFIO"))
 
 zeromq = Info(
-    url='https://github.com/zeromq/zeromq4-1/releases/download/v4.1.4/zeromq-4.1.4.tar.gz',
+    url='https://github.com/zeromq/zeromq4-1/releases/download/v4.1.5/zeromq-4.1.5.tar.gz',
     description=' ZeroMQ',
     default_path=join(QP_ROOT_LIB, "libzmq.a"))
 
@@ -166,7 +166,7 @@ d_info = dict()
 
 for m in ["ocaml", "m4", "curl", "zlib", "patch", "irpf90", "docopt",
           "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz",
-          "zeromq", "f77zmq","bats"]:
+          "zeromq", "f77zmq","bats" ]:
     exec ("d_info['{0}']={0}".format(m))
 
 
@@ -543,7 +543,6 @@ def recommendation():
     print ""
     print "Finally :"
     print "   ninja"
-    print "   make -C ocaml"
     print ""
     print "You can install more plugin with the qp_module.py install command" 
     print "PS : For more info on compiling the code, read the README.md"

data/basis/aug-cc-pcv5z

@@ -705,3 +705,5 @@ H   1
   1     21.1040000              1.0000000        
 H   1
   1      0.7420000              1.0000000
+
+

data/basis/aug-cc-pcvdz

@@ -893,3 +893,5 @@ D   1
   1     11.4590000              1.0000000        
 D   1
   1      0.2400000              1.0000000
+
+

data/basis/aug-cc-pcvqz

@@ -1594,3 +1594,5 @@ G   1
   1     17.2430000              1.0000000        
 G   1
   1      0.4590000              1.0000000
+
+

data/basis/aug-cc-pcvtz

@@ -1224,3 +1224,5 @@ F   1
   1     13.6740000              1.0000000        
 F   1
   1      0.4060000              1.0000000
+
+

data/basis/aug-cc-pv5z

@@ -7065,3 +7065,5 @@ H   1
   1      0.9303000              1.0000000        
 H   1
   1      0.5800000              1.0000000
+
+

data/basis/aug-cc-pv6z

@@ -1515,3 +1515,5 @@ I   1
   1      1.5066000              1.0000000        
 I   1
   1      0.9926000              1.0000000
+
+

data/basis/aug-cc-pvdz

@@ -3485,3 +3485,5 @@ D   1
   1      0.5030000              1.0000000        
 D   1
   1      0.2155000              1.0000000
+
+

data/basis/aug-cc-pvqz

@@ -5685,3 +5685,5 @@ G   1
   1      0.7395000              1.0000000        
 G   1
   1      0.3590000              1.0000000
+
+

data/basis/aug-cc-pvtz

@@ -4421,3 +4421,5 @@ F   1
   1      0.6622000              1.0000000        
 F   1
   1      0.3280000              1.0000000
+
+

data/basis/cc-pcv5z

@@ -1614,3 +1614,5 @@ G   1
   1      0.3023000              1.0000000        
 H   1
   1      0.2534000              1.0000000
+
+

data/basis/cc-pcv6z

@@ -1515,3 +1515,5 @@ I   1
   1      1.5066000              1.0000000        
 I   1
   1     24.5369000              1.0000000
+
+

data/basis/cc-pcvdz

@@ -905,3 +905,5 @@ D   1
   1      0.0537000              1.0000000        
 D   1
   1      1.3743000              1.0000000
+
+

data/basis/cc-pcvqz

@@ -1611,3 +1611,5 @@ G   1
   1      0.1466000              1.0000000        
 G   1
   1      1.5908000              1.0000000
+
+

data/basis/cc-pcvtz

@@ -1246,3 +1246,5 @@ F   1
   1      0.1509000              1.0000000        
 F   1
   1      1.3909000              1.0000000
+
+

data/basis/cc-pv5z

@@ -7212,3 +7212,5 @@ G   1
   1      1.1040000              1.0000000        
 H   1
   1      0.9303000              1.0000000
+
+

data/basis/cc-pv6z

@@ -1323,3 +1323,5 @@ H   1
   1      0.8871000              1.0000000        
 I   1
   1      1.5066000              1.0000000
+
+

data/basis/cc-pvdz

@@ -3367,3 +3367,5 @@ D   5
   5      1.5075240              0.2667560        
 D   1
   1      0.5030000              1.0000000
+
+

data/basis/cc-pvqz

@@ -5482,3 +5482,5 @@ F   1
   1      0.9557000              1.0000000        
 G   1
   1      0.7395000              1.0000000
+
+

data/basis/cc-pvtz

@@ -4269,3 +4269,5 @@ D   1
   1      0.3006000              1.0000000        
 F   1
   1      0.6622000              1.0000000
+
+

data/basis/chipman-dzp

@@ -166,3 +166,5 @@ D   1
   1      0.6650000              1.0000000        
 D   1
   1      2.6600000              1.0000000        
+
+

data/basis/v5z-bfd

@@ -1017,3 +1017,5 @@ F 1
 G 1
   1 0.623669      1.000000
 
+
+

data/basis/vdz-ano-bfd

@@ -586,3 +586,5 @@ S 1
 P 1
   1 1.275000      1.000000
 
+
+

data/basis/vdz-bfd

@@ -997,3 +997,5 @@ P 8
   7 9.063386      -0.224631
   8 16.737180      0.098422
 
+
+

data/basis/vqz-ano-bfd

@@ -2409,3 +2409,5 @@ G 1
 H 1
   1 3.164456      1.000000
 
+
+

data/basis/vqz-bfd

@@ -809,3 +809,5 @@ D 1
 F 1
   1 1.021427      1.000000
 
+
+

data/basis/vtz-ano-bfd

@@ -1850,3 +1850,5 @@ F 1
 G 1
   1 2.775762      1.000000
 
+
+

data/basis/vtz-bfd

@@ -1279,3 +1279,5 @@ P 1
 D 1
   1 1.913792      1.000000
 
+
+

data/pseudo/tn_df

@@ -780,7 +780,7 @@ Ar GEN 10 2
         -1386.79918148    2    4.23753203    
          1350.57102634    2    6.12344921    
 
-Ag GEN 36 2 
+Ag GEN 36 2
   6
            11.00000000    1    7.02317516    
           178.71479273    2    1.36779344    

ocaml/Element.mli

@@ -19,3 +19,4 @@ val to_charge : t -> Charge.t
 val of_charge : Charge.t -> t
 val covalent_radius : t -> Qptypes.Positive_float.t
 val vdw_radius : t -> Qptypes.Positive_float.t
+val mass : t -> Qptypes.Positive_float.t

ocaml/Id.ml

@@ -1,26 +1,22 @@
-open Core.Std
-
-module Id : sig
-  type t
-  val of_int : int -> t
-  val to_int : t -> int
-  val of_string : string -> t
-  val to_string : t -> string
-  val increment : t -> t
-  val decrement : t -> t
-end
-= struct
+module Id = struct
   type t = int
+
   let of_int x = 
     assert (x>0); x
+
   let to_int x = x
+
   let of_string x = 
-    Int.of_string x
+    int_of_string x
     |> of_int
+
   let to_string x =
-    Int.to_string x
+    string_of_int x
+
   let increment x = x + 1
   let decrement x = x - 1
+
+  let compare = compare
 end
 
 module Task = struct

ocaml/Id.mli

@@ -0,0 +1,23 @@
+module Id :
+  sig
+    type t
+    val of_int : int -> t
+    val to_int : t -> int
+    val of_string : string -> t
+    val to_string : t -> string
+    val increment : t -> t
+    val decrement : t -> t
+    val compare   : t -> t -> int
+  end
+
+
+module Task :
+  sig
+    include (module type of Id)
+  end
+
+
+module Client :
+  sig
+    include (module type of Id)
+  end

ocaml/Input_determinants_by_hand.ml

@@ -93,23 +93,6 @@ end = struct
   ;;
 
 
-  let read_n_states_diag () =
-    if not (Ezfio.has_determinants_n_states_diag ()) then
-      read_n_states ()
-      |> States_number.to_int
-      |> Ezfio.set_determinants_n_states_diag 
-    ;
-    Ezfio.get_determinants_n_states_diag ()
-    |> States_number.of_int
-  ;;
-
-  let write_n_states_diag ~n_states n =
-    let n_states = States_number.to_int n_states
-    and n = States_number.to_int n
-    in
-    Ezfio.set_determinants_n_states_diag (max n_states n)
-  ;;
-
   let read_expected_s2 () =
     if not (Ezfio.has_determinants_expected_s2 ()) then
       begin

ocaml/Makefile

@@ -79,7 +79,7 @@ git:
 ${QP_ROOT}/install/EZFIO/Ocaml/ezfio.ml: 
 	$(NINJA) -C ${QP_ROOT}/install/EZFIO
 
-Input_auto_generated.ml qp_edit.ml:
+Input_auto_generated.ml qp_edit.ml: $(filter-out Input_auto_generated.ml, $(wildcard Input_*.ml))
 	ei_handler.py ocaml_global
 
 clean: 

ocaml/Message.ml

@@ -248,16 +248,20 @@ end
 (** GetTaskReply : Reply to the GetTask message *)
 module GetTaskReply_msg : sig
   type t  
-  val create : task_id:Id.Task.t -> task:string -> t
+  val create : task_id:Id.Task.t option -> task:string option -> t
   val to_string : t -> string
 end = struct
   type t =
-  { task_id: Id.Task.t ;
-    task   : string ;
+  { task_id: Id.Task.t option ;
+    task   : string option ;
   }
   let create ~task_id ~task = { task_id ; task }
   let to_string x =
-    Printf.sprintf "get_task_reply %d %s" (Id.Task.to_int x.task_id) x.task
+    match x.task_id, x.task with
+    | Some task_id, Some task -> 
+      Printf.sprintf "get_task_reply %d %s" (Id.Task.to_int task_id) task
+    | _ ->
+      Printf.sprintf "get_task_reply 0"
 end
 
 (** GetPsi : get the current variational wave function *)
@@ -288,13 +292,14 @@ module Psi : sig
       n_det_selectors : Strictly_positive_int.t option;
       psi_det   :  string                    ;
       psi_coef  :  string                    ;
+      energy    :  string;
   }
   val create : n_state:Strictly_positive_int.t
      -> n_det:Strictly_positive_int.t 
      -> psi_det_size:Strictly_positive_int.t 
      -> n_det_generators:Strictly_positive_int.t option
      -> n_det_selectors:Strictly_positive_int.t option
-     -> psi_det:string -> psi_coef:string -> t
+     -> psi_det:string -> psi_coef:string -> energy:string -> t
 end = struct
   type t = 
   {
@@ -305,14 +310,16 @@ end = struct
       n_det_selectors : Strictly_positive_int.t option;
       psi_det   :  string                    ;
       psi_coef  :  string                    ;
+      energy    :  string                    ;
   }
   let create ~n_state ~n_det ~psi_det_size
-    ~n_det_generators ~n_det_selectors ~psi_det ~psi_coef =
+    ~n_det_generators ~n_det_selectors ~psi_det ~psi_coef 
+    ~energy =
     assert (Strictly_positive_int.to_int n_det <=
             Strictly_positive_int.to_int psi_det_size );
     {  n_state; n_det ; psi_det_size ;
        n_det_generators ; n_det_selectors ;
-       psi_det ; psi_coef }
+       psi_det ; psi_coef ; energy }
 end
 
 (** GetPsiReply_msg : Reply to the GetPsi message *)
@@ -329,19 +336,6 @@ end = struct
     psi       :  Psi.t }
   let create ~client_id ~psi =
     {  client_id ; psi }
-  let to_string_list x =
-    let g, s = 
-      match x.psi.Psi.n_det_generators, x.psi.Psi.n_det_selectors with
-      | Some g, Some s -> Strictly_positive_int.to_int g, Strictly_positive_int.to_int s
-      | _ -> -1, -1
-    in
-    [ Printf.sprintf "get_psi_reply %d %d %d %d %d %d"
-      (Id.Client.to_int x.client_id)
-      (Strictly_positive_int.to_int x.psi.Psi.n_state)
-      (Strictly_positive_int.to_int x.psi.Psi.n_det) 
-      (Strictly_positive_int.to_int x.psi.Psi.psi_det_size)
-      g s ;
-      x.psi.Psi.psi_det ; x.psi.Psi.psi_coef ]
   let to_string x =
     let g, s = 
       match x.psi.Psi.n_det_generators, x.psi.Psi.n_det_selectors with
@@ -354,6 +348,9 @@ end = struct
       (Strictly_positive_int.to_int x.psi.Psi.n_det) 
       (Strictly_positive_int.to_int x.psi.Psi.psi_det_size) 
       g s
+  let to_string_list x =
+    [ to_string x ;
+      x.psi.Psi.psi_det ; x.psi.Psi.psi_coef ; x.psi.Psi.energy ]
 end
 
 
@@ -375,7 +372,8 @@ module PutPsi_msg : sig
      psi_det:string option ->
      psi_coef:string option ->
      n_det_generators: string option -> 
-     n_det_selectors:string option ->  t
+     n_det_selectors:string option ->
+     energy:string option -> t
   val to_string_list : t -> string list
   val to_string : t -> string 
 end = struct
@@ -388,7 +386,7 @@ end = struct
     n_det_selectors  : Strictly_positive_int.t option;
     psi       :  Psi.t option }
   let create ~client_id ~n_state ~n_det ~psi_det_size ~psi_det ~psi_coef 
-    ~n_det_generators ~n_det_selectors  =
+    ~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 ,
@@ -407,45 +405,19 @@ end = struct
       | _ -> None, None
     in
     let psi =
-      match (psi_det, psi_coef) with
-      | (Some psi_det, Some psi_coef) ->
+      match (psi_det, psi_coef, energy) with
+      | (Some psi_det, Some psi_coef, Some energy) ->
         Some (Psi.create ~n_state ~n_det ~psi_det_size ~psi_det
-          ~psi_coef ~n_det_generators ~n_det_selectors)
+          ~psi_coef ~n_det_generators ~n_det_selectors ~energy)
       | _ -> None
     in
     { client_id = Id.Client.of_string client_id ;
       n_state ; n_det ; psi_det_size ; n_det_generators ;
       n_det_selectors ; psi }
-  let to_string_list x =
-    match x.n_det_generators, x.n_det_selectors, x.psi with
-    | Some g, Some s, Some psi ->
-      [ Printf.sprintf "put_psi %d %d %d %d %d %d"
-        (Id.Client.to_int x.client_id)
-        (Strictly_positive_int.to_int x.n_state)
-        (Strictly_positive_int.to_int x.n_det) 
-        (Strictly_positive_int.to_int x.psi_det_size)  
-        (Strictly_positive_int.to_int g)  
-        (Strictly_positive_int.to_int s) ; 
-          psi.Psi.psi_det ; psi.Psi.psi_coef ]
-    | Some g, Some s, None ->
-      [ Printf.sprintf "put_psi %d %d %d %d %d %d"
-        (Id.Client.to_int x.client_id)
-        (Strictly_positive_int.to_int x.n_state)
-        (Strictly_positive_int.to_int x.n_det) 
-        (Strictly_positive_int.to_int x.psi_det_size)  
-        (Strictly_positive_int.to_int g)  
-        (Strictly_positive_int.to_int s) ; 
-          "None" ; "None" ]
-    | _ ->
-      [ Printf.sprintf "put_psi %d %d %d %d -1 -1"
-        (Id.Client.to_int x.client_id)
-        (Strictly_positive_int.to_int x.n_state)
-        (Strictly_positive_int.to_int x.n_det) 
-        (Strictly_positive_int.to_int x.psi_det_size) ;
-          "None" ; "None" ]
+
   let to_string x =
-    match x.n_det_generators, x.n_det_selectors, x.psi with
-    | Some g, Some s, Some psi ->
+    match x.n_det_generators, x.n_det_selectors with
+    | Some g, Some s ->
       Printf.sprintf "put_psi %d %d %d %d %d %d"
         (Id.Client.to_int x.client_id)
         (Strictly_positive_int.to_int x.n_state)
@@ -453,21 +425,20 @@ end = struct
         (Strictly_positive_int.to_int x.psi_det_size)  
         (Strictly_positive_int.to_int g)  
         (Strictly_positive_int.to_int s) 
-    | Some g, Some s, None ->
-      Printf.sprintf "put_psi %d %d %d %d %d %d"
-        (Id.Client.to_int x.client_id)
-        (Strictly_positive_int.to_int x.n_state)
-        (Strictly_positive_int.to_int x.n_det) 
-        (Strictly_positive_int.to_int x.psi_det_size)  
-        (Strictly_positive_int.to_int g)  
-        (Strictly_positive_int.to_int s) 
-    | _, _, _ ->
+    | _, _ ->
       Printf.sprintf "put_psi %d %d %d %d %d %d"
         (Id.Client.to_int x.client_id)
         (Strictly_positive_int.to_int x.n_state)
         (Strictly_positive_int.to_int x.n_det) 
         (Strictly_positive_int.to_int x.psi_det_size)  
         (-1) (-1)
+
+  let to_string_list x =
+    match x.psi with
+    | Some psi ->
+      [ to_string x ; psi.Psi.psi_det ; psi.Psi.psi_coef ; psi.Psi.energy ]
+    | None ->
+      [ to_string x ; "None" ; "None" ; "None" ]
 end
 
 (** PutPsiReply_msg : Reply to the PutPsi message *)
@@ -574,6 +545,9 @@ type t =
 | Terminate           of  Terminate_msg.t
 | Ok                  of  Ok_msg.t
 | Error               of  Error_msg.t
+| SetStopped 
+| SetWaiting 
+| SetRunning 
 
 
 let of_string s = 
@@ -606,14 +580,15 @@ let of_string s =
   | "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 )
+                 ~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 )
-  | "ok"         :: [] ->
-       Ok (Ok_msg.create ())
-  | "error"      :: rest ->
-       Error (Error_msg.create (String.concat ~sep:" " rest))
+                ~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"
     
 
@@ -638,6 +613,9 @@ let to_string = function
 | Error               x -> Error_msg.to_string              x
 | PutPsi              x -> PutPsi_msg.to_string             x
 | GetPsiReply         x -> GetPsiReply_msg.to_string        x
+| SetStopped            -> "set_stopped"
+| SetRunning            -> "set_running"
+| SetWaiting            -> "set_waiting"
 
 
 let to_string_list = function

ocaml/Molecule.ml

@@ -147,10 +147,28 @@ let of_xyz_file
   let (_,buffer) = In_channel.read_all filename 
   |> String.lsplit2_exn ~on:'\n' in
   let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in
-  of_xyz_string ~charge:charge ~multiplicity:multiplicity 
-    ~units:units buffer
+  of_xyz_string ~charge ~multiplicity ~units buffer
 
 
+let of_zmt_file  
+    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
+    ?(units=Units.Angstrom)
+    filename = 
+  In_channel.read_all filename
+  |> Zmatrix.of_string
+  |> Zmatrix.to_xyz_string
+  |> of_xyz_string ~charge ~multiplicity ~units
+
+
+let of_file 
+    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
+    ?(units=Units.Angstrom)
+    filename = 
+  try
+    of_xyz_file ~charge ~multiplicity ~units filename
+  with _ ->
+    of_zmt_file ~charge ~multiplicity ~units filename
+
 
 let distance_matrix molecule =
   let coord =

ocaml/Molecule.mli

@@ -29,6 +29,18 @@ val of_xyz_file :
   ?multiplicity:Multiplicity.t ->
   ?units:Units.units -> string -> t
 
+(** Creates a molecule from a zmt file *)
+val of_zmt_file :
+  ?charge:Charge.t -> 
+  ?multiplicity:Multiplicity.t ->
+  ?units:Units.units -> string -> t
+
+(** Creates a molecule from a file (xyz or zmt) *)
+val of_file :
+  ?charge:Charge.t -> 
+  ?multiplicity:Multiplicity.t ->
+  ?units:Units.units -> string -> t
+
 (** Creates a molecule from an xyz file in a string *)
 val of_xyz_string :
   ?charge:Charge.t ->

ocaml/Progress_bar.ml

@@ -14,13 +14,13 @@ type t =
 
 let init ?(bar_length=20) ?(start_value=0.) ?(end_value=1.) ~title =
   { title ; start_value ; end_value ; bar_length ; cur_value=start_value ;
-    init_time= Time.now () ; dirty = true ; next = Time.now () } 
+    init_time= Time.now () ; dirty = false ; next = Time.now () } 
 
 let update ~cur_value bar =
   { bar with cur_value ; dirty=true }
 
 let increment_end bar =
-  { bar with end_value=(bar.end_value +. 1.) ; dirty=true }
+  { bar with end_value=(bar.end_value +. 1.) ; dirty=false }
 
 let increment_cur bar =
   { bar with cur_value=(bar.cur_value +. 1.) ; dirty=true }

ocaml/Qpackage.ml

@@ -127,3 +127,14 @@ let get_ezfio_default directory data =
   |> aux 
 ;;
 
+let ezfio_work ezfio_file = 
+  let result = 
+    Filename.concat ezfio_file  "work"
+  in
+  begin
+    match Sys.is_directory result with
+    | `Yes -> ()
+    | _ -> Unix.mkdir result
+  end;
+  result
+;;

ocaml/Queuing_system.ml

@@ -1,25 +1,35 @@
-open Core.Std
-open Qptypes
-
+module RunningMap = Map.Make (Id.Task)
+module TasksMap   = Map.Make (Id.Task)
+module ClientsSet = Set.Make (Id.Client)
 
 type t =
-{ queued         : Id.Task.t list ;
-  running        : (Id.Task.t, Id.Client.t) Map.Poly.t ;
-  tasks          : (Id.Task.t, string)      Map.Poly.t;
-  clients        : Id.Client.t Set.Poly.t;
+{ queued_front   : Id.Task.t list ;
+  queued_back    : Id.Task.t list ;
+  running        : Id.Client.t RunningMap.t;
+  tasks          : string TasksMap.t;
+  clients        : ClientsSet.t;
   next_client_id : Id.Client.t;
   next_task_id   : Id.Task.t;
+  number_of_queued  : int;
+  number_of_running : int;
+  number_of_tasks   : int;
+  number_of_clients : int;
 }
 
 
 
 let create () =
-  { queued         = [] ; 
-    running        = Map.Poly.empty ;
-    tasks          = Map.Poly.empty;
-    clients        = Set.Poly.empty;
+  { queued_front   = [] ; 
+    queued_back    = [] ; 
+    running        = RunningMap.empty ;
+    tasks          = TasksMap.empty;
+    clients        = ClientsSet.empty;
     next_client_id = Id.Client.of_int 1;
     next_task_id   = Id.Task.of_int 1;
+    number_of_queued  = 0;
+    number_of_running = 0;
+    number_of_tasks   = 0;
+    number_of_clients = 0;
   }
 
 
@@ -30,9 +40,11 @@ let add_task ~task q =
     q.next_task_id 
   in
   { q with
-    queued = task_id :: q.queued ;
-    tasks  = Map.add q.tasks ~key:task_id ~data:task ;
+    queued_front = task_id :: q.queued_front ;
+    tasks  = TasksMap.add task_id task q.tasks;
     next_task_id = Id.Task.increment task_id ;
+    number_of_queued = q.number_of_queued + 1;
+    number_of_tasks  = q.number_of_tasks  + 1;
   }
 
 
@@ -43,55 +55,73 @@ let add_client q =
     q.next_client_id 
   in
   { q with
-    clients = Set.add q.clients client_id;
+    clients = ClientsSet.add client_id q.clients;
     next_client_id = Id.Client.increment client_id;
+    number_of_clients = q.number_of_clients + 1;
   }, client_id
 
 
 let pop_task ~client_id q = 
-  let { queued ; running ; _ } =
+  let { queued_front ; queued_back ; running ; _ } =
     q
   in
-  assert (Set.mem q.clients client_id);
-  match queued with
+  assert (ClientsSet.mem client_id q.clients);
+  let queued_front', queued_back' =
+    match queued_front, queued_back with
+    | (l, []) -> ( [], List.rev l)
+    | t -> t
+  in
+  match queued_back' with
   | task_id :: new_queue ->
     let new_q =
       { q with
-        queued = new_queue ;
-        running = Map.add running ~key:task_id ~data:client_id ;
+        queued_front= queued_front' ;
+        queued_back = new_queue ;
+        running = RunningMap.add task_id client_id running;
+        number_of_queued  = q.number_of_queued  - 1;
+        number_of_running = q.number_of_running + 1;
       }
-    in new_q, Some task_id, (Map.find q.tasks task_id)
+    and found =
+      try Some (TasksMap.find task_id q.tasks)
+      with Not_found -> None
+    in new_q, Some task_id, found 
   | [] -> q, None, None
     
 
 let del_client ~client_id q =
-  assert (Set.mem q.clients client_id);
+  assert (ClientsSet.mem client_id q.clients);
   { q with 
-    clients = Set.remove q.clients client_id }
+    clients = ClientsSet.remove client_id q.clients;
+    number_of_clients = q.number_of_clients - 1
+  }
 
 
 let end_task ~task_id ~client_id q = 
   let { running ; tasks ; _ } =
     q
   in
-  assert (Set.mem q.clients client_id);
-  let () = 
-    match Map.Poly.find running task_id with
-    | None -> failwith "Task already finished"
-    | Some client_id_check -> assert (client_id_check = client_id)
+  assert (ClientsSet.mem client_id q.clients);
+  let () =
+    let client_id_check = 
+      try RunningMap.find task_id running with
+      Not_found -> failwith "Task already finished"
+    in
+    assert (client_id_check = client_id)
   in
   { q with
-    running  = Map.remove running task_id ;
+    running = RunningMap.remove task_id running ;
+    number_of_running = q.number_of_running - 1
   }
-    
+
 let del_task ~task_id q = 
   let { tasks ; _ } =
     q
   in
   
-  if (Map.mem tasks task_id) then
+  if (TasksMap.mem task_id tasks) then
       { q with
-        tasks    = Map.remove tasks task_id ;
+        tasks = TasksMap.remove task_id tasks;
+        number_of_tasks = q.number_of_tasks - 1;
       }
   else
       Printf.sprintf "Task %d is already deleted" (Id.Task.to_int task_id)
@@ -99,33 +129,81 @@ let del_task ~task_id q =
 
     
 
+let number_of_tasks q =
+  assert (q.number_of_tasks >= 0);
+  q.number_of_tasks
+
 let number_of_queued q =
-  Map.length q.tasks
+  assert (q.number_of_queued >= 0);
+  q.number_of_queued
 
 let number_of_running q =
-  Map.length q.running
+  assert (q.number_of_running >= 0);
+  q.number_of_running
+
+let number_of_clients q =
+  assert (q.number_of_clients >= 0);
+  q.number_of_clients
 
 
-let to_string { queued ; running ; tasks ; _ } =
+let to_string qs =
+  let { queued_back ; queued_front ; running ; tasks ; _ } = qs in
   let q =
-    List.map ~f:Id.Task.to_string queued
-    |> String.concat ~sep:" ; "
+     (List.map Id.Task.to_string queued_front) @
+     (List.map Id.Task.to_string @@ List.rev queued_back)
+    |> String.concat " ; "
   and r =
-    Map.Poly.to_alist running
-    |> List.map ~f:(fun (t,c) -> "("^(Id.Task.to_string t)^", "
+    RunningMap.bindings running
+    |> List.map (fun (t,c) -> "("^(Id.Task.to_string t)^", "
        ^(Id.Client.to_string c)^")") 
-    |> String.concat ~sep:" ; "
+    |> String.concat " ; "
   and t = 
-    Map.Poly.to_alist tasks
-    |> List.map ~f:(fun (t,c) -> "("^(Id.Task.to_string t)^", \""
+    TasksMap.bindings tasks
+    |> List.map (fun (t,c) -> "("^(Id.Task.to_string t)^", \""
        ^c^"\")") 
-    |> String.concat ~sep:" ; "
+    |> String.concat " ; "
   in
   Printf.sprintf "{
+Tasks : %d   Queued : %d   Running : %d   Clients : %d
 queued   : { %s }
 running  : { %s }
 tasks    : [ %s
            ]
-}" q r t
+}"
+(number_of_tasks qs) (number_of_queued qs) (number_of_running qs) (number_of_clients qs)
+q r t
 
   
+
+let test () =
+  let q = 
+    create ()
+    |> add_task ~task:"First Task"
+    |> add_task ~task:"Second Task"
+  in
+  let q, client_id = 
+    add_client q
+  in
+  let q, task_id, task_content =
+    match pop_task ~client_id q with
+    | q, Some x, Some y -> q, Id.Task.to_int x, y
+    | _ -> assert false
+  in
+  Printf.printf "Task_id : %d \t\t Task : %s\n" task_id task_content;
+  to_string q |> print_endline  ;
+  let q, task_id, task_content =
+    match pop_task ~client_id q with
+    | q, Some x, Some y -> q, Id.Task.to_int x, y
+    | _ -> assert false
+  in
+  Printf.printf "Task_id : %d \t\t Task : %s\n" task_id task_content;
+  let q, task_id, task_content =
+    match pop_task ~client_id q with
+    | q, None, None -> q, 0, "None"
+    | _ -> assert false
+  in
+  Printf.printf "Task_id : %d \t\t Task : %s\n"  task_id task_content;
+  q
+  |> to_string
+  |> print_endline 
+

ocaml/Queuing_system.mli

@@ -0,0 +1,63 @@
+module RunningMap : Map.S with type key = Id.Task.t
+module TasksMap   : Map.S with type key = Id.Task.t
+module ClientsSet : Set.S with type elt = Id.Client.t
+
+type t = {
+  queued_front      : Id.Task.t list           ;
+  queued_back       : Id.Task.t list           ;
+  running           : Id.Client.t RunningMap.t ;
+  tasks             : string TasksMap.t        ;          
+  clients           : ClientsSet.t             ;          
+  next_client_id    : Id.Client.t              ;
+  next_task_id      : Id.Task.t                ;
+  number_of_queued  : int                      ;
+  number_of_running : int                      ;
+  number_of_tasks   : int                      ;
+  number_of_clients : int                      ;
+}
+
+(** Creates a new queuing system. Returns the new queue. *)
+val create : unit -> t
+
+(** Add a new task represented as a string. Returns the queue with the added task. *)
+val add_task : task:string -> t -> t
+
+(** Add a new client. Returns the queue and a new client_id. *)
+val add_client : t -> t * Id.Client.t
+
+(** Pops a task from the queue. The task is set as running on client client_id.
+    Returns the queue, a task_id and the content of the task. If the queue contains
+    no task, the task_id and the task content are None. *)
+val pop_task :
+  client_id:ClientsSet.elt -> t -> t * Id.Task.t option * string option
+
+(** Deletes a client from the queuing system *)
+val del_client : client_id:ClientsSet.elt -> t -> t
+
+(** Deletes a client from the queuing system. The client is assumed to be a member
+    of the set of clients. Returns the queue without the removed client. *)
+val end_task : task_id:RunningMap.key -> client_id:ClientsSet.elt -> t -> t
+
+(** Deletes a task from the queuing system. The task is assumed to be a member
+    of the map of tasks. Returns the queue without the removed task. *)
+val del_task : task_id:TasksMap.key -> t -> t
+
+(** Returns the number of tasks, assumed >= 0 *)
+val number_of_tasks : t -> int
+
+(** Returns the number of queued tasks, assumed >= 0 *)
+val number_of_queued : t -> int
+
+(** Returns the number of running tasks, assumed >= 0 *)
+val number_of_running : t -> int
+
+(** Returns the number of connected clients, assumed >= 0 *)
+val number_of_clients : t -> int
+
+(** Prints the content of the queue *)
+val to_string : t -> string
+
+(** Test function for debug *)
+val test : unit -> unit
+
+

ocaml/TaskServer.ml

@@ -2,6 +2,23 @@ open Core.Std
 open Qptypes
 
 
+type pub_state =
+| Waiting
+| Running of string
+| Stopped
+
+let pub_state_of_string = function
+| "Waiting"  -> Waiting
+| "Stopped"  -> Stopped
+| s -> Running s
+
+let string_of_pub_state = function
+| Waiting  -> "Waiting"
+| Stopped  -> "Stopped"
+| Running s -> s
+
+
+
 type t =
 {
     queue           : Queuing_system.t ;
@@ -31,20 +48,21 @@ let zmq_context =
   ZMQ.Context.create ()
 
 
-let bind_socket ~socket_type ~socket ~address =
+let bind_socket ~socket_type ~socket ~port =
   let rec loop = function
   | 0 -> failwith @@ Printf.sprintf
-        "Unable to bind the %s socket : %s "
-        socket_type address
+        "Unable to bind the %s socket to port : %d "
+        socket_type port
   | -1 -> ()
   | i -> 
       try
-        ZMQ.Socket.bind socket address;
+        ZMQ.Socket.bind socket @@ Printf.sprintf "tcp://*:%d" port;
         loop (-1)
       with
       | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) )
       | other_exception -> raise other_exception
-  in loop 10
+  in loop 60;
+  ZMQ.Socket.bind socket @@ Printf.sprintf "ipc:///tmp/qp_run:%d" port
 
 
 let hostname = lazy (
@@ -98,7 +116,7 @@ let stop ~port =
     let req_socket =
       ZMQ.Socket.create zmq_context ZMQ.Socket.req
     and address =
-      Printf.sprintf "tcp://%s:%d" (Lazy.force ip_address) port
+      Printf.sprintf "ipc:///tmp/qp_run:%d" port
     in
     ZMQ.Socket.set_linger_period req_socket 1_000_000;
     ZMQ.Socket.connect req_socket address;
@@ -120,7 +138,7 @@ let stop ~port =
     ZMQ.Socket.close req_socket
 
 
-let new_job msg program_state rep_socket =
+let new_job msg program_state rep_socket pair_socket =
 
     let state = 
         msg.Message.Newjob_msg.state
@@ -143,10 +161,32 @@ let new_job msg program_state rep_socket =
       }
     in
     reply_ok rep_socket;
+    string_of_pub_state Waiting
+    |> ZMQ.Socket.send pair_socket ;
     result
 
+let change_pub_state msg program_state rep_socket pair_socket =
+  let msg = 
+    match msg with
+    | `Waiting -> Waiting 
+    | `Stopped -> Stopped
+    | `Running ->
+      begin
+	let state =
+	   match program_state.state with
+	   | Some x -> x
+	   | None  -> failwith "Trying to change pub state while no job is ready"
+        in
+        Running (Message.State.to_string state)
+      end
+  in
+  reply_ok rep_socket;
+  string_of_pub_state msg
+  |> ZMQ.Socket.send pair_socket ;
 
-let end_job msg program_state rep_socket =
+  program_state
+
+let end_job msg program_state rep_socket pair_socket =
 
     let failure () =
         reply_wrong_state rep_socket;
@@ -165,7 +205,11 @@ let end_job msg program_state rep_socket =
     | Some state -> 
       begin
         if (msg.Message.Endjob_msg.state = state) then
-          success state
+          begin
+            string_of_pub_state Waiting
+            |> ZMQ.Socket.send pair_socket ;
+            success state
+          end
         else
           failure ()
       end
@@ -262,8 +306,7 @@ let del_task msg program_state rep_socket =
             }
         in
         let more = 
-            (Queuing_system.number_of_queued new_program_state.queue +
-             Queuing_system.number_of_running new_program_state.queue) > 0
+            (Queuing_system.number_of_tasks new_program_state.queue > 0)
         in
         Message.DelTaskReply (Message.DelTaskReply_msg.create ~task_id ~more)
         |> Message.to_string
@@ -355,7 +398,7 @@ let add_task msg program_state rep_socket =
 
 
 
-let get_task msg program_state rep_socket =
+let get_task msg program_state rep_socket pair_socket =
 
     let state, client_id =
       msg.Message.GetTask_msg.state, 
@@ -371,6 +414,12 @@ let get_task msg program_state rep_socket =
         let new_queue, task_id, task =
             Queuing_system.pop_task ~client_id program_state.queue
         in
+        if (Queuing_system.number_of_queued new_queue = 0) then
+          string_of_pub_state Waiting 
+          |> ZMQ.Socket.send pair_socket
+        else
+          string_of_pub_state (Running (Message.State.to_string state))
+          |> ZMQ.Socket.send pair_socket;
 
         let new_program_state = 
             { program_state with
@@ -378,21 +427,10 @@ let get_task msg program_state rep_socket =
             }
         in
 
-        match (task, task_id) with
-          | Some task, Some task_id -> 
-              begin
-                Message.GetTaskReply (Message.GetTaskReply_msg.create ~task ~task_id)
-                |> Message.to_string
-                |> ZMQ.Socket.send rep_socket ;
-                new_program_state
-              end
-          | _ -> 
-              begin
-                Message.Terminate (Message.Terminate_msg.create ())
-                |> Message.to_string
-                |> ZMQ.Socket.send rep_socket ;
-                program_state
-              end
+        Message.GetTaskReply (Message.GetTaskReply_msg.create ~task ~task_id)
+        |> Message.to_string
+        |> ZMQ.Socket.send rep_socket ;
+        new_program_state
 
     in
 
@@ -454,9 +492,9 @@ let put_psi msg rest_of_msg program_state rep_socket =
         | Some x -> x
         | None ->
           begin
-            let psi_det, psi_coef =
+            let psi_det, psi_coef, energy =
               match rest_of_msg with
-              | [ x ; y ] -> x, y
+              | [ x ; y ; e ] -> x, y, e
               | _ -> failwith "Badly formed put_psi message"
             in
             Message.Psi.create
@@ -467,6 +505,7 @@ let put_psi msg rest_of_msg program_state rep_socket =
               ~n_det_selectors:msg.Message.PutPsi_msg.n_det_selectors
               ~psi_det
               ~psi_coef
+              ~energy
           end
     in
     let new_program_state =
@@ -501,29 +540,85 @@ let get_psi msg program_state rep_socket =
 let terminate program_state rep_socket =
     reply_ok rep_socket;
     { program_state with
+      psi = None;
+      address_tcp = None;
+      address_inproc = None;
       running = false
     }
 
   
 let error msg program_state rep_socket =
-    Printf.printf "%s\n%!" msg;
     Message.Error (Message.Error_msg.create msg)
     |> Message.to_string
     |> ZMQ.Socket.send rep_socket ;
     program_state
 
+let start_pub_thread ~port = 
+    Thread.create (fun () ->
+      let timeout =
+        1000
+      in
 
+      let pair_socket =
+        ZMQ.Socket.create zmq_context ZMQ.Socket.pair
+      and address = 
+        "inproc://pair" 
+      in
+      ZMQ.Socket.connect pair_socket address;
+
+      let pub_socket =
+        ZMQ.Socket.create zmq_context ZMQ.Socket.pub
+      in
+      bind_socket ~socket_type:"PUB" ~socket:pub_socket ~port;
+
+      let pollitem =
+        ZMQ.Poll.mask_of 
+          [| (pair_socket, ZMQ.Poll.In) |]
+      in
+
+      let rec run state = 
+        let new_state =
+          let polling = 
+            ZMQ.Poll.poll ~timeout pollitem
+          in
+          if (polling.(0) = Some ZMQ.Poll.In) then
+            ZMQ.Socket.recv ~block:false pair_socket
+            |> pub_state_of_string
+          else
+            state
+        in
+        ZMQ.Socket.send pub_socket @@ string_of_pub_state new_state;
+        match state with
+        | Stopped -> ()
+        | _ -> run new_state
+      in
+      run Waiting;
+      ZMQ.Socket.set_linger_period pair_socket 1000 ;
+      ZMQ.Socket.close pair_socket;
+      ZMQ.Socket.set_linger_period pub_socket 1000 ;
+      ZMQ.Socket.close pub_socket;
+    )
 
 let run ~port =
 
+    (** Bind inproc socket for changing state of pub *)
+    let pair_socket =
+      ZMQ.Socket.create zmq_context ZMQ.Socket.pair
+    and address = 
+      "inproc://pair" 
+    in
+    ZMQ.Socket.bind pair_socket address;
+
+    let pub_thread = 
+      start_pub_thread ~port:(port+1) ()
+    in
+
     (** Bind REP socket *)
     let rep_socket =
       ZMQ.Socket.create zmq_context ZMQ.Socket.rep
-    and address =
-      Printf.sprintf "tcp://%s:%d" (Lazy.force ip_address) port
     in
-    bind_socket "REP" rep_socket address;
     ZMQ.Socket.set_linger_period rep_socket 1_000_000;
+    bind_socket "REP" rep_socket port;
 
     let initial_program_state =
     {   queue = Queuing_system.create () ;
@@ -542,6 +637,9 @@ let run ~port =
         [| (rep_socket, ZMQ.Poll.In) |]
     in
 
+    let address =
+      Printf.sprintf "tcp://%s:%d" (Lazy.force ip_address) port
+    in
     Printf.printf "Task server running : %s\n%!" address;
 
 
@@ -579,9 +677,10 @@ let run ~port =
               in
 
               (** Debug input *)
-              Printf.sprintf "%d %d : %s\n%!"
-              (Queuing_system.number_of_queued program_state.queue)
+              Printf.sprintf "q:%d  r:%d  n:%d  : %s\n%!"
+              (Queuing_system.number_of_queued  program_state.queue)
               (Queuing_system.number_of_running program_state.queue)
+              (Queuing_system.number_of_tasks   program_state.queue)
               (Message.to_string message)
               |> debug;
 
@@ -591,15 +690,18 @@ let run ~port =
                   | _     , Message.Terminate   _ -> terminate program_state rep_socket
                   | _     , Message.PutPsi      x -> put_psi x rest program_state rep_socket
                   | _     , Message.GetPsi      x -> get_psi x program_state rep_socket
-                  | None  , Message.Newjob      x -> new_job x program_state rep_socket
+                  | None  , Message.Newjob      x -> new_job x program_state rep_socket pair_socket
                   | _     , Message.Newjob      _ -> error "A job is already running" program_state rep_socket
-                  | Some _, Message.Endjob      x -> end_job x program_state rep_socket
+                  | Some _, Message.Endjob      x -> end_job x program_state rep_socket pair_socket
+                  | Some _, Message.SetRunning    -> change_pub_state `Running program_state rep_socket pair_socket
+                  | _, Message.SetWaiting    -> change_pub_state `Waiting program_state rep_socket pair_socket
+                  | _, Message.SetStopped    -> change_pub_state `Stopped program_state rep_socket pair_socket
                   | None  ,                     _ -> error "No job is running" program_state rep_socket
                   | Some _, Message.Connect     x -> connect x program_state rep_socket
                   | Some _, Message.Disconnect  x -> disconnect x program_state rep_socket
                   | Some _, Message.AddTask     x -> add_task x program_state rep_socket
                   | Some _, Message.DelTask     x -> del_task x program_state rep_socket
-                  | Some _, Message.GetTask     x -> get_task x program_state rep_socket
+                  | Some _, Message.GetTask     x -> get_task x program_state rep_socket pair_socket
                   | Some _, Message.TaskDone    x -> task_done x program_state rep_socket
                   | _     , _                     ->
                     error ("Invalid message : "^(Message.to_string message))  program_state rep_socket
@@ -614,6 +716,11 @@ let run ~port =
           end
     in main_loop initial_program_state true;
 
+    ZMQ.Socket.send pair_socket @@ string_of_pub_state Stopped;
+    Thread.join pub_thread;
+    ZMQ.Socket.close rep_socket
+
+
 
 
     

ocaml/TaskServer.mli

@@ -23,9 +23,9 @@ val debug : string -> unit
 (** ZeroMQ context *)
 val zmq_context : ZMQ.Context.t
 
-(** Bind a ZMQ socket *)
+(** Bind a ZMQ socket to a TCP port and to an IPC file /tmp/qp_run.<port> *)
 val bind_socket :
-  socket_type:string -> socket:'a ZMQ.Socket.t -> address:string -> unit
+  socket_type:string -> socket:'a ZMQ.Socket.t -> port:int -> unit
 
 (** Name of the host on which the server runs *)
 val hostname : string lazy_t
@@ -43,10 +43,10 @@ val stop : port:int -> unit
 (** {1} Server functions *)
 
 (** Create a new job *)
-val new_job : Message.Newjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+val new_job : Message.Newjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> [> `Pair] ZMQ.Socket.t -> t
 
 (** Finish a running job *)
-val end_job : Message.Endjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+val end_job : Message.Endjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> [> `Pair] ZMQ.Socket.t -> t
 
 (** Connect a client *)
 val connect: Message.Connect_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
@@ -64,7 +64,7 @@ val task_done: Message.TaskDone_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
 val del_task: Message.DelTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
 
 (** The client get a new task to execute *)
-val get_task: Message.GetTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+val get_task: Message.GetTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t ->  [> `Pair] ZMQ.Socket.t -> t
 
 (** Terminate server *)
 val terminate : t -> [> `Req ] ZMQ.Socket.t -> t

ocaml/Zmatrix.ml

@@ -0,0 +1,326 @@
+open Qptypes
+
+module StringMap = Map.Make(String)
+
+type atom_id  = int
+type angle    = Label of string | Value of float
+type distance = Label of string | Value of float
+type dihedral = Label of string | Value of float
+
+let pi = acos (-1.)
+let to_radian = pi /. 180.
+
+let rec in_range (xmin, xmax) x =
+   if (x <= xmin) then
+     in_range (xmin, xmax) (x -. xmin +. xmax )
+   else if (x > xmax) then
+     in_range (xmin, xmax) (x -. xmax +. xmin )
+   else
+     x
+
+let atom_id_of_int : int -> atom_id = 
+  fun x -> ( assert (x>0) ; x)
+
+let distance_of_float : float -> distance = 
+  fun x -> ( assert (x>=0.) ; Value x)
+
+let angle_of_float : float -> angle = 
+  fun x -> Value (in_range (-180., 180.) x)
+
+let dihedral_of_float : float -> dihedral = 
+  fun x -> Value (in_range (-360., 360.) x)
+
+
+let atom_id_of_string : string -> atom_id = 
+  fun i -> atom_id_of_int @@ int_of_string  i
+
+let distance_of_string : string -> distance =
+  fun s -> 
+    try
+      distance_of_float @@ float_of_string s
+    with _ -> Label s
+
+let angle_of_string : string -> angle =
+  fun s -> 
+    try
+      angle_of_float @@ float_of_string s
+    with _ -> Label s
+
+let dihedral_of_string : string -> dihedral =
+  fun s -> 
+    try
+      dihedral_of_float @@ float_of_string s
+    with _ -> Label s
+
+
+let int_of_atom_id : atom_id -> int = fun x -> x
+
+let float_of_distance : float StringMap.t -> distance -> float = 
+  fun map -> function
+  | Value x -> x
+  | Label s -> StringMap.find s map
+
+let float_of_angle : float StringMap.t -> angle -> float =
+  fun map -> function
+  | Value x -> x
+  | Label s -> StringMap.find s map
+
+let float_of_dihedral : float StringMap.t -> dihedral -> float =
+  fun map -> function
+  | Value x -> x
+  | Label s -> StringMap.find s map
+
+
+type line = 
+| First  of  Element.t
+| Second of (Element.t * distance)
+| Third  of (Element.t * atom_id * distance * atom_id * angle)
+| Other  of (Element.t * atom_id * distance * atom_id * angle * atom_id * dihedral )
+| Coord  of (string * float)
+
+
+let string_of_line map = 
+  let f_r = float_of_distance map
+  and f_a = float_of_angle    map
+  and f_d = float_of_dihedral map
+  and i_i = int_of_atom_id
+  in function
+| First  e ->  Printf.sprintf "%-3s" (Element.to_string e)
+| Second (e, r) -> Printf.sprintf "%-3s %5d %f" (Element.to_string e) 1 (f_r r)
+| Third  (e, i, r, j, a) -> Printf.sprintf "%-3s %5d %f %5d %f" (Element.to_string e) (i_i i) (f_r r) (i_i j) (f_a a)
+| Other  (e, i, r, j, a, k, d) -> Printf.sprintf "%-3s %5d %f %5d %f %5d %f" (Element.to_string e) (i_i i) (f_r r) (i_i j) (f_a a) (i_i k) (f_d d)
+| Coord  (c, f) -> Printf.sprintf "%s  %f" c f
+
+
+let line_of_string l =
+  let line_clean =
+    Str.split (Str.regexp " ") l
+    |> List.filter (fun x -> x <> "")
+  in
+  match line_clean with
+  | e :: [] -> First (Element.of_string e)
+  | e :: i :: r :: [] -> Second
+    (Element.of_string e,
+     distance_of_string r)
+  | e :: i :: r :: j :: a :: [] -> Third 
+    (Element.of_string e,
+     atom_id_of_string i,
+     distance_of_string r,
+     atom_id_of_string j,
+     angle_of_string a)
+  | e :: i :: r :: j :: a :: k :: d :: [] -> Other
+    (Element.of_string e,
+     atom_id_of_string i,
+     distance_of_string r,
+     atom_id_of_string j,
+     angle_of_string a,
+     atom_id_of_string k,
+     dihedral_of_string d)
+  | c :: f :: [] -> Coord (c, float_of_string f)
+  | _ -> failwith ("Syntax error: "^l)
+
+
+type t = (line array * float StringMap.t)
+
+let of_string t =
+  let l =
+    Str.split (Str.regexp "\n") t
+    |> List.map String.trim
+    |> List.filter (fun x -> x <> "")
+    |> List.map line_of_string
+  in
+
+  let l = 
+    match l with
+    | First _ :: Second _ :: Third _ :: _
+    | First _ :: Second _ :: Coord _ :: []
+    | First _ :: Second _ :: []
+    | First _ :: [] -> l
+    | _ -> failwith "Syntax error"
+  in
+
+  let (l, m) =
+    let rec work lst map = function
+    | (First  _ as x) :: rest
+    | (Second _ as x) :: rest
+    | (Third  _ as x) :: rest
+    | (Other  _ as x) :: rest -> work (x::lst) map rest
+    | (Coord  (c,f)) :: rest -> work lst (StringMap.add c f map) rest
+    | [] -> (List.rev lst, map)
+    in
+    work [] (StringMap.empty) l
+  in
+  (Array.of_list l, m)
+ 
+
+(** Linear algebra *)
+
+let (|-) (x,y,z) (x',y',z') =
+  ( x-.x', y-.y', z-.z' )
+
+let (|+) (x,y,z) (x',y',z') =
+  ( x+.x', y+.y', z+.z' )
+
+let (|.) s (x,y,z) =
+  ( s*.x, s*.y, s*.z )
+
+let dot (x,y,z) (x',y',z') =
+  x*.x' +. y*.y' +. z*.z'
+
+let norm u =
+  sqrt @@ dot u u
+
+let normalized u =
+  1. /. (norm u) |. u
+
+let cross (x,y,z) (x',y',z') =
+  ((y *. z' -. z *. y'), -. (x *. z' -. z *. x'), (x *. y' -. y *. x'))
+
+let rotation_matrix axis angle = 
+   (* Euler-Rodrigues formula for rotation matrix, taken from
+      https://github.com/jevandezande/zmatrix/blob/master/converter.py
+   *)
+   let a = 
+      (cos (angle *. to_radian *. 0.5))
+   in
+   let (b, c, d) = 
+      (-. sin (angle *. to_radian *. 0.5)) |. (normalized axis)
+   in
+   Array.of_list @@ 
+     [(a *. a +. b *. b -. c *. c -. d *. d,
+       2. *. (b *. c -. a *. d),
+       2. *. (b *. d +. a *. c));
+      (2. *. (b *. c +. a *. d),
+       a *. a +. c *. c -.b *. b -. d *. d,
+       2. *. (c *. d -. a *. b));
+      (2. *. (b *. d -. a *. c),
+       2. *. (c *. d +. a *. b),
+       a *. a +. d *. d -. b *. b -. c *. c)]
+(*
+     [(a *. a +. b *. b -. c *. c -. d *. d,
+       2. *. (b *. c +. a *. d),
+       2. *. (b *. d -. a *. c));
+      (2. *. (b *. c -. a *. d),
+       a *. a +. c *. c -.b *. b -. d *. d,
+       2. *. (c *. d +. a *. b));
+      (2. *. (b *. d +. a *. c),
+       2. *. (c *. d -. a *. b),
+       a *. a +. d *. d -. b *. b -. c *. c)]
+*)
+      
+
+
+let apply_rotation_matrix rot u =
+  (dot rot.(0) u, dot rot.(1) u, dot rot.(2) u)
+  
+let center_of_mass l =
+let (x,y,z) =
+  let sum_mass, com =
+    Array.fold_left (fun (s,com) (e,x,y,z) ->
+      let mass =
+         Positive_float.to_float @@ Element.mass e
+      in
+      (s +. mass, ( mass |. (x,y,z) ) |+ com) ) 
+      (0., (0.,0.,0.)) l
+  in
+  (1. /. sum_mass) |. com
+in
+Printf.printf "%f %f %f\n" x y z ; (x,y,z)
+
+let to_xyz (z,map) =
+  let result =
+    Array.make (Array.length z) None
+  in 
+
+  let get_cartesian_coord i =
+    match result.(i-1) with
+    | None -> failwith @@ Printf.sprintf "Atom %d is defined in the future" i
+    | Some (_, x, y, z) -> (x, y, z)
+  in
+
+
+  let append_line i' =
+    match z.(i') with
+    | First e -> 
+        result.(i') <- Some (e, 0., 0., 0.)
+    | Second (e, r) -> 
+        let r =
+          float_of_distance map r
+        in
+        result.(i') <- Some (e, 0., 0., r)
+    | Third  (e, i, r, j, a) -> 
+      begin
+        let i, r, j, a =
+          int_of_atom_id i,
+          float_of_distance map r,
+          int_of_atom_id j,
+          float_of_angle map a
+        in
+        let ui, uj =
+          get_cartesian_coord i, 
+          get_cartesian_coord j
+        in
+        let u_ij = 
+          (uj |- ui)
+        in
+        let rot = 
+          rotation_matrix (0., 1., 0.) a
+        in
+        let new_vec =
+          apply_rotation_matrix rot ( r |. (normalized u_ij))
+        in
+        let (x, y, z) =
+          new_vec |+ ui
+        in
+        result.(i') <- Some (e, x, y, z)
+      end
+    | Other  (e, i, r, j, a, k, d) -> 
+      begin
+        let i, r, j, a, k, d =
+          int_of_atom_id i,
+          float_of_distance map r,
+          int_of_atom_id j,
+          float_of_angle map a,
+          int_of_atom_id k,
+          float_of_dihedral map d
+        in
+        let ui, uj, uk =
+          get_cartesian_coord i, 
+          get_cartesian_coord j,
+          get_cartesian_coord k
+        in
+        let u_ij, u_kj =
+          (uj |- ui) , (uj |- uk)
+        in
+        let normal =
+          cross u_ij u_kj
+        in
+        let new_vec = 
+          r |. (normalized u_ij)
+          |> apply_rotation_matrix (rotation_matrix normal a)
+          |> apply_rotation_matrix (rotation_matrix u_ij d)
+        in
+        let (x, y, z) =
+          new_vec |+ ui
+        in
+        result.(i') <- Some (e, x, y, z)
+      end
+    | Coord _ -> ()
+  in
+  Array.iteri (fun i _ -> append_line i) z;
+  let result = 
+    Array.map (function
+    | Some x -> x
+    | None -> failwith "Some atoms were not defined" ) result
+  in
+  Array.to_list result
+
+
+let to_xyz_string (l,map) =
+  String.concat "\n" 
+    ( to_xyz (l,map) 
+      |> List.map (fun (e,x,y,z) -> 
+        Printf.sprintf "%s %f %f %f\n" (Element.to_string e) x y z)  )
+
+
+

ocaml/_tags

@@ -1,3 +1,3 @@
-true: package(core,cryptokit,ZMQ,sexplib.syntax)
+true: package(core,cryptokit,ZMQ,sexplib.syntax,str)
 true: thread
 false: profile

ocaml/qp_create_ezfio_from_xyz.ml

@@ -19,7 +19,7 @@ let spec =
      ~doc:"string Name of the pseudopotential"
   +> flag "cart" no_arg
      ~doc:" Compute AOs in the Cartesian basis set (6d, 10f, ...)"
-  +> anon ("xyz_file" %: file )
+  +> anon ("(xyz_file|zmt_file)" %: file )
 
 
 (** Handle dummy atoms placed on bonds *)
@@ -93,7 +93,7 @@ let run ?o b c d m p cart xyz_file =
 
   (* Read molecule *)
   let molecule =
-    (Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c)
+    (Molecule.of_file xyz_file ~charge:(Charge.of_int c)
       ~multiplicity:(Multiplicity.of_int m) )
   in
   let dummy =
@@ -309,7 +309,8 @@ let run ?o b c d m p cart xyz_file =
     | None ->
       begin
         match String.rsplit2 ~on:'.' xyz_file with
-        | Some (x,"xyz") -> x^".ezfio"
+        | Some (x,"xyz") 
+        | Some (x,"zmt") -> x^".ezfio"
         | _ -> xyz_file^".ezfio"
       end
   in
@@ -640,9 +641,10 @@ let command =
 
 ============================
 
-Creates an EZFIO directory from a standard xyz file.  The basis set is defined
-as a single string if all the atoms are taken from the same basis set,
-otherwise specific elements can be defined as follows:
+Creates an EZFIO directory from a standard xyz file or from a z-matrix file
+in Gaussian format. The basis set is defined as a single string if all the
+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\"
 

ocaml/qp_create_guess.ml

@@ -0,0 +1,141 @@
+open Qputils
+open Qptypes
+open Core.Std
+
+let run ~multiplicity ezfio_file =
+  if (not (Sys.file_exists_exn ezfio_file)) then
+    failwith ("EZFIO directory "^ezfio_file^" not found");
+  Ezfio.set_file ezfio_file;
+  let d =
+    Input.Determinants_by_hand.read ()
+  in
+  let m =
+    Multiplicity.of_int multiplicity
+  in
+  let ne = 
+    Ezfio.get_electrons_elec_alpha_num () +
+    Ezfio.get_electrons_elec_beta_num () 
+    |> Elec_number.of_int
+  in
+  let alpha, beta = 
+    let (a,b) = 
+      Multiplicity.to_alpha_beta ne m 
+    in
+    (Elec_alpha_number.to_int a, Elec_beta_number.to_int b)
+  in
+  let n_open_shells =
+    alpha - beta
+  in
+  let mo_tot_num = 
+    Ezfio.get_mo_basis_mo_tot_num ()
+  in
+  let build_list_of_dets ne n_closed n_open =
+    let init = 
+      Array.create ~len:n_closed Bit.One 
+      |> Array.to_list
+    in
+    let rec set_electron accu = function
+    | 1 ->  [ Bit.One :: accu ]
+    | i -> 
+          assert (i>1);
+          let rest =
+            set_electron (Bit.Zero :: accu) (i-1)
+          in
+          (Bit.One::accu) :: rest
+    in
+    let rec extend accu = function 
+    | 0 -> List.rev accu
+    | i -> extend (Bit.Zero::accu) (i-1)
+    in
+    let rec set_n_electrons accu imax = function
+    | 0 -> []
+    | 1 -> set_electron accu imax
+    | i ->
+          assert (i>1);
+          let l =
+            set_electron accu (imax-1)
+          in
+          List.map ~f:(fun x -> set_n_electrons x (imax-1) (i-1)) l
+          |> List.concat
+    in
+    set_n_electrons init n_open ne
+          |> List.filter ~f:(fun x -> List.length x <= n_closed+n_open) 
+          |> List.map ~f:(fun x -> extend x (((mo_tot_num-1)/64+1)*64 - List.length x))
+  in
+
+  let alpha_new = 
+    (Elec_number.to_int ne + 1)/2
+  and beta_new = 
+    Elec_number.to_int ne/2
+  in
+  let l_alpha = 
+    build_list_of_dets ((alpha-beta+1)/2) beta n_open_shells
+  in 
+  let l_beta = 
+    if alpha_new = beta_new then
+      l_alpha
+    else
+      build_list_of_dets ((alpha-beta)/2)beta n_open_shells 
+  in
+
+  let n_int = 
+    Bitlist.n_int_of_mo_tot_num mo_tot_num
+  in
+  let determinants = 
+    List.map l_alpha ~f:(fun x -> List.map l_beta ~f:(fun y -> (x,y) ))
+    |> List.concat
+    |> List.map ~f:(fun pair -> Determinant.of_bitlist_couple ~n_int
+      ~alpha:(Elec_alpha_number.of_int alpha_new)
+      ~beta:(Elec_beta_number.of_int beta_new) pair )
+  in
+  let c = 
+    Array.create ~len:(List.length determinants) (Det_coef.of_float 1.)
+  in
+ 
+  determinants 
+  |> List.map ~f:(fun x -> Determinant.to_string ~mo_tot_num:(MO_number.of_int mo_tot_num) x) 
+  |> List.iter ~f:(fun x -> Printf.printf "%s\n\n%!" x);
+
+  let l =
+    List.length determinants
+  in
+  if l > 0 then
+    begin
+      let d = 
+        let s = (Float.of_int (alpha - beta)) *. 0.5 in
+        let open Input.Determinants_by_hand in
+        { d with n_int ; 
+          n_det = Det_number.of_int ~min:1 ~max:l  l;
+          expected_s2 = Positive_float.of_float (s *. (s +. 1.)) ;
+          psi_coef = c;
+          psi_det = Array.of_list determinants;
+        }
+      in
+      Input.Determinants_by_hand.write d;
+      Ezfio.set_determinants_read_wf true
+    end
+  else
+    Ezfio.set_determinants_read_wf false
+  
+
+
+let spec =
+  let open Command.Spec in
+  empty
+  +> flag "m" (required int)
+     ~doc:"int Spin multiplicity"
+  +> anon ("ezfio_file" %: string)
+
+let () = 
+  Command.basic
+  ~summary: "Quantum Package command"
+  ~readme:( fun () -> "
+Creates an open-shell multiplet initial guess\n\n" )
+  spec
+  (fun multiplicity ezfio_file () ->
+    run ~multiplicity ezfio_file
+  )
+  |> Command.run   ~version: Git.sha1   ~build_info: Git.message
+
+
+

ocaml/qp_overlap_of_wf.ml

@@ -0,0 +1,66 @@
+open Input_determinants_by_hand
+open Qptypes
+
+let () =
+  let ezfio, ezfio' =
+    try
+      Sys.argv.(1), Sys.argv.(2)
+    with Invalid_argument _ ->
+      raise (Invalid_argument (Printf.sprintf
+         "Syntax : %s EZFIO1 EZFIO2" Sys.argv.(0)))
+  in
+
+  let fetch_wf filename =
+    Ezfio.set_file filename;
+    let mo_tot_num =
+      Ezfio.get_mo_basis_mo_tot_num ()
+      |> MO_number.of_int
+    in
+    let d =
+      Determinants_by_hand.read ()
+    in
+    let n_det =
+      Det_number.to_int d.Determinants_by_hand.n_det
+    in
+    let keys = 
+      Array.map (Determinant.to_string ~mo_tot_num) 
+        d.Determinants_by_hand.psi_det
+    and values =
+      Array.map Det_coef.to_float
+        d.Determinants_by_hand.psi_coef
+    in
+    let hash = 
+      Hashtbl.create n_det
+    in
+    for i=0 to n_det-1
+    do
+      Hashtbl.add hash keys.(i) values.(i);
+    done;
+    hash
+  in
+
+  let overlap wf wf' =
+    let result, norm, norm' = 
+      Hashtbl.fold (fun k c (accu,norm,norm') ->
+        let c' =
+          try  Hashtbl.find wf' k 
+          with Not_found -> 0.
+        in
+        (accu +. c *. c' , 
+        norm +. c *. c  , 
+        norm'+. c'*. c' ) 
+      ) wf (0.,0.,0.)
+    in 
+    result /. (norm *. norm')
+  in
+
+  let wf, wf' = 
+     fetch_wf ezfio,
+     fetch_wf ezfio'
+  in
+
+  let o = 
+    overlap wf wf'
+  in
+  print_float (abs_float o)
+  

ocaml/qp_run.ml

@@ -15,7 +15,7 @@ let print_list () =
 let () = 
   Random.self_init ()
 
-let run ~master exe ezfio_file =
+let run slave exe ezfio_file =
 
 
   (** Check availability of the ports *)
@@ -28,7 +28,7 @@ let run ~master exe ezfio_file =
     in
     let rec try_new_port port_number =
       try 
-        List.iter [ 0;1;2;3;4 ] ~f:(fun i ->
+        List.iter [ 0;1;2;3;4;5;6;7;8;9 ] ~f:(fun i ->
             let address = 
               Printf.sprintf "tcp://%s:%d" (Lazy.force TaskServer.ip_address) (port_number+i)
             in
@@ -43,6 +43,7 @@ let run ~master exe ezfio_file =
       try_new_port 41279
     in
     ZMQ.Socket.close dummy_socket;
+    ZMQ.Context.terminate zmq_context;
     result
   in
   let time_start = 
@@ -74,16 +75,23 @@ let run ~master exe ezfio_file =
     | 0 -> ()
     | i -> failwith "Error: Input inconsistent\n"
   end;
-  begin
-    match master with
-    | Some address -> Unix.putenv ~key:"QP_RUN_ADDRESS_MASTER" ~data:address
-    | None -> ()
-  end;
 
-  (** Start task server *)
-  let address =
-    Printf.sprintf "tcp://%s:%d" (Lazy.force TaskServer.ip_address) port_number
+  let qp_run_address_filename = 
+   Filename.concat (Qpackage.ezfio_work ezfio_file) "qp_run_address"
   in
+
+  let () = 
+    if slave then
+      try
+        let address = 
+          In_channel.read_all qp_run_address_filename
+          |> String.strip
+        in
+        Unix.putenv ~key:"QP_RUN_ADDRESS_MASTER" ~data:address
+      with Sys_error _ -> failwith "No master is not running"
+  in
+       
+  (** Start task server *)
   let task_thread =
      let thread = 
       Thread.create ( fun () -> 
@@ -91,7 +99,16 @@ let run ~master exe ezfio_file =
      in
      thread ();
   in
+  let address =
+    Printf.sprintf "tcp://%s:%d" (Lazy.force TaskServer.ip_address) port_number
+  in
   Unix.putenv ~key:"QP_RUN_ADDRESS" ~data:address;
+  let () = 
+    if (not slave) then
+      Out_channel.with_file qp_run_address_filename  ~f:(
+        fun oc -> Out_channel.output_lines oc [address]) 
+  in
+
 
   (** Run executable *)
   let prefix = 
@@ -110,6 +127,8 @@ let run ~master exe ezfio_file =
 
   TaskServer.stop ~port:port_number;
   Thread.join task_thread;
+  if (not slave) then
+    Sys.remove qp_run_address_filename;
 
   let duration = Time.diff (Time.now()) time_start 
   |> Core.Span.to_string in
@@ -118,8 +137,8 @@ let run ~master exe ezfio_file =
 let spec = 
   let open Command.Spec in
   empty
-  +> flag "master" (optional string)
-     ~doc:("address Address of the master process")
+  +> flag "slave" no_arg
+     ~doc:(" Needed for slave tasks")
   +> anon ("executable" %: string)
   +> anon ("ezfio_file" %: string)
 ;;
@@ -137,8 +156,8 @@ Executes a Quantum Package binary file among these:\n\n"
     )
   )
   spec
-  (fun master exe ezfio_file () ->
-    run ~master exe ezfio_file
+  (fun slave exe ezfio_file () ->
+    run slave exe ezfio_file
   )
   |> Command.run   ~version: Git.sha1   ~build_info: Git.message
 

ocaml/qptypes_generator.ml

@@ -47,12 +47,8 @@ let input_data = "
 
 * States_number : int 
   assert (x > 0) ; 
-  if (x > 100) then
-    warning \"More than 100 states\";
-  if (Ezfio.has_determinants_n_states_diag ()) then
-    assert (x <= (Ezfio.get_determinants_n_states_diag ()))
-  else if (Ezfio.has_determinants_n_states ()) then
-    assert (x <= (Ezfio.get_determinants_n_states ()));
+  if (x > 1000) then
+    warning \"More than 1000 states\";
 
 * Bit_kind_size : int  
   begin match x with

plugins/All_singles/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Generators_restart Perturbation Properties Selectors_no_sorted Utils
+Generators_restart Perturbation Properties Selectors_no_sorted Utils Davidson

plugins/All_singles/README.rst

@@ -6,7 +6,77 @@ Needed Modules
 ==============
 .. Do not edit this section It was auto-generated
 .. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `Generators_restart <http://github.com/LCPQ/quantum_package/tree/master/plugins/Generators_restart>`_
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation>`_
+* `Properties <http://github.com/LCPQ/quantum_package/tree/master/plugins/Properties>`_
+* `Selectors_no_sorted <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_no_sorted>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+
 Documentation
 =============
 .. Do not edit this section It was auto-generated
 .. by the `update_README.py` script.
+
+
+h_apply_just_1h_1p
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_just_1h_1p_diexc
+  Undocumented
+
+
+h_apply_just_1h_1p_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_just_1h_1p_diexcp
+  Undocumented
+
+
+h_apply_just_1h_1p_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_just_mono
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_just_mono_diexc
+  Undocumented
+
+
+h_apply_just_mono_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_just_mono_diexcp
+  Undocumented
+
+
+h_apply_just_mono_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+`restart_more_singles <http://github.com/LCPQ/quantum_package/tree/master/plugins/All_singles/all_singles.irp.f#L1>`_
+  Generates and select single excitations
+  on the top of a given restart wave function
+
+
+`routine <http://github.com/LCPQ/quantum_package/tree/master/plugins/All_singles/all_singles.irp.f#L11>`_
+  Undocumented
+

plugins/CAS_SD/.gitignore

@@ -22,6 +22,9 @@ Properties
 Pseudo
 Selectors_full
 Utils
+ZMQ
+cas_s
+cas_s_selected
 cas_sd
 cas_sd_selected
 ezfio_interface.irp.f

plugins/CAS_SD/H_apply.irp.f

@@ -3,6 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import *
 
 s = H_apply("CAS_SD")
+s.unset_skip()
 print s
 
 s = H_apply("CAS_SD_selected_no_skip")
@@ -12,6 +13,7 @@ print s
 
 s = H_apply("CAS_SD_selected")
 s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_skip()
 print s
 
 s = H_apply("CAS_SD_PT2")
@@ -22,13 +24,9 @@ print s
 s = H_apply("CAS_S",do_double_exc=False)
 print s
 
-s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
-s.set_selection_pt2("epstein_nesbet_2x2")
-s.unset_skip()
-print s
-
 s = H_apply("CAS_S_selected",do_double_exc=False)
 s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_skip()
 print s
 
 s = H_apply("CAS_S_PT2",do_double_exc=False)

plugins/CAS_SD/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_full Generators_CAS
+Perturbation Selectors_full Generators_CAS Davidson

plugins/CAS_SD/README.rst

@@ -118,6 +118,106 @@ Documentation
   Undocumented
 
 
+h_apply_cas_s
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cas_s_diexc
+  Undocumented
+
+
+h_apply_cas_s_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cas_s_diexcp
+  Undocumented
+
+
+h_apply_cas_s_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cas_s_pt2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cas_s_pt2_diexc
+  Undocumented
+
+
+h_apply_cas_s_pt2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cas_s_pt2_diexcp
+  Undocumented
+
+
+h_apply_cas_s_pt2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cas_s_selected
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cas_s_selected_diexc
+  Undocumented
+
+
+h_apply_cas_s_selected_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cas_s_selected_diexcp
+  Undocumented
+
+
+h_apply_cas_s_selected_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cas_s_selected_no_skip
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cas_s_selected_no_skip_diexc
+  Undocumented
+
+
+h_apply_cas_s_selected_no_skip_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cas_s_selected_no_skip_diexcp
+  Undocumented
+
+
+h_apply_cas_s_selected_no_skip_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
 h_apply_cas_sd
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.

plugins/CAS_SD/cas_s_selected.irp.f

@@ -12,6 +12,7 @@ program full_ci
   
   pt2 = 1.d0
   diag_algorithm = "Lapack"
+
   if (N_det > N_det_max) then
     call diagonalize_CI
     call save_wavefunction
@@ -28,49 +29,84 @@ program full_ci
     print *,  'E+PT2    = ', CI_energy+pt2
     print *,  '-----'
   endif
+  double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
+  double precision :: E_CI_before(N_states)
+  if(read_wf)then
+   call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
+   h = diag_H_mat_elem(psi_det(1,1,N_det),N_int)
+   selection_criterion = dabs(psi_coef(N_det,1) *  (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0
+   soft_touch selection_criterion
+  endif
 
+
+  integer :: n_det_before
+  print*,'Beginning the selection ...'
+  E_CI_before(1:N_states) = CI_energy(1:N_states)
   do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
-    call H_apply_CAS_S_selected_no_skip(pt2, norm_pert, H_pert_diag,  N_st)
+    n_det_before = N_det
+    call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag,  N_st)
 
     PROVIDE  psi_coef
     PROVIDE  psi_det
     PROVIDE  psi_det_sorted
 
-    if (N_det > N_det_max) then
-       psi_det = psi_det_sorted
-       psi_coef = psi_coef_sorted
-       N_det = N_det_max
-       soft_touch N_det psi_det psi_coef
-    endif
     call diagonalize_CI
+
+    if (N_det > N_det_max) then
+      N_det = N_det_max
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
+    endif
+
+
     call save_wavefunction
-    print *,  'N_det    = ', N_det
-    print *,  'N_states = ', N_states
-    print *,  'PT2      = ', pt2
-    print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
+    if(n_det_before == N_det)then
+     selection_criterion = selection_criterion * 0.5d0
+    endif
+    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 *,  '-----'
+    if(N_states.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
+     enddo
+    endif
+    E_CI_before(1:N_states) = CI_energy(1:N_states)
     call ezfio_set_cas_sd_energy(CI_energy(1))
   enddo
-  call diagonalize_CI
-
+   N_det = min(N_det_max,N_det)
+   touch N_det psi_det psi_coef
+   call diagonalize_CI
    if(do_pt2_end)then
     print*,'Last iteration only to compute the PT2'
     threshold_selectors = 1.d0
     threshold_generators = 0.999d0
-    call H_apply_CAS_S_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+    call H_apply_CAS_SD_PT2(pt2, norm_pert, H_pert_diag,  N_st)
 
     print *,  'Final step'
     print *,  'N_det    = ', N_det
     print *,  'N_states = ', N_states
     print *,  'PT2      = ', pt2
-    print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  'E        = ', CI_energy(1:N_states)
+    print *,  'E+PT2    = ', CI_energy(1:N_states)+pt2(1:N_states)
     print *,  '-----'
     call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
    endif
 
-
   integer :: exc_max, degree_min
   exc_max = 0
   print *,  'CAS determinants : ', N_det_cas
@@ -79,6 +115,7 @@ program full_ci
       call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
       exc_max = max(exc_max,degree)
     enddo
+    print *,  psi_coef_cas_diagonalized(i,:)
     call debug_det(psi_cas(1,1,i),N_int)
     print *,  ''
   enddo

plugins/CAS_SD/cas_sd.irp.f

@@ -1,7 +1,6 @@
 program full_ci
   implicit none
   integer                        :: i,k
-  integer                        :: N_det_old
 
   
   double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
@@ -11,9 +10,9 @@ program full_ci
   character*(64)                 :: perturbation
   PROVIDE N_det_cas
   
-  N_det_old = 0
   pt2 = 1.d0
   diag_algorithm = "Lapack"
+
   if (N_det > N_det_max) then
     call diagonalize_CI
     call save_wavefunction
@@ -30,36 +29,68 @@ program full_ci
     print *,  'E+PT2    = ', CI_energy+pt2
     print *,  '-----'
   endif
+  double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
+  double precision :: E_CI_before(N_states)
+  if(read_wf)then
+   call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
+   h = diag_H_mat_elem(psi_det(1,1,N_det),N_int)
+   selection_criterion = dabs(psi_coef(N_det,1) *  (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0
+   soft_touch selection_criterion
+  endif
 
+
+  integer :: n_det_before
+  print*,'Beginning the selection ...'
+  E_CI_before(1:N_states) = CI_energy(1:N_states)
   do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
-    N_det_old = N_det
+    n_det_before = N_det
     call H_apply_CAS_SD(pt2, norm_pert, H_pert_diag,  N_st)
 
     PROVIDE  psi_coef
     PROVIDE  psi_det
     PROVIDE  psi_det_sorted
 
-    if (N_det > N_det_max) then
-       psi_det = psi_det_sorted
-       psi_coef = psi_coef_sorted
-       N_det = N_det_max
-       soft_touch N_det psi_det psi_coef
-    endif
     call diagonalize_CI
-    call save_wavefunction
-    print *,  'N_det    = ', N_det
-    print *,  'N_states = ', N_states
-    print *,  'PT2      = ', pt2
-    print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
-    print *,  '-----'
-    call ezfio_set_cas_sd_energy(CI_energy(1))
-    if (N_det == N_det_old) then
-      exit
-    endif
-  enddo
-  call diagonalize_CI
 
+    if (N_det > N_det_max) then
+      N_det = N_det_max
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
+    endif
+
+
+    call save_wavefunction
+    if(n_det_before == N_det)then
+     selection_criterion = selection_criterion * 0.5d0
+    endif
+    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 *,  '-----'
+    if(N_states.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
+     enddo
+    endif
+    E_CI_before(1:N_states) = CI_energy(1:N_states)
+    call ezfio_set_cas_sd_energy(CI_energy(1))
+  enddo
+   N_det = min(N_det_max,N_det)
+   touch N_det psi_det psi_coef
+   call diagonalize_CI
    if(do_pt2_end)then
     print*,'Last iteration only to compute the PT2'
     threshold_selectors = 1.d0
@@ -70,13 +101,12 @@ program full_ci
     print *,  'N_det    = ', N_det
     print *,  'N_states = ', N_states
     print *,  'PT2      = ', pt2
-    print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  'E        = ', CI_energy(1:N_states)
+    print *,  'E+PT2    = ', CI_energy(1:N_states)+pt2(1:N_states)
     print *,  '-----'
     call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
    endif
 
-
   integer :: exc_max, degree_min
   exc_max = 0
   print *,  'CAS determinants : ', N_det_cas
@@ -85,6 +115,7 @@ program full_ci
       call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
       exc_max = max(exc_max,degree)
     enddo
+    print *,  psi_coef_cas_diagonalized(i,:)
     call debug_det(psi_cas(1,1,i),N_int)
     print *,  ''
   enddo

plugins/CAS_SD/cas_sd_selected.irp.f

@@ -12,6 +12,7 @@ program full_ci
   
   pt2 = 1.d0
   diag_algorithm = "Lapack"
+
   if (N_det > N_det_max) then
     call diagonalize_CI
     call save_wavefunction
@@ -28,32 +29,68 @@ program full_ci
     print *,  'E+PT2    = ', CI_energy+pt2
     print *,  '-----'
   endif
+  double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
+  double precision :: E_CI_before(N_states)
+  if(read_wf)then
+   call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
+   h = diag_H_mat_elem(psi_det(1,1,N_det),N_int)
+   selection_criterion = dabs(psi_coef(N_det,1) *  (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0
+   soft_touch selection_criterion
+  endif
 
+
+  integer :: n_det_before
+  print*,'Beginning the selection ...'
+  E_CI_before(1:N_states) = CI_energy(1:N_states)
   do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    n_det_before = N_det
     call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag,  N_st)
 
     PROVIDE  psi_coef
     PROVIDE  psi_det
     PROVIDE  psi_det_sorted
 
-    if (N_det > N_det_max) then
-       psi_det = psi_det_sorted
-       psi_coef = psi_coef_sorted
-       N_det = N_det_max
-       soft_touch N_det psi_det psi_coef
-    endif
     call diagonalize_CI
+
+    if (N_det > N_det_max) then
+      N_det = N_det_max
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
+    endif
+
+
     call save_wavefunction
-    print *,  'N_det    = ', N_det
-    print *,  'N_states = ', N_states
-    print *,  'PT2      = ', pt2
-    print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
+    if(n_det_before == N_det)then
+     selection_criterion = selection_criterion * 0.5d0
+    endif
+    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 *,  '-----'
+    if(N_states.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
+     enddo
+    endif
+    E_CI_before(1:N_states) = CI_energy(1:N_states)
     call ezfio_set_cas_sd_energy(CI_energy(1))
   enddo
-  call diagonalize_CI
-
+   N_det = min(N_det_max,N_det)
+   touch N_det psi_det psi_coef
+   call diagonalize_CI
    if(do_pt2_end)then
     print*,'Last iteration only to compute the PT2'
     threshold_selectors = 1.d0
@@ -64,13 +101,12 @@ program full_ci
     print *,  'N_det    = ', N_det
     print *,  'N_states = ', N_states
     print *,  'PT2      = ', pt2
-    print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  'E        = ', CI_energy(1:N_states)
+    print *,  'E+PT2    = ', CI_energy(1:N_states)+pt2(1:N_states)
     print *,  '-----'
     call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
    endif
 
-
   integer :: exc_max, degree_min
   exc_max = 0
   print *,  'CAS determinants : ', N_det_cas
@@ -79,6 +115,7 @@ program full_ci
       call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
       exc_max = max(exc_max,degree)
     enddo
+    print *,  psi_cas_coef(i,:)
     call debug_det(psi_cas(1,1,i),N_int)
     print *,  ''
   enddo

plugins/CID/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Selectors_full SingleRefMethod
+Selectors_full SingleRefMethod Davidson

plugins/CIS/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Selectors_full SingleRefMethod
+Selectors_full SingleRefMethod Davidson

plugins/CISD/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Selectors_full SingleRefMethod
+Selectors_full SingleRefMethod Davidson

plugins/Casino/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Determinants
+Determinants Davidson

plugins/DDCI_selected/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_full Generators_CAS
+Perturbation Selectors_full Generators_CAS Davidson

plugins/FCIdump/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Determinants
+Determinants Davidson

plugins/FOBOCI/H_apply_dressed_autonom.irp.f

@@ -273,7 +273,7 @@ subroutine H_apply_dressed_pert_monoexc(key_in, hole_1,particl_1,i_generator,ipr
   integer,parameter              :: size_max = 3072
   
     integer, intent(in) :: Ndet_generators
-    double precision, intent(in) :: E_ref
+    double precision, intent(inout) :: E_ref
     double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
     integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
 
@@ -438,7 +438,7 @@ subroutine H_apply_dressed_pert(delta_ij_generators_,  Ndet_generators,psi_det_g
   
   
     integer, intent(in) :: Ndet_generators
-    double precision, intent(in) :: E_ref
+    double precision, intent(inout) :: E_ref
     double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
     integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
 

plugins/FOBOCI/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_no_sorted Hartree_Fock
+Perturbation Selectors_no_sorted Hartree_Fock Davidson CISD

plugins/FOBOCI/dress_simple.irp.f

@@ -207,16 +207,16 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
 
  call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators)  ! Diagonalize the Dressed_H_matrix
  
- double precision :: s2,E_ref(N_states)
+ double precision :: s2(N_det_generators),E_ref(N_states)
  integer :: i_state(N_states)
  integer :: n_state_good
  n_state_good = 0
  if(s2_eig)then
+  call u_0_S2_u_0(s2,eigvectors,Ndet_generators,psi_det_generators_input,N_int,N_det_generators,size(eigvectors,1))
   do i = 1, Ndet_generators
-    call get_s2_u0(psi_det_generators_input,eigvectors(1,i),Ndet_generators,Ndet_generators,s2)
-    print*,'s2 = ',s2
-    print*,dabs(s2-expected_s2)
-    if(dabs(s2-expected_s2).le.0.3d0)then
+    print*,'s2 = ',s2(i)
+    print*,dabs(s2(i)-expected_s2)
+    if(dabs(s2(i)-expected_s2).le.0.3d0)then
      n_state_good +=1
      i_state(n_state_good) = i
      E_ref(n_state_good) = eigvalues(i)
@@ -274,7 +274,6 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
  integer :: i_good_state(0:N_states)
  i_good_state(0) = 0
   do i = 1, Ndet_generators
-    call get_s2_u0(psi_det_generators_input,eigvectors(1,i),Ndet_generators,Ndet_generators,s2)
     ! State following
     do k = 1, N_states 
      accu = 0.d0

plugins/FOBOCI/hcc_1h1p.irp.f

@@ -15,11 +15,10 @@ subroutine routine
   call diagonalize_CI
   call test_hcc
   call test_mulliken
-! call SC2_1h1p(psi_det,psi_coef,energies, &
-!       diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
   allocate(H_matrix(N_det,N_det))
-  call SC2_1h1p_full(psi_det,psi_coef,energies, &
-        H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
+  stop 'SC2_1h1p_full is not in the git!'
+!  call SC2_1h1p_full(psi_det,psi_coef,energies, &
+!        H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
   deallocate(H_matrix)
   integer :: i,j
   double precision :: accu,coef_hf

plugins/FOBOCI/routines_foboci.irp.f

@@ -799,7 +799,7 @@ end
   call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det)
   call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,i_hole)
   call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,i_hole)
-  call davidson_diag_hjj(psi_det,psi_coef,dressing_H_mat_elem,energies,size(psi_coef,1),N_det,N_states_diag,N_int,output_determinants)
+  call davidson_diag_hjj(psi_det,psi_coef,dressing_H_mat_elem,energies,size(psi_coef,1),N_det,N_states,N_states_diag,N_int,output_determinants)
   do i = 1, 2
    print*,'psi_coef = ',psi_coef(i,1)
   enddo

plugins/Full_CI/.gitignore

@@ -28,6 +28,7 @@ full_ci
 full_ci_no_skip
 irpf90.make
 irpf90_entities
+micro_pt2
 tags
 target_pt2
 var_pt2_ratio

plugins/Full_CI/H_apply.irp.f

@@ -23,6 +23,11 @@ s.unset_skip()
 #s.unset_openmp()
 print s
 
+s = H_apply("FCI_no_selection")
+s.set_selection_pt2("dummy")
+s.unset_skip()
+print s
+
 s = H_apply("FCI_mono")
 s.set_selection_pt2("epstein_nesbet_2x2")
 s.unset_double_excitations()
@@ -30,28 +35,6 @@ s.unset_openmp()
 print s
 
 
-s = H_apply("select_mono_delta_rho")
-s.unset_double_excitations()
-s.set_selection_pt2("delta_rho_one_point")
-s.unset_openmp()
-print s
-
-s = H_apply("pt2_mono_delta_rho")
-s.unset_double_excitations()
-s.set_perturbation("delta_rho_one_point")
-s.unset_openmp()
-print s
-
-s = H_apply("select_mono_di_delta_rho")
-s.set_selection_pt2("delta_rho_one_point")
-s.unset_openmp()
-print s
-
-s = H_apply("pt2_mono_di_delta_rho")
-s.set_perturbation("delta_rho_one_point")
-s.unset_openmp()
-print s
-
 
 END_SHELL
 

plugins/Full_CI/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_full Generators_full
+Perturbation Selectors_full Generators_full Davidson

plugins/Full_CI/README.rst

@@ -107,6 +107,10 @@ h_apply_fci_pt2
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
+h_apply_fci_pt2_collector
+  Collects results from the selection in an array of generators
+
+
 h_apply_fci_pt2_diexc
   Undocumented
 
@@ -127,6 +131,19 @@ h_apply_fci_pt2_monoexc
   Assume N_int is already provided.
 
 
+h_apply_fci_pt2_slave
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_fci_pt2_slave_inproc
+  Computes a buffer using threads
+
+
+h_apply_fci_pt2_slave_tcp
+  Computes a buffer over the network
+
+
 h_apply_pt2_mono_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
@@ -227,6 +244,18 @@ h_apply_select_mono_di_delta_rho_monoexc
   Assume N_int is already provided.
 
 
+`micro_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/micro_pt2.irp.f#L1>`_
+  Helper program to compute the PT2 in distributed mode.
+
+
+`provide_everything <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/micro_pt2.irp.f#L15>`_
+  Undocumented
+
+
+`run_wf <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/micro_pt2.irp.f#L19>`_
+  Undocumented
+
+
 `var_pt2_ratio_run <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/var_pt2_ratio.irp.f#L1>`_
   Undocumented
 

plugins/Full_CI/full_ci.irp.f

@@ -11,7 +11,7 @@ program full_ci
   
   pt2 = 1.d0
   diag_algorithm = "Lapack"
-  
+
   if (N_det > N_det_max) then
     call diagonalize_CI
     call save_wavefunction
@@ -40,7 +40,7 @@ program full_ci
 
   integer :: n_det_before
   print*,'Beginning the selection ...'
-  E_CI_before = CI_energy
+  E_CI_before(1:N_states) = CI_energy(1:N_states)
   do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
     n_det_before = N_det
     call H_apply_FCI(pt2, norm_pert, H_pert_diag,  N_st)
@@ -49,13 +49,16 @@ program full_ci
     PROVIDE  psi_det
     PROVIDE  psi_det_sorted
 
-    if (N_det > N_det_max) then
-       psi_det = psi_det_sorted
-       psi_coef = psi_coef_sorted
-       N_det = N_det_max
-       soft_touch N_det psi_det psi_coef
-    endif
     call diagonalize_CI
+
+    if (N_det > N_det_max) then
+      N_det = N_det_max
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
+    endif
+
+
     call save_wavefunction
     if(n_det_before == N_det)then
      selection_criterion = selection_criterion * 0.5d0
@@ -69,7 +72,6 @@ program full_ci
     print *,  'E(before)+PT2  = ', E_CI_before(k)+pt2(k)
     enddo
     print *,  '-----'
-    E_CI_before = CI_energy
     if(N_states.gt.1)then
      print*,'Variational Energy difference'
      do i = 2, N_states
@@ -82,8 +84,8 @@ program full_ci
       print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
      enddo
     endif
-    E_CI_before = CI_energy
-    call ezfio_set_full_ci_energy(CI_energy)
+    E_CI_before(1:N_states) = CI_energy(1:N_states)
+    call ezfio_set_full_ci_energy(CI_energy(1))
   enddo
    N_det = min(N_det_max,N_det)
    touch N_det psi_det psi_coef
@@ -98,10 +100,10 @@ program full_ci
     print *,  'N_det    = ', N_det
     print *,  'N_states = ', N_states
     print *,  'PT2      = ', pt2
-    print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  'E        = ', CI_energy(1:N_states)
+    print *,  'E+PT2    = ', CI_energy(1:N_states)+pt2(1:N_states)
     print *,  '-----'
-    call ezfio_set_full_ci_energy_pt2(CI_energy+pt2)
+    call ezfio_set_full_ci_energy_pt2(CI_energy(1)+pt2(1))
    endif
    call save_wavefunction
   deallocate(pt2,norm_pert)

plugins/Full_CI/micro_pt2.irp.f

@@ -1,46 +0,0 @@
-program micro_pt2
-  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
-
-  print *,  'Getting wave function'
-  zmq_context = f77_zmq_ctx_new ()
-
-  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
-
-  call zmq_get_psi(zmq_to_qp_run_socket, 1)
-  call write_double(6,ci_energy,'Energy')
-  zmq_state = 'h_apply_fci_pt2'
-
-  call provide_everything
-  integer :: rc, i
-
-  print *,  'Contribution to PT2 running'
-
-  !$OMP PARALLEL PRIVATE(i)
-  i = omp_get_thread_num()
-  call H_apply_FCI_PT2_slave_tcp(i)
-  !$OMP END PARALLEL
-
-
-end

plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Perturbation Selectors_full Generators_full ZMQ Full_CI

plugins/Full_CI_ZMQ/fci_zmq.irp.f

@@ -0,0 +1,220 @@
+program fci_zmq
+  implicit none
+  integer                        :: i,j,k
+  logical, external              :: detEq
+  
+  double precision, allocatable  :: pt2(:)
+  integer                        :: degree
+  
+  allocate (pt2(N_states))
+  
+  pt2 = 1.d0
+  diag_algorithm = "Lapack"
+  
+  if (N_det > N_det_max) then
+    call diagonalize_CI
+    call save_wavefunction
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = N_det_max
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    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+PT2    = ', CI_energy(k) + pt2(k)
+      print *,  '-----'
+    enddo
+  endif
+  double precision               :: E_CI_before(N_states)
+  
+  
+  integer                        :: n_det_before
+  print*,'Beginning the selection ...'
+  E_CI_before(1:N_states) = CI_energy(1:N_states)
+  
+  do while ( (N_det < N_det_max) .and. (maxval(abs(pt2(1:N_states))) > pt2_max) )
+    n_det_before = N_det
+    call ZMQ_selection(max(1024-N_det, N_det), pt2)
+    
+    PROVIDE  psi_coef
+    PROVIDE  psi_det
+    PROVIDE  psi_det_sorted
+
+    call diagonalize_CI
+    call save_wavefunction
+    
+    if (N_det > N_det_max) then
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      N_det = N_det_max
+      soft_touch N_det psi_det psi_coef
+      call diagonalize_CI
+      call save_wavefunction
+    endif
+    
+    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 *,  '-----'
+    if(N_states.gt.1)then
+      print*,'Variational Energy difference'
+      do i = 2, N_states
+        print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+      enddo
+    endif
+    if(N_states.gt.1)then
+      print*,'Variational + perturbative Energy difference'
+      do i = 2, N_states
+        print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
+      enddo
+    endif
+    E_CI_before(1:N_states) = CI_energy(1:N_states)
+    call ezfio_set_full_ci_energy(CI_energy)
+  enddo
+
+  if(do_pt2_end)then
+    print*,'Last iteration only to compute the PT2'
+    threshold_selectors = 1.d0
+    threshold_generators = 0.9999d0
+    E_CI_before(1:N_states) = CI_energy(1:N_states)
+    call ZMQ_selection(0, pt2)
+    print *,  'Final step'
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    do k=1,N_states
+      print *, 'State', k
+      print *,  'PT2      = ', pt2
+      print *,  'E        = ', E_CI_before
+      print *,  'E+PT2    = ', E_CI_before+pt2
+      print *,  '-----'
+    enddo
+    call ezfio_set_full_ci_energy_pt2(E_CI_before+pt2)
+  endif
+  call save_wavefunction
+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)
+  
+  
+  N = max(N_in,1)
+  provide nproc
+  provide ci_electronic_energy
+  call new_parallel_job(zmq_to_qp_run_socket,"selection")
+  call zmq_put_psi(zmq_to_qp_run_socket,1,ci_electronic_energy,size(ci_electronic_energy))
+  call zmq_set_running(zmq_to_qp_run_socket)
+  call create_selection_buffer(N, N*2, b)
+
+  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 = max(i-step+1,1)
+    i_generator_max = i
+    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(none)  SHARED(b, pt2)  PRIVATE(i) NUM_THREADS(nproc+1) shared(ci_electronic_energy_is_built, n_det_generators_is_built, n_states_is_built, n_int_is_built, nproc_is_built)
+      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()
+  endif
+end subroutine
+
+
+subroutine selection_slave_inproc(i)
+  implicit none
+  integer, intent(in)            :: i
+
+  call run_selection_slave(1,i,ci_electronic_energy)
+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/run_selection_slave.irp.f

@@ -0,0 +1,156 @@
+
+subroutine run_selection_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)
+
+  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 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
+  end if
+  buf%N = 0
+  ctask = 1
+  pt2 = 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_generator_start, i_generator_max, step, N
+      read (task,*) i_generator_start, i_generator_max, step, 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)
+      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
+    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_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))
+        enddo
+        call sort_selection_buffer(buf2)
+        buf%mini = buf2%mini
+        pt2 = 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_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
+  use f77_zmq
+  use selection_types
+  implicit none
+
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_push
+  double precision, intent(in)   :: pt2(N_states)
+  type(selection_buffer), intent(inout) :: b
+  integer, intent(in) :: ntask, task_id(*)
+  integer :: rc
+
+  call sort_selection_buffer(b)
+
+  rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE)
+  if(rc /= 4) stop "push"
+  rc = f77_zmq_send( zmq_socket_push, pt2, 8*N_states, ZMQ_SNDMORE)
+  if(rc /= 8*N_states) stop "push"
+
+  rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE)
+  if(rc /= 8*b%cur) stop "push"
+
+  rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE)
+  if(rc /= bit_kind*N_int*2*b%cur) 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(1), ntask*4, 0)
+  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)
+end subroutine
+
+
+subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, ntask)
+  use f77_zmq
+  use selection_types
+  implicit none
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
+  double precision, intent(inout) :: pt2(N_states)
+  double precision, intent(out) :: val(*)
+  integer(bit_kind), intent(out) :: det(N_int, 2, *)
+  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, pt2, N_states*8, 0)
+  if(rc /= 8*N_states) stop "pull"
+
+  rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0)
+  if(rc /= 8*N) stop "pull"
+
+  rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0)
+  if(rc /= bit_kind*N_int*2*N) stop "pull"
+
+  rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
+  if(rc /= 4) stop "pull"
+
+  rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)
+  if(rc /= 4*ntask) stop "pull"
+
+! Activate is zmq_socket_pull is a REP
+!  rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
+end subroutine
+ 
+ 
+

plugins/Full_CI_ZMQ/selection.irp.f

@@ -0,0 +1,106 @@
+use bitmasks
+
+
+double precision function integral8(i,j,k,l)
+  implicit none
+  
+  integer, intent(in) :: i,j,k,l
+  double precision, external :: get_mo_bielec_integral
+  
+  integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
+end function
+
+
+BEGIN_PROVIDER [ integer(1), psi_phasemask, (N_int*bit_kind_size, 2, N_det)]
+  use bitmasks
+  implicit none
+  
+  integer :: i
+  do i=1, N_det
+    call get_mask_phase(psi_det_sorted(1,1,i), psi_phasemask(1,1,i))
+  end do
+END_PROVIDER
+
+
+subroutine assert(cond, msg)
+  character(*), intent(in) :: msg
+  logical, intent(in) :: cond
+  
+  if(.not. cond) then
+    print *, "assert fail: "//msg
+    stop
+  end if
+end subroutine
+
+
+subroutine get_mask_phase(det, phasemask)
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in) :: det(N_int, 2)
+  integer(1), intent(out) :: phasemask(N_int*bit_kind_size, 2)
+  integer :: s, ni, i
+  logical :: change
+
+  phasemask = 0_1
+  do s=1,2
+    change = .false.
+    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
+      end do
+    end do
+  end do
+end subroutine
+
+
+subroutine select_connected(i_generator,E0,pt2,b)
+  use bitmasks
+  use selection_types
+  implicit none
+  integer, intent(in)            :: i_generator
+  type(selection_buffer), intent(inout) :: b
+  double precision, intent(inout)  :: pt2(N_states)
+  integer :: k,l
+  double precision, intent(in)   :: E0(N_states)
+
+  integer(bit_kind)              :: hole_mask(N_int,2), particle_mask(N_int,2)
+  double precision               :: fock_diag_tmp(2,mo_tot_num+1)
+  
+  call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
+
+  do l=1,N_generators_bitmask
+    do k=1,N_int
+      hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole,l), psi_det_generators(k,1,i_generator))
+      hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole,l), psi_det_generators(k,2,i_generator))
+      particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) )
+      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)
+  enddo
+end subroutine
+
+
+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, intent(in) :: s1, s2, h1, h2, p1, p2
+  logical :: change
+  integer(1) :: np
+  double precision, parameter :: 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
+  
+  if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1_1
+  get_phase_bi = res(iand(np,1_1))
+end subroutine
+
+
+

plugins/Full_CI_ZMQ/selection_buffer.irp.f

@@ -0,0 +1,70 @@
+
+subroutine create_selection_buffer(N, siz, res)
+  use selection_types
+  implicit none
+
+  integer, intent(in) :: N, siz
+  type(selection_buffer), intent(out) :: res
+
+  allocate(res%det(N_int, 2, siz), res%val(siz))
+
+  res%val = 0d0
+  res%det = 0_8
+  res%N = N
+  res%mini = 0d0
+  res%cur = 0
+end subroutine
+
+  
+subroutine add_to_selection_buffer(b, det, val)
+  use selection_types
+  implicit none
+
+  type(selection_buffer), intent(inout) :: b
+  integer(bit_kind), intent(in) :: det(N_int, 2)
+  double precision, intent(in) :: val
+  integer :: i
+
+  if(dabs(val) >= b%mini) then
+    b%cur += 1
+    b%det(:,:,b%cur) = det(:,:)
+    b%val(b%cur) = val
+    if(b%cur == size(b%val)) then
+      call sort_selection_buffer(b)
+    end if
+  end if
+end subroutine
+
+
+subroutine sort_selection_buffer(b)
+  use selection_types
+  implicit none
+
+  type(selection_buffer), intent(inout) :: b
+  double precision, allocatable :: vals(:), absval(:)
+  integer, allocatable :: iorder(:)
+  integer(bit_kind), allocatable :: 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))
+  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))
+    vals(i) = b%val(iorder(i))
+  end do
+  b%det(:,:,:nmwen) = detmp(:,:,:)
+  b%det(:,:,nmwen+1:) = 0_bit_kind
+  b%val(:nmwen) = vals(:)
+  b%val(nmwen+1:) = 0d0
+  b%mini = max(b%mini,dabs(b%val(b%N)))
+  b%cur = nmwen
+end subroutine
+

plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f

@@ -0,0 +1,107 @@
+program selection_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 mo_mono_elec_integral
+!   PROVIDE ci_electronic_energy mo_tot_num N_int
+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(2)
+  integer :: rc, i
+  
+  call provide_everything
+  
+  zmq_context = f77_zmq_ctx_new ()
+  states(1) = 'selection'
+  states(2) = 'davidson'
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+
+  do
+
+    call wait_for_states(states,zmq_state,2)
+
+    if(trim(zmq_state) == 'Stopped') then
+
+      exit
+
+    else if (trim(zmq_state) == 'selection') then
+
+      ! Selection
+      ! ---------
+
+      print *,  'Selection'
+      call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states_diag)
+  
+      !$OMP PARALLEL PRIVATE(i)
+      i = omp_get_thread_num()
+      call selection_slave_tcp(i, energy)
+      !$OMP END PARALLEL
+      print *,  'Selection done'
+
+    else if (trim(zmq_state) == 'davidson') then
+
+      ! Davidson
+      ! --------
+
+      print *,  'Davidson'
+      call davidson_miniserver_get()
+      !$OMP PARALLEL PRIVATE(i)
+      i = omp_get_thread_num()
+      call davidson_slave_tcp(i)
+      !$OMP END PARALLEL
+      print *,  'Davidson done'
+
+    endif
+
+  end do
+end
+
+subroutine update_energy(energy)
+  implicit none
+  double precision, intent(in) :: energy(N_states_diag)
+  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,size(ci_electronic_energy)
+      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_diag)
+  integer, intent(in)            :: i
+
+  call run_selection_slave(0,i,energy)
+end
+

plugins/Full_CI_ZMQ/selection_double.irp.f

@@ -0,0 +1,726 @@
+ 
+subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf)
+  use bitmasks
+  use selection_types
+  implicit none
+  
+  integer, intent(in)            :: i_generator
+  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                :: mat(N_states, mo_tot_num, mo_tot_num)
+  integer                         :: h1,h2,s1,s2,s3,i1,i2,ib,sp,k,i,j,nt,ii
+  integer(bit_kind)               :: hole(N_int,2), particle(N_int,2), mask(N_int, 2), pmask(N_int, 2)
+  logical                         :: fullMatch, ok
+  
+  integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2)
+  integer,allocatable               :: preinteresting(:), prefullinteresting(:), interesting(:), fullinteresting(:)
+  integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
+  
+  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))
+  
+  do k=1,N_int
+    hole    (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1))
+    hole    (k,2) = iand(psi_det_generators(k,2,i_generator), hole_mask(k,2))
+    particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), particle_mask(k,1))
+    particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2))
+  enddo
+
+  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)
+
+  
+  preinteresting(0) = 0
+  prefullinteresting(0) = 0
+  
+  do i=1,N_int
+    negMask(i,1) = not(psi_det_generators(i,1,i_generator))
+    negMask(i,2) = not(psi_det_generators(i,2,i_generator))
+  end do
+  
+  do i=1,N_det
+    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))
+    end do
+
+    if(nt <= 4) then
+      if(i <= N_det_selectors) then
+        preinteresting(0) += 1
+        preinteresting(preinteresting(0)) = i
+      else if(nt <= 2) then
+        prefullinteresting(0) += 1
+        prefullinteresting(prefullinteresting(0)) = i
+      end if
+    end if
+  end do
+  
+        
+  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
+      
+      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))
+        end do
+        
+        if(nt <= 4) then
+          interesting(0) += 1
+          interesting(interesting(0)) = i
+          minilist(:,:,interesting(0)) = psi_det_sorted(:,:,i)
+          if(nt <= 2) then
+            fullinteresting(0) += 1
+            fullinteresting(fullinteresting(0)) = i
+            fullminilist(:,:,fullinteresting(0)) = psi_det_sorted(:,:,i)
+          end if
+        end if
+      end do
+      
+      do ii=1,prefullinteresting(0)
+        i = prefullinteresting(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))
+        end do
+
+        if(nt <= 2) then
+          fullinteresting(0) += 1
+          fullinteresting(fullinteresting(0)) = i
+          fullminilist(:,:,fullinteresting(0)) = psi_det_sorted(:,:,i)
+        end if
+      end do
+      
+      do s2=s1,2
+        sp = s1
+        if(s1 /= s2) sp = 3
+        
+        ib = 1
+        if(s1 == s2) ib = i1+1
+        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)
+          
+          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 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)
+        enddo
+      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
+  implicit none
+  
+  integer, intent(in) :: i_generator, sp, h1, h2
+  double precision, intent(in) :: mat(N_states, mo_tot_num, mo_tot_num)
+  logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, 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, j, istate
+  integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
+  double precision :: e_pert, delta_E, val, Hii, max_e_pert
+  double precision, external :: diag_H_mat_elem_fock
+  
+  logical, external :: detEq
+  
+  
+  if(sp == 3) then
+    s1 = 1
+    s2 = 2
+  else
+    s1 = sp
+    s2 = sp
+  end if
+  
+  call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
+  
+  do p1=1,mo_tot_num
+    if(bannedOrb(p1, s1)) cycle
+    ib = 1
+    if(sp /= 3) ib = p1+1
+    do p2=ib,mo_tot_num
+      if(bannedOrb(p2, s2)) cycle
+      if(banned(p1,p2)) cycle
+      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
+      
+      do istate=1,N_states
+        delta_E = E0(istate) - Hii
+        val = mat(istate, p1, p2)
+        if (delta_E < 0.d0) then
+          e_pert = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * val * val) - delta_E)
+        else
+          e_pert = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * val * val) - delta_E)
+        endif
+        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) then
+        call add_to_selection_buffer(buf, det, max_e_pert)
+      end if
+    end do
+  end do
+end subroutine
+
+
+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)
+  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
+  integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
+!   logical :: bandon
+!   
+!   bandon = .false.
+  mat = 0d0
+  
+  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 ! interesting(0)
+    !i = interesting(ii)
+    
+    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 > 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(sp == 3) then
+          banned(:,:,2) = transpose(banned(:,:,1))
+        else
+          do k=1,mo_tot_num
+          do l=k+1,mo_tot_num
+            banned(l,k,1) = banned(k,l,1)
+          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
+  end do
+end subroutine
+
+
+subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
+  use bitmasks
+  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)
+  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)
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+  
+  double precision, external :: get_phase_bi, integral8
+  
+  integer :: i, j, tip, ma, mi, puti, putj
+  integer :: h1, h2, p1, p2, i1, i2
+  double precision :: hij, phase
+  
+  integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0,  3,4, 0,0, 0,0,  2,4, 1,4, 0,0,  2,3, 1,3, 1,2 /), (/2,3,4/))
+  integer, parameter :: turn2(2) = (/2, 1/)
+  integer, parameter :: turn3(2,3) = reshape((/2,3,  1,3, 1,2/), (/2,3/))
+  
+  integer :: bant
+  bant = 1
+
+  tip = p(0,1) * p(0,2)
+  
+  ma = sp
+  if(p(0,1) > p(0,2)) ma = 1
+  if(p(0,1) < p(0,2)) ma = 2
+  mi = mod(ma, 2) + 1
+  
+  if(sp == 3) then
+    if(ma == 2) bant = 2
+    
+    if(tip == 3) then
+      puti = p(1, mi)
+      do i = 1, 3
+        putj = p(i, ma)
+        if(banned(putj,puti,bant)) cycle
+        i1 = turn3(1,i)
+        i2 = turn3(2,i)
+        p1 = p(i1, ma)
+        p2 = p(i2, ma)
+        h1 = h(1, ma)
+        h2 = h(2, ma)
+        
+        hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
+        if(ma == 1) then
+          mat(:, putj, puti) += coefs * hij
+        else
+          mat(:, puti, putj) += coefs * hij
+        end if
+      end do
+    else
+      do i = 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
+      end do
+    end if
+
+  else
+    if(tip == 0) then
+      h1 = h(1, ma)
+      h2 = h(2, ma)
+      do i=1,3
+      puti = p(i, ma)
+      do j=i+1,4
+        putj = p(j, ma)
+        if(banned(puti,putj,1)) cycle
+        
+        i1 = turn2d(1, i, j)
+        i2 = turn2d(2, i, j)
+        p1 = p(i1, ma)
+        p2 = p(i2, ma)
+        hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
+        mat(:, puti, putj) += coefs * hij
+      end do
+      end do
+    else if(tip == 3) then
+      h1 = h(1, mi)
+      h2 = h(1, ma)
+      p1 = p(1, mi)
+      do i=1,3
+        puti = p(turn3(1,i), ma)
+        putj = p(turn3(2,i), ma)
+        if(banned(puti,putj,1)) cycle
+        p2 = p(i, ma)
+        
+        hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2)
+        mat(:, min(puti, putj), max(puti, putj)) += coefs * hij
+      end do
+    else ! tip == 4
+      puti = p(1, sp)
+      putj = p(2, sp)
+      if(.not. banned(puti,putj,1)) then
+        p1 = p(1, mi)
+        p2 = p(2, mi)
+        h1 = h(1, mi)
+        h2 = h(2, mi)
+        hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2)
+        mat(:, puti, putj) += coefs * hij
+      end if
+    end if
+  end if
+end subroutine
+
+
+subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
+  use bitmasks
+  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)
+  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)
+  double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
+  double precision :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num)
+  double precision, external :: get_phase_bi, integral8
+  
+  logical :: lbanned(mo_tot_num, 2), ok
+  integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j, hfix, pfix, h1, h2, p1, p2, ib
+  
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+  
+  integer, parameter :: turn2(2) = (/2,1/)
+  integer, parameter :: turn3(2,3) = reshape((/2,3,  1,3, 1,2/), (/2,3/))
+  
+  integer :: bant
+  
+  
+  lbanned = bannedOrb
+    
+  do i=1, p(0,1)
+    lbanned(p(i,1), 1) = .true.
+  end do
+  do i=1, p(0,2)
+    lbanned(p(i,2), 2) = .true.
+  end do
+  
+  ma = 1
+  if(p(0,2) >= 2) ma = 2
+  mi = turn2(ma)
+  
+  bant = 1
+
+  if(sp == 3) then
+    !move MA
+    if(ma == 2) bant = 2
+    puti = p(1,mi)
+    hfix = h(1,ma)
+    p1 = p(1,ma)
+    p2 = p(2,ma)
+    if(.not. bannedOrb(puti, mi)) then
+      tmp_row = 0d0
+      do putj=1, hfix-1
+        if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
+        hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
+        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+      end do
+      do putj=hfix+1, mo_tot_num
+        if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
+        hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
+        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+      end do
+
+      if(ma == 1) then           
+        mat(1:N_states,1:mo_tot_num,puti) += tmp_row(1:N_states,1:mo_tot_num)
+      else
+        mat(1:N_states,puti,1:mo_tot_num) += tmp_row(1:N_states,1:mo_tot_num)
+      end if
+    end if
+
+    !MOVE MI
+    pfix = p(1,mi)
+    tmp_row = 0d0
+    tmp_row2 = 0d0
+    do puti=1,mo_tot_num
+      if(lbanned(puti,mi)) cycle
+      !p1 fixed
+      putj = p1
+      if(.not. banned(putj,puti,bant)) then
+        hij = integral8(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix)
+        tmp_row(:,puti) += hij * coefs
+      end if
+      
+      putj = p2
+      if(.not. banned(putj,puti,bant)) then
+        hij = integral8(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
+        tmp_row2(:,puti) += hij * coefs
+      end if
+    end do
+    
+    if(mi == 1) then
+      mat(:,:,p1) += tmp_row(:,:)
+      mat(:,:,p2) += tmp_row2(:,:)
+    else
+      mat(:,p1,:) += tmp_row(:,:)
+      mat(:,p2,:) += tmp_row2(:,:)
+    end if
+  else
+    if(p(0,ma) == 3) then
+      do i=1,3
+        hfix = h(1,ma)
+        puti = p(i, ma)
+        p1 = p(turn3(1,i), ma)
+        p2 = p(turn3(2,i), ma)
+        tmp_row = 0d0
+        do putj=1,hfix-1
+          if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
+          hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
+          tmp_row(:,putj) += hij * coefs
+        end do
+        do putj=hfix+1,mo_tot_num
+          if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
+          hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
+          tmp_row(:,putj) += hij * coefs
+        end do
+
+        mat(:, :puti-1, puti) += tmp_row(:,:puti-1)
+        mat(:, puti, puti:) += tmp_row(:,puti:)
+      end do
+    else
+      hfix = h(1,mi)
+      pfix = p(1,mi)
+      p1 = p(1,ma)
+      p2 = p(2,ma)
+      tmp_row = 0d0
+      tmp_row2 = 0d0
+      do puti=1,mo_tot_num
+        if(lbanned(puti,ma)) cycle
+        putj = p2
+        if(.not. banned(puti,putj,1)) then
+          hij = integral8(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1)
+          tmp_row(:,puti) += hij * coefs
+        end if
+        
+        putj = p1
+        if(.not. banned(puti,putj,1)) then
+          hij = integral8(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
+          tmp_row2(:,puti) += hij * coefs
+        end if
+      end do
+      mat(:,:p2-1,p2) += tmp_row(:,:p2-1)
+      mat(:,p2,p2:) += tmp_row(:,p2:)
+      mat(:,:p1-1,p1) += tmp_row2(:,:p1-1)
+      mat(:,p1,p1:) += tmp_row2(:,p1:)
+    end if
+  end if
+
+ !! MONO
+    if(sp == 3) then
+      s1 = 1
+      s2 = 2
+    else
+      s1 = sp
+      s2 = sp
+    end if
+
+    do i1=1,p(0,s1)
+      ib = 1
+      if(s1 == s2) ib = i1+1
+      do i2=ib,p(0,s2)
+        p1 = p(i1,s1)
+        p2 = p(i2,s2)
+        if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
+        call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
+        call i_h_j(gen, det, N_int, hij)
+        mat(:, p1, p2) += coefs * hij
+      end do
+    end do
+end subroutine
+
+
+
+
+subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, 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)
+  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)
+  double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+  
+  integer :: i, j, s, h1, h2, p1, p2, puti, putj
+  double precision :: hij, phase
+  double precision, external :: get_phase_bi, integral8
+  logical :: ok
+  
+  integer :: bant
+  bant = 1
+  
+
+  if(sp == 3) then ! AB
+    h1 = p(1,1)
+    h2 = p(1,2)
+    do p1=1, mo_tot_num
+      if(bannedOrb(p1, 1)) cycle
+      do p2=1, mo_tot_num
+        if(bannedOrb(p2,2)) cycle
+        if(banned(p1, p2, bant)) cycle ! rentable?
+        if(p1 == h1 .or. p2 == h2) then
+          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)
+        end if
+        mat(:, p1, p2) += coefs(:) * hij
+      end do
+    end do
+  else ! AA BB
+    p1 = p(1,sp)
+    p2 = p(2,sp)
+    do puti=1, mo_tot_num
+      if(bannedOrb(puti, sp)) cycle
+      do putj=puti+1, mo_tot_num
+        if(bannedOrb(putj, sp)) cycle
+        if(banned(puti, putj, bant)) cycle ! rentable?
+        if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
+          call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
+          call i_h_j(gen, det, N_int, hij)
+        else
+          hij = (integral8(p1, p2, puti, putj) -  integral8(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2)
+        end if
+        mat(:, puti, putj) += coefs(:) * hij
+      end do
+    end do
+  end if
+end subroutine
+ 
+
+subroutine past_d1(bannedOrb, p)
+  use bitmasks
+  implicit none
+
+  logical, intent(inout) :: bannedOrb(mo_tot_num, 2)
+  integer, intent(in) :: p(0:4, 2)
+  integer :: i,s
+
+  do s = 1, 2
+    do i = 1, p(0, s)
+      bannedOrb(p(i, s), s) = .true.
+    end do
+  end do
+end subroutine
+
+
+subroutine past_d2(banned, p, sp)
+  use bitmasks
+  implicit none
+
+  logical, intent(inout) :: banned(mo_tot_num, mo_tot_num)
+  integer, intent(in) :: p(0:4, 2), sp
+  integer :: i,j
+
+  if(sp == 3) then
+    do i=1,p(0,1)
+      do j=1,p(0,2)
+        banned(p(i,1), p(j,2)) = .true.
+      end do
+    end do
+  else
+    do i=1,p(0, sp)
+      do j=1,i-1
+        banned(p(j,sp), p(i,sp)) = .true.
+        banned(p(i,sp), p(j,sp)) = .true.
+      end do
+    end do
+  end if
+end subroutine
+
+
+
+subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
+  use bitmasks
+  implicit none
+  
+  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
+
+
+  integer :: i, j, na, nb, list(3)
+  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
+    do j=1, N_int
+      if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl
+      if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl
+    end do
+
+    if(interesting(i) < i_gen) then
+      fullMatch = .true.
+      return
+    end if
+
+    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))
+    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)
+    banned(list(1), list(2)) = .true.
+  end do genl
+end subroutine
+

plugins/Full_CI_ZMQ/selection_single.irp.f

@@ -0,0 +1,354 @@
+ 
+
+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
+  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)
+      delta_E = E0(istate) - Hii
+      if (delta_E < 0.d0) then
+        e_pert = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * val * val) - delta_E)
+      else
+        e_pert = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * val * val) - delta_E)
+      endif
+      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)
+  logical, intent(in) :: bannedOrb(mo_tot_num)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: vect(N_states, mo_tot_num)
+  integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2)
+  integer :: i, j, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti
+  double precision :: hij
+  double precision, external :: get_phase_bi, integral8
+  
+  integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
+  integer, parameter :: turn2(2) = (/2,1/) 
+  
+  if(h(0,sp) == 2) then
+    h1 = h(1, sp)
+    h2 = h(2, sp)
+    do i=1,3
+      puti = p(i, sp)
+      if(bannedOrb(puti)) cycle
+      p1 = p(turn3_2(1,i), sp)
+      p2 = p(turn3_2(2,i), sp)
+      hij = integral8(p1, p2, h1, h2) - integral8(p2, p1, h1, h2)
+      hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2)
+      vect(:, puti) += hij * coefs
+    end do
+  else if(h(0,sp) == 1) then
+    sfix = turn2(sp)
+    hfix = h(1,sfix)
+    pfix = p(1,sfix)
+    hmob = h(1,sp)
+    do j=1,2
+      puti = p(j, sp)
+      if(bannedOrb(puti)) cycle
+      pmob = p(turn2(j), sp)
+      hij = integral8(pfix, pmob, hfix, hmob)
+      hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix)
+      vect(:, puti) += hij * coefs
+    end do
+  else
+    puti = p(1,sp)
+    if(.not. bannedOrb(puti)) then
+      sfix = turn2(sp)
+      p1 = p(1,sfix)
+      p2 = p(2,sfix)
+      h1 = h(1,sfix)
+      h2 = h(2,sfix)
+      hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2))
+      hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2)
+      vect(:, puti) += hij * coefs
+    end if
+  end if
+end subroutine
+
+
+
+subroutine get_m1(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)
+  logical, intent(in) :: bannedOrb(mo_tot_num)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: vect(N_states, mo_tot_num)
+  integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2)
+  integer :: i, hole, p1, p2, sh
+  logical :: ok, lbanned(mo_tot_num)
+  integer(bit_kind) :: det(N_int, 2)
+  double precision :: hij
+  double precision, external :: get_phase_bi, integral8
+  
+  lbanned = bannedOrb
+  sh = 1
+  if(h(0,2) == 1) sh = 2
+  hole = h(1, sh)
+  lbanned(p(1,sp)) = .true.
+  if(p(0,sp) == 2) lbanned(p(2,sp)) = .true.
+  !print *, "SPm1", sp, sh
+  
+  p1 = p(1, sp)
+  
+  if(sp == sh) then
+    p2 = p(2, sp)
+    lbanned(p2) = .true.
+    
+    do i=1,hole-1
+      if(lbanned(i)) cycle
+      hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole))
+      hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2)
+      vect(:,i) += hij * coefs
+    end do
+    do i=hole+1,mo_tot_num
+      if(lbanned(i)) cycle
+      hij = (integral8(p1, p2, hole, i) - integral8(p2, p1, hole, i))
+      hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2)
+      vect(:,i) += hij * coefs
+    end do
+
+    call apply_particle(mask, sp, p2, det, ok,  N_int)
+    call i_h_j(gen, det, N_int, hij)
+    vect(:, p2) += hij * coefs
+  else
+    p2 = p(1, sh)
+    do i=1,mo_tot_num
+      if(lbanned(i)) cycle
+      hij = integral8(p1, p2, i, hole)
+      hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2)
+      vect(:,i) += hij * coefs
+    end do
+  end if
+
+  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
+
+
+subroutine get_m0(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)
+  logical, intent(in) :: bannedOrb(mo_tot_num)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: vect(N_states, mo_tot_num)
+  integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2)
+  integer :: i
+  logical :: ok, lbanned(mo_tot_num)
+  integer(bit_kind) :: det(N_int, 2)
+  double precision :: hij
+  
+  lbanned = bannedOrb
+  lbanned(p(1,sp)) = .true.
+  do i=1,mo_tot_num
+    if(lbanned(i)) cycle
+    call apply_particle(mask, sp, i, det, ok, N_int)
+    call i_h_j(gen, det, N_int, hij)
+    vect(:, i) += hij * coefs
+  end do
+end subroutine
+
+
+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
+
+
+

plugins/Full_CI_ZMQ/selection_slave.irp.f

@@ -0,0 +1,93 @@
+program selection_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
+!   PROVIDE ci_electronic_energy mo_tot_num N_int
+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) = 'selection'
+
+  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) == 'selection') then
+
+      ! Selection
+      ! ---------
+
+      print *,  'Selection'
+      call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states_diag)
+  
+      !$OMP PARALLEL PRIVATE(i)
+      i = omp_get_thread_num()
+      call selection_slave_tcp(i, energy)
+      !$OMP END PARALLEL
+      print *,  'Selection done'
+
+    endif
+
+  end do
+end
+
+subroutine update_energy(energy)
+  implicit none
+  double precision, intent(in) :: energy(N_states_diag)
+  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,size(ci_electronic_energy)
+      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_diag)
+  integer, intent(in)            :: i
+
+  call run_selection_slave(0,i,energy)
+end
+

plugins/Full_CI_ZMQ/selection_types.f90

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

plugins/Generators_full/generators.irp.f

@@ -30,7 +30,9 @@ END_PROVIDER
  ! Hartree-Fock determinant
  END_DOC
  integer                        :: i, k
- do i=1,N_det
+ psi_coef_generators = 0.d0
+ psi_det_generators  = 0_bit_kind
+ do i=1,N_det_generators
    do k=1,N_int
      psi_det_generators(k,1,i) = psi_det_sorted(k,1,i)
      psi_det_generators(k,2,i) = psi_det_sorted(k,2,i)