Merge branch 'master' into bugfix

* Update biblio

* Update paper

* Journal missing in research.bib

* Added paper

* Fixed Pierre Francois

INSTALL.rst

@@ -143,7 +143,7 @@ IRPF90
 to Parameters (IRP) method. 
 
 * Download the latest version of IRPF90
-  here : `<https://github.com/scemama/irpf90/releases/latest>`_ and move
+  here : `<https://gitlab.com/scemama/irpf90/-/archive/v1.7.2/irpf90-v1.7.2.tar.gz>`_ and move
   the downloaded archive in the :file:`${QP_ROOT}/external` directory
 
 * Extract the archive and go into the :file:`irpf90-*` directory to run

REPLACE

@@ -834,3 +834,4 @@ qp_name potential_sr_c_alpha_ao_pbe --rename=potential_c_alpha_ao_sr_pbe
 qp_name potential_sr_c_beta_ao_pbe --rename=potential_c_beta_ao_sr_pbe
 qp_name potential_sr_xc_alpha_ao_pbe --rename=potential_xc_alpha_ao_sr_pbe
 qp_name potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe
+qp_name disk_access_nuclear_repulsion --rename=io_nuclear_repulsion

bin/qp_tunnel

@@ -1,71 +0,0 @@
-#!/usr/bin/env python2
-
-
-"""
-Creates an ssh tunnel for using slaves on another network.
-Launch a server on the front-end node of the cluster on which the master
-process runs. Then start a client ont the front-end node of the distant
-cluster.
-
-Usage:
-      qp_tunnel server EZFIO_DIR
-      qp_tunnel client <address> EZFIO_DIR
-
-Options:
-    -h --help
-
-"""
-
-import os
-import sys
-import zmq
-
-try:
-    import qp_path
-except ImportError:
-    print "source .quantum_package.rc"
-    raise
-
-from docopt import docopt
-from ezfio import ezfio
-
-
-def get_address(filename):
-  with open(os.path.join(filename,'work','qp_run_address'),'r') as f:
-    a = f.readlines()[0].strip()
-  return a
-
-
-def set_address(filename,address):
-  with open(os.path.join(filename,'work','qp_run_address'),'r') as f:
-    backup = f.readlines()
-
-  with open(os.path.join(filename,'work','qp_run_address'),'w') as f:
-    f.write('\n'.join([address]+backup))
-
-
-def main_server(arguments,filename):
-    destination = get_address(filename)
-    print destination
-
-
-def main_client(arguments,filename):
-    destination = arguments["<address>"]
-    print destination
-
-
-def main(arguments):
-    """Main function"""
-
-    print arguments
-    filename = arguments["EZFIO_DIR"]
-
-    if arguments["server"]:
-        return main_server(arguments, filename)
-    if arguments["client"]:
-        return main_client(arguments, filename)
-
-    
-if __name__ == '__main__':
-    ARGUMENTS = docopt(__doc__)
-    main(ARGUMENTS)

configure

@@ -438,18 +438,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then
         fail
 fi
 
-BWRAP=$(find_exe bwrap)
-if [[ ${BWRAP} = $(not_found) ]] ; then
-        error "Bubblewrap (bwrap) is not installed."
-        fail
-fi
-
 LIBCAP=$(find_lib -lcap)
 if [[ ${LIBCAP} = $(not_found) ]] ; then
         error "Libcap (libcap) is not installed."
         fail
 fi
 
+BWRAP=$(find_exe bwrap)
+if [[ ${BWRAP} = $(not_found) ]] ; then
+        error "Bubblewrap (bwrap) is not installed."
+        fail
+fi
+
 OPAM=$(find_exe opam)
 if [[ ${OPAM} = $(not_found) ]] ; then
         error "OPAM (ocaml) package manager is not installed."

docs/source/research.bib

@@ -1,4 +1,14 @@
 %%% ARXIV TO BE UPDATED %%%
+@article{Hollett2019Aug,
+        author = {Hollett, Joshua W. and Loos, Pierre-Fran{\c{c}}ois},
+        title = {{Capturing static and dynamic correlation with $\Delta \text{NO}$-MP2 and $\Delta \text{NO}$-CCSD}},
+        journal = {arXiv},
+        year = {2019},
+        month = {Aug},
+        eprint = {1908.09914},
+        url = {https://arxiv.org/abs/1908.09914}
+}
+
 @article{Giner2019Jul,
         author = {Giner, Emmanuel and Scemama, Anthony and Toulouse, Julien and Loos, Pierre-Fran{\c{c}}ois},
         title = {{Chemically Accurate Excitation Energies With Small Basis Sets}},

docs/source/users_guide/qp_tunnel.rst

@@ -0,0 +1,89 @@
+.. _qp_tunnel:
+
+=========
+qp_tunnel
+=========
+
+.. TODO
+
+.. program:: qp_tunnel
+
+Establishes a tunnel to allow communications between machines within
+different networks, for example multiple MPI slave jobs running on
+different clusters.
+
+
+Usage
+-----
+
+.. code:: bash
+
+  qp_tunnel [-g] (ADDRESS|EZFIO_DIR)
+
+``EZFIO_DIR`` is the name of the |EZFIO| directory containing the data,
+and ``ADDRESS`` is the address of another tunnel.
+
+
+.. option:: -h, --help
+
+   Displays the help message
+
+
+.. option:: -g, --get-input
+
+   Download the EZFIO directory from the remote instance of qp_tunnel.
+
+
+Example
+-------
+
+.. code:: text
+
+    +-------------------+    +------------------+
+    |                   |    |                  |
+    |   N1_1 N1_2 N1_3  |    |  N2_1 N2_2 N2_3  |
+    |    |    |    |    |    |   |    |    |    |
+    |    +----+----+    |    |   +----+----+    |
+    |         |         |    |        |         |
+    | C1      F1        |    |        F2     C2 |
+    |         +---------=----=--------+         |
+    |                   |    |                  |
+    +-------------------+    +------------------+
+
+
+Imagine you have two clusters, C1 and C2. Each cluster is accessible via SSH
+on a front-end named respectively F1 and F2. Groups of nodes N1 and N2 have
+been reserved by the batch scheduling system on both clusters.
+Each node in N1 is on the same network as the other nodes of N1, but they
+can't access the network on which the nodes of N2 are.
+
+1) Start a parallel simulation on the cluster C1, running on nodes N1.
+   We assume that there is a shared file system, such that F1 can access
+   the EZFIO directory. We also assume that F1 can communicate with the
+   nodes of N1.
+
+2) Run a tunnel on the front-end F1 and keep it running:
+
+.. code:: bash
+
+    me@f1 $ qp_tunnel my_directory.ezfio
+    Connect to:
+    tcp://31.122.230.47:42379
+    Ready   
+
+3) On the front-end F2, run another instance connecting to the other one,
+   which will fetch the |EZFIO| directory:
+
+.. code:: bash
+
+    me@f2 $ qp_tunnel --get-input tcp://31.122.230.47:42379
+    Connect to:
+    tcp://31.122.209.139:42379
+    Communication [ OK ]
+    Getting input... my_directory.ezfio ...done
+    Ready
+
+4) Keep the tunnel running, and you can now run a slave simulation within the
+   nodes N2.
+
+

man/qp_tunnel.1

@@ -0,0 +1,139 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_TUNNEL" "1" "Jun 15, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_tunnel \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+Establishes a tunnel to allow communications between machines within
+different networks, for example multiple MPI slave jobs running on
+different clusters.
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_tunnel [\-g] (ADDRESS|EZFIO_DIR)
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+\fBEZFIO_DIR\fP is the name of the \fI\%EZFIO\fP directory containing the data,
+and \fBADDRESS\fP is the address of another tunnel.
+.INDENT 0.0
+.TP
+.B \-h, \-\-help
+Displays the help message
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-g, \-\-get\-input
+Download the EZFIO directory from the remote instance of qp_tunnel.
+.UNINDENT
+.SH EXAMPLE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
++\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+    +\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+
+|                   |    |                  |
+|   N1_1 N1_2 N1_3  |    |  N2_1 N2_2 N2_3  |
+|    |    |    |    |    |   |    |    |    |
+|    +\-\-\-\-+\-\-\-\-+    |    |   +\-\-\-\-+\-\-\-\-+    |
+|         |         |    |        |         |
+| C1      F1        |    |        F2     C2 |
+|         +\-\-\-\-\-\-\-\-\-=\-\-\-\-=\-\-\-\-\-\-\-\-+         |
+|                   |    |                  |
++\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+    +\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+Imagine you have two clusters, C1 and C2. Each cluster is accessible via SSH
+on a front\-end named respectively F1 and F2. Groups of nodes N1 and N2 have
+been reserved by the batch scheduling system on both clusters.
+Each node in N1 is on the same network as the other nodes of N1, but they
+can’t access the network on which the nodes of N2 are.
+.INDENT 0.0
+.IP 1. 3
+Start a parallel simulation on the cluster C1, running on nodes N1.
+We assume that there is a shared file system, such that F1 can access
+the EZFIO directory. We also assume that F1 can communicate with the
+nodes of N1.
+.IP 2. 3
+Run a tunnel on the front\-end F1 and keep it running:
+.UNINDENT
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+me@f1 $ qp_tunnel my_directory.ezfio
+Connect to:
+tcp://31.122.230.47:42379
+Ready
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.IP 3. 3
+On the front\-end F2, run another instance connecting to the other one,
+which will fetch the \fI\%EZFIO\fP directory:
+.UNINDENT
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+me@f2 $ qp_tunnel \-\-get\-input tcp://31.122.230.47:42379
+Connect to:
+tcp://31.122.209.139:42379
+Communication [ OK ]
+Getting input... my_directory.ezfio ...done
+Ready
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.IP 4. 3
+Keep the tunnel running, and you can now run a slave simulation within the
+nodes N2.
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2019, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.

man/rotate_mos.1

@@ -34,7 +34,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
 .INDENT 3.5
 Rotates molecular orbitals i and j by combining them as
 $1/sqrt{2} ( phi_i + phi_j )$ and
-$1/sqrt{2} ( phi_i \- phi_j )$.
+$1/sqrt{2} ( phi_i - phi_j )$.
 .sp
 Needs:
 .INDENT 0.0

man/test.1

@@ -0,0 +1,55 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "TEST" "1" "Jun 15, 2019" "2.0" "Quantum Package"
+.SH NAME
+test \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.INDENT 0.0
+.INDENT 3.5
+Calls:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBtwo_e_integrals_index()\fP
+.UNINDENT
+.INDENT 2.0
+.IP \(bu 2
+\fBtwo_e_integrals_index_reverse()\fP
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2019, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.

ocaml/.merlin

@@ -1,4 +1,4 @@
-PKG core ZMQ cryptokit
+PKG core zmq cryptokit
 B _build/
 
 

ocaml/Bitlist.ml

@@ -7,82 +7,61 @@ Type for bits strings
 list of Bits
 *)
 
-type t = Bit.t list
+type t = int64 array
+
+
+let n_int = Array.length
+
+(* Create a zero bit list *)
+let zero n_int =
+  Array.make (N_int_number.to_int n_int) 0L
 
 (* String representation *)
 let to_string b =
-  let rec do_work accu = function
+  let int64_to_string x = 
-    | [] -> accu
+    String.init 64 (fun i ->
-    | head :: tail ->
+      if Int64.logand x @@ Int64.shift_left 1L i <> 0L then
-      let new_accu = (Bit.to_string head) ^ accu
+        '+'
-      in do_work new_accu tail
+      else
+        '-')
   in
-  do_work "" b
+  Array.map int64_to_string b
+  |> Array.to_list
+  |> String.concat ""
 
 
 let of_string ?(zero='0') ?(one='1') s =
-  List.init (String.length s) (String.get s)
+  let n_int = ( (String.length s - 1) lsr 6 ) + 1 in
-  |> List.rev_map ( fun c ->
+  let result = Array.make n_int 0L in
-     if (c = zero) then Bit.Zero
+  String.iteri (fun i c ->
-     else if (c = one) then Bit.One
+    if c = one then
-     else (failwith ("Error in bitstring ") ) )
+      begin
+        let iint = i lsr 6 in (* i / 64 *)
+        let k = i - (iint lsl 6) in 
+        result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k;
+      end) s;
+  result
 
-let of_string_mp s =
+let of_string_mp = of_string ~zero:'-' ~one:'+'
-  List.init (String.length s) (String.get s)
-  |> List.rev_map (function
-    | '-' -> Bit.Zero
-    | '+' -> Bit.One
-    | _   -> failwith ("Error in bitstring ") )
 
 
 (* Create a bit list from an int64 *)
-let of_int64 i =
+let of_int64 i = [| i |]
-
-  let rec do_work accu = function
-  | 0L -> Bit.Zero :: accu |> List.rev
-  | 1L -> Bit.One  :: accu |> List.rev
-  | i  ->
-    let b =
-      match (Int64.logand i 1L ) with
-      | 0L -> Bit.Zero
-      | 1L -> Bit.One
-      | _ -> raise (Failure "i land 1 not in (0,1)")
-    in
-    do_work (b :: accu) (Int64.shift_right_logical i 1)
-  in
-
-  let adjust_length result =
-    let rec do_work accu = function
-    | 64 -> List.rev accu
-    | i when i>64 -> raise (Failure "Error in of_int64 > 64")
-    | i when i<0 -> raise (Failure "Error in of_int64 < 0")
-    | i -> do_work (Bit.Zero :: accu)  (i+1)
-    in
-    do_work (List.rev result) (List.length result)
-  in
-  adjust_length (do_work [] i)
-
 
 (* Create an int64 from a bit list *)
-let to_int64 l =
+let to_int64 = function
-  assert ( (List.length l) <= 64) ;
+  | [| i |] -> i
-  let rec do_work accu = function
+  | _ -> failwith "N_int > 1"
-    | [] -> accu
+
-    | Bit.Zero::tail -> do_work Int64.(shift_left accu 1) tail
+
-    | Bit.One::tail  -> do_work Int64.(logor one (shift_left accu 1)) tail
+(* Create a bit list from an array of int64 *)
-  in do_work Int64.zero (List.rev l)
+external of_int64_array : int64 array -> t = "%identity"
+external to_int64_array : t -> int64 array = "%identity"
 
  
 (* Create a bit list from a list of int64 *)
 let of_int64_list l =
-  List.map of_int64 l
+  Array.of_list l |> of_int64_array
-  |> List.concat
-
-(* Create a bit list from an array of int64 *)
-let of_int64_array l =
-  Array.map of_int64 l
-  |> Array.to_list
-  |> List.concat
 
 
 (* Compute n_int *)
@@ -91,101 +70,64 @@ let n_int_of_mo_num mo_num =
   N_int_number.of_int ( (mo_num-1)/bit_kind_size + 1 )
 
 
-(* Create a zero bit list *)
-let zero n_int =
-  let n_int = N_int_number.to_int n_int in
-  let a = Array.init n_int (fun i-> 0L)  in
-  of_int64_list ( Array.to_list a )
-
 
 (* Create an int64 list from a bit list *)
 let to_int64_list l =
-  let rec do_work accu buf counter = function
+  to_int64_array l |> Array.to_list
-    | [] ->
-        begin
-          match buf with
-          | [] -> accu
-          | _  -> (List.rev buf)::accu
-        end
-    | i::tail ->
-        if (counter < 64) then
-          do_work accu (i::buf) (counter+1) tail
-        else
-          do_work ( (List.rev (i::buf))::accu) [] 1 tail
-  in
-  let l = do_work [] [] 1 l
-  in
-  List.rev_map to_int64 l
 
-(* Create an array of int64 from a bit list *)
-let to_int64_array l =
-    to_int64_list l
-    |> Array.of_list
 
 (* Create a bit list from a list of MO indices *)
 let of_mo_number_list n_int l =
-  let n_int = N_int_number.to_int n_int in
+  let result = zero n_int in
-  let length = n_int*64 in
+  List.iter (fun j ->
-  let a = Array.make length (Bit.Zero) in
+        let i = (MO_number.to_int j) - 1 in
-  List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l;
+        let iint = i lsr 6 in (* i / 64 *)
-  Array.to_list a
+        let k = i - (iint lsl 6) in
+        result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k;
+      ) l;
+  result
 
 
 let to_mo_number_list l =
-  let a = Array.of_list l in
+  let rec aux_one x shift accu = function
-  let mo_num = MO_number.get_max () in
+    | -1 -> accu
-  let rec do_work accu = function
+    | i  -> if Int64.logand x (Int64.shift_left 1L i) <> 0L then
-  | 0 -> accu
+              aux_one x shift ( (i+shift) ::accu) (i-1)
-  | i ->
+            else
-      begin
+              aux_one x shift accu (i-1)
-        let new_accu =
-        match a.(i-1) with
-        | Bit.One  -> (MO_number.of_int ~max:mo_num i)::accu
-        | Bit.Zero -> accu
-        in
-        do_work new_accu (i-1)
-      end
   in
-  do_work [] (List.length l)
+  Array.mapi (fun i x ->
-
+    let shift = (i lsr 6) lsl 6 + 1 in
+    aux_one x shift [] 63
+    ) l
+  |> Array.to_list
+  |> List.concat
+  |> List.map MO_number.of_int
 
 
 
 (* logical operations on bit_list *)
-let logical_operator2 op a b =
+let and_operator a b = Array.map2 Int64.logand a b
-  let rec do_work_binary result a b =
+let xor_operator a b = Array.map2 Int64.logxor a b
-  match a, b with
+let  or_operator a b = Array.map2 Int64.logor  a b
-  | [], [] -> result
+let not_operator   b = Array.map  Int64.lognot   b
-  | [], _  | _ , [] -> raise (Failure "Lists should have same length")
+
-  | (ha::ta), (hb::tb) ->
+
-    let newbit = op ha hb
+
-    in  do_work_binary (newbit::result) ta tb
+let pop_sign = 
+  let mask =
+      (Int64.pred (Int64.shift_left 1L 63))
   in
-  List.rev (do_work_binary [] a b)
+  fun x -> Int64.logand mask x
-
-
-let logical_operator1 op b =
-  let rec do_work_unary result b =
-  match b with
-  | [] -> result
-  | (hb::tb) ->
-    let newbit = op hb
-    in  do_work_unary (newbit::result) tb
-  in
-  List.rev (do_work_unary [] b)
-
-
-let and_operator a b = logical_operator2 Bit.and_operator a b
-let xor_operator a b = logical_operator2 Bit.xor_operator a b
-let  or_operator a b = logical_operator2  Bit.or_operator a b
-let not_operator   b = logical_operator1 Bit.not_operator   b
 
 
 let popcnt b =
-  List.fold_left (fun accu -> function
+  Array.fold_left (fun accu x ->
-    | Bit.One -> accu+1
+    if x >= 0L then
-    | Bit.Zero -> accu
+      accu + (Z.popcount @@ Z.of_int64 x)
-  ) 0 b
+    else
+      accu + 1 + (Z.popcount @@ Z.of_int64 (pop_sign x))
+      ) 0 b
 
 
 

ocaml/Bitlist.mli

@@ -1,4 +1,4 @@
-type t = Bit.t list
+type t
 
 (** The zero bit list *)
 val zero : Qptypes.N_int_number.t -> t

ocaml/Determinant.ml

@@ -25,19 +25,6 @@ let to_bitlist_couple x =
   in (xa,xb)
 
 
-let bitlist_to_string ~mo_num x =
-  let len =
-    MO_number.to_int mo_num
-  in
-  let s =
-    List.map (function
-        | Bit.Zero -> "-"
-        | Bit.One  -> "+"
-    ) x
-    |> String.concat ""
-  in
-  String.sub s 0 len
-
 
 
 let of_int64_array ~n_int ~alpha ~beta x =
@@ -48,37 +35,29 @@ let of_int64_array ~n_int ~alpha ~beta x =
    in
    if ( (Bitlist.popcnt a) <> alpha) then
      begin
-       let mo_num = MO_number.get_max () in
-       let mo_num = MO_number.of_int mo_num ~max:mo_num  in
        failwith (Printf.sprintf "Expected %d electrons in alpha determinant
-%s" alpha (bitlist_to_string ~mo_num:mo_num a) )
+%s" alpha (Bitlist.to_string a) )
      end;
    if ( (Bitlist.popcnt b) <> beta ) then
      begin
-       let mo_num = MO_number.get_max () in
-       let mo_num = MO_number.of_int mo_num ~max:mo_num  in
        failwith (Printf.sprintf "Expected %d electrons in beta determinant
-%s" beta (bitlist_to_string ~mo_num:mo_num b) )
+%s" beta (Bitlist.to_string b) )
      end;
    x
 
 
-let of_bitlist_couple ?n_int ~alpha ~beta (xa,xb) =
+let of_bitlist_couple ~n_int ~alpha ~beta (xa,xb) =
   let ba, bb =
     Bitlist.to_int64_array xa ,
     Bitlist.to_int64_array xb
-  and n_int =
-    match n_int with
-    | Some x -> x
-    | None -> Bitlist.n_int_of_mo_num (List.length xa)
   in
   of_int64_array ~n_int ~alpha ~beta (Array.concat [ba;bb])
 
 
 let to_string ~mo_num x =
   let (xa,xb) = to_bitlist_couple x in
-  [ "  " ; bitlist_to_string ~mo_num xa ; "\n" ;
+  [ "  " ; Bitlist.to_string xa ; "\n" ;
-    "  " ; bitlist_to_string ~mo_num xb ]
+    "  " ; Bitlist.to_string xb ]
   |> String.concat ""
 
 

ocaml/Determinant.mli

@@ -24,7 +24,7 @@ val to_alpha_beta : t -> (int64 array)*(int64 array)
 val to_bitlist_couple : t -> Bitlist.t * Bitlist.t
 
 (** Create from a bit list *)
-val of_bitlist_couple : ?n_int:Qptypes.N_int_number.t ->
+val of_bitlist_couple : n_int:Qptypes.N_int_number.t ->
   alpha:Qptypes.Elec_alpha_number.t ->
   beta:Qptypes.Elec_beta_number.t ->
   Bitlist.t * Bitlist.t -> t

ocaml/Input_determinants_by_hand.ml

@@ -7,14 +7,14 @@ module Determinants_by_hand : sig
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
+      n_det_qp_edit          : Det_number.t;
       n_states               : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
       state_average_weight   : Positive_float.t array;
     } [@@deriving sexp]
-  val read  : unit -> t
+  val read : ?full:bool -> unit -> t option
-  val read_maybe  : unit -> t option
   val write : t -> unit
   val to_string : t -> string
   val to_rst : t -> Rst_string.t
@@ -28,6 +28,7 @@ end = struct
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
+      n_det_qp_edit          : Det_number.t;
       n_states               : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
@@ -38,8 +39,6 @@ end = struct
 
   let get_default = Qpackage.get_ezfio_default "determinants";;
 
-  let n_det_read_max = 10_000 ;;
-
   let read_n_int () =
     if not (Ezfio.has_determinants_n_int()) then
        Ezfio.get_mo_basis_mo_num ()
@@ -80,11 +79,27 @@ end = struct
     |> Det_number.of_int
   ;;
 
+  let read_n_det_qp_edit () =
+    if not (Ezfio.has_determinants_n_det_qp_edit ()) then
+      begin
+        let n_det = read_n_det () |> Det_number.to_int in
+        Ezfio.set_determinants_n_det_qp_edit n_det
+      end;
+    Ezfio.get_determinants_n_det_qp_edit ()
+    |> Det_number.of_int
+  ;;
+
   let write_n_det n =
     Det_number.to_int n
     |> Ezfio.set_determinants_n_det
   ;;
 
+  let write_n_det_qp_edit n =
+    let n_det = read_n_det () |> Det_number.to_int in
+    min n_det (Det_number.to_int n)
+    |> Ezfio.set_determinants_n_det_qp_edit
+  ;;
+
   let read_n_states () =
     if not (Ezfio.has_determinants_n_states ()) then
       Ezfio.set_determinants_n_states 1
@@ -178,7 +193,7 @@ end = struct
     |> Ezfio.set_determinants_expected_s2
   ;;
 
-  let read_psi_coef () =
+  let read_psi_coef ~read_only () =
     if not (Ezfio.has_determinants_psi_coef ()) then
       begin
         let n_states =
@@ -189,7 +204,12 @@ end = struct
           ~data:(List.init n_states (fun i -> if (i=0) then 1. else 0. ))
           |> Ezfio.set_determinants_psi_coef
       end;
-    Ezfio.get_determinants_psi_coef ()
+    begin
+      if read_only then
+        Ezfio.get_determinants_psi_coef_qp_edit ()
+      else
+        Ezfio.get_determinants_psi_coef ()
+    end
     |> Ezfio.flattened_ezfio
     |> Array.map Det_coef.of_float
   ;;
@@ -202,12 +222,15 @@ end = struct
     and n_states =
       States_number.to_int n_states
     in
-    Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
+    let r = 
-    |> Ezfio.set_determinants_psi_coef
+      Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
+    in
+    Ezfio.set_determinants_psi_coef r;
+    Ezfio.set_determinants_psi_coef_qp_edit r
   ;;
 
 
-  let read_psi_det () =
+  let read_psi_det ~read_only () =
     let n_int = read_n_int ()
     and n_alpha = Ezfio.get_electrons_elec_alpha_num ()
         |> Elec_alpha_number.of_int
@@ -233,13 +256,21 @@ end = struct
           |> Ezfio.set_determinants_psi_det ;
       end  ;
     let n_int = N_int_number.to_int n_int in
-    let psi_det_array = Ezfio.get_determinants_psi_det () in
+    let psi_det_array =
+      if read_only then
+        Ezfio.get_determinants_psi_det_qp_edit ()
+      else
+        Ezfio.get_determinants_psi_det ()
+    in
     let dim = psi_det_array.Ezfio.dim
     and data =  Ezfio.flattened_ezfio psi_det_array
     in
     assert (n_int = dim.(0));
     assert (dim.(1) = 2);
-    assert (dim.(2) = (Det_number.to_int (read_n_det ())));
+    if read_only then
+      assert (dim.(2) = (Det_number.to_int (read_n_det_qp_edit ())))
+    else
+      assert (dim.(2) = (Det_number.to_int (read_n_det ())));
     Array.init dim.(2) (fun i ->
       Array.sub data (2*n_int*i) (2*n_int) )
     |> Array.map (Determinant.of_int64_array
@@ -252,40 +283,45 @@ end = struct
       |> Array.concat
       |> Array.to_list
     in
-    Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data
+    let r = 
-    |> Ezfio.set_determinants_psi_det
+      Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data
+    in
+    Ezfio.set_determinants_psi_det r;
+    Ezfio.set_determinants_psi_det_qp_edit r
   ;;
 
 
-  let read () =
+  let read ?(full=true) () =
+
+    let n_det_qp_edit = read_n_det_qp_edit () in
+    let n_det         = read_n_det         () in
+    let read_only = 
+      if full then false else n_det_qp_edit <> n_det
+    in
+
     if (Ezfio.has_mo_basis_mo_num ()) then
+      try
+        Some
         { n_int                  = read_n_int ()                ;
           bit_kind               = read_bit_kind ()             ;
           n_det                  = read_n_det ()                ;
+          n_det_qp_edit          = read_n_det_qp_edit ()        ;
           expected_s2            = read_expected_s2 ()          ;
-          psi_coef               = read_psi_coef ()             ;
+          psi_coef               = read_psi_coef ~read_only ()  ;
-          psi_det                = read_psi_det ()              ;
+          psi_det                = read_psi_det ~read_only  ()  ;
           n_states               = read_n_states ()             ;
           state_average_weight   = read_state_average_weight () ;
         }
+      with _ -> None
     else
-      failwith "No molecular orbitals, so no determinants"
+      (* No molecular orbitals, so no determinants *)
-  ;;
-
-  let read_maybe () =
-    let n_det =
-       read_n_det ()
-    in
-    if ( (Det_number.to_int n_det) < n_det_read_max ) then
-      try Some (read ()) with
-      | Failure _ -> None
-    else
       None
   ;;
 
   let write { n_int                ;
               bit_kind             ;
               n_det                ;
+              n_det_qp_edit        ;
               expected_s2          ;
               psi_coef             ;
               psi_det              ;
@@ -297,9 +333,13 @@ end = struct
      write_n_det n_det;
      write_n_states n_states;
      write_expected_s2 expected_s2;
-     write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
+     if n_det <= n_det_qp_edit then
-     write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
+        begin
-     write_state_average_weight state_average_weight;
+          write_n_det_qp_edit n_det;
+          write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
+          write_psi_det ~n_int:n_int ~n_det:n_det psi_det
+        end;
+     write_state_average_weight state_average_weight
   ;;
 
 
@@ -316,11 +356,13 @@ end = struct
         |> States_number.to_int
       and ndet =
         Det_number.to_int b.n_det
+      and ndet_qp_edit =
+        Det_number.to_int b.n_det_qp_edit
       in
       let coefs_string i =
         Array.init nstates (fun j ->
           let ishift =
-            j*ndet
+            j*ndet_qp_edit
           in
           if (ishift < Array.length b.psi_coef) then
             b.psi_coef.(i+ishift)
@@ -331,7 +373,7 @@ end = struct
         )
         |> Array.to_list |> String.concat "\t"
       in
-      Array.init ndet (fun i ->
+      Array.init ndet_qp_edit (fun i ->
         Printf.sprintf "  %s\n%s\n"
           (coefs_string i)
           (Determinant.to_string ~mo_num:mo_num b.psi_det.(i)
@@ -472,6 +514,7 @@ psi_det                = %s
 
     (* Handle determinants *)
     let psi_det =
+      let n_int = N_int_number.of_int @@ (MO_number.get_max () - 1) / 64 + 1 in
       let n_alpha = Ezfio.get_electrons_elec_alpha_num ()
         |> Elec_alpha_number.of_int
       and n_beta = Ezfio.get_electrons_elec_beta_num ()
@@ -483,8 +526,8 @@ psi_det                = %s
           begin
             let newdet =
                (Bitlist.of_string ~zero:'-' ~one:'+' alpha ,
-               Bitlist.of_string ~zero:'-' ~one:'+' beta)
+                Bitlist.of_string ~zero:'-' ~one:'+' beta)
-               |> Determinant.of_bitlist_couple  ~alpha:n_alpha ~beta:n_beta
+               |> Determinant.of_bitlist_couple ~n_int ~alpha:n_alpha ~beta:n_beta
                |> Determinant.sexp_of_t
                |> Sexplib.Sexp.to_string
             in
@@ -492,9 +535,6 @@ psi_det                = %s
           end
       | _::tail -> read_dets accu tail
       in
-      let dets =
-        List.map String_ext.rev dets
-      in
       let a =
         read_dets [] dets
         |> String.concat ""
@@ -510,9 +550,11 @@ psi_det                = %s
       Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
     and n_states =
       Printf.sprintf "(n_states %d)" (States_number.to_int @@ read_n_states ())
+    and n_det_qp_edit =
+      Printf.sprintf "(n_det_qp_edit %d)" (Det_number.to_int @@ read_n_det_qp_edit ())
     in
     let s =
-       String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det]
+       String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det ; n_det_qp_edit ]
     in
 
 
@@ -527,7 +569,9 @@ psi_det                = %s
       Det_number.to_int n_det_new
     in
     let det =
-      read ()
+      match read () with
+      | Some x -> x
+      | None   -> failwith "No determinants in file"
     in
     let n_det_old, n_states =
       Det_number.to_int det.n_det,
@@ -558,7 +602,9 @@ psi_det                = %s
   let extract_state istate =
     Printf.printf "Extracting state %d\n" (States_number.to_int istate);
     let det =
-      read ()
+      match read () with
+      | Some x -> x
+      | None   -> failwith "No determinants in file"
     in
     let n_det, n_states =
       Det_number.to_int det.n_det,
@@ -588,7 +634,9 @@ psi_det                = %s
   let extract_states range =
     Printf.printf "Extracting states %s\n" (Range.to_string range);
     let det =
-      read ()
+      match read () with
+      | Some x -> x
+      | None   -> failwith "No determinants in file"
     in
     let n_det, n_states =
       Det_number.to_int det.n_det,
@@ -614,7 +662,8 @@ psi_det                = %s
             j*n_det
           in
           for i=0 to (n_det-1) do
-            det.psi_coef.(!state_shift+i) <- det.psi_coef.(i+ishift)
+            det.psi_coef.(!state_shift+i) <-
+            det.psi_coef.(i+ishift)
           done
       end;
       state_shift := !state_shift + n_det

ocaml/Input_mo_basis.ml

@@ -65,8 +65,15 @@ end = struct
 
 
   let read_mo_num () =
-    Ezfio.get_mo_basis_mo_num ()
+    let elec_alpha_num =
-    |> MO_number.of_int
+      Ezfio.get_electrons_elec_alpha_num ()
+    in
+    let result = 
+      Ezfio.get_mo_basis_mo_num ()
+    in
+    if result < elec_alpha_num then
+      failwith "More alpha electrons than MOs";
+    MO_number.of_int result
 
 
   let read_mo_class () =

ocaml/Makefile

@@ -43,7 +43,7 @@ $(QP_ROOT)/data/executables: remake_executables element_create_db.byte Qptypes.m
 	$(QP_ROOT)/ocaml/element_create_db.byte
 
 external_libs:
-	opam install cryptokit core
+	opam install cryptokit sexplib
 
 qpackage.odocl: $(MLIFILES)
 	ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl

ocaml/_tags

@@ -1,4 +1,4 @@
-true: package(cryptokit,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
+true: package(cryptokit,zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
 true: thread
 false: profile
 <*byte>   : linkdep(c_bindings.o), custom

ocaml/qp_print_basis.ml

@@ -44,9 +44,12 @@ let psi_det () =
   let psi_det =
     Input.Determinants_by_hand.read ()
   in
-  Input.Determinants_by_hand.to_rst psi_det
+  match psi_det with
-  |> Rst_string.to_string
+  | Some psi_det -> 
-  |> print_endline
+      Input.Determinants_by_hand.to_rst psi_det
+      |> Rst_string.to_string
+      |> print_endline
+  | None -> ()
 
 
 

ocaml/qp_set_mo_class.ml

@@ -112,9 +112,9 @@ let set ~core ~inact ~act ~virt ~del =
   and av = Excitation.create_single act virt
   in
   let single_excitations = [ ia ; aa ; av ]
-                           |> List.map (fun x ->
+                           |> List.map (fun z ->
                                let open Excitation in
-                               match x with
+                               match z with
                                | Single (x,y) ->
                                  ( MO_class.to_bitlist n_int (Hole.to_mo_class x),
                                    MO_class.to_bitlist n_int (Particle.to_mo_class y) )
@@ -187,9 +187,10 @@ let set ~core ~inact ~act ~virt ~del =
     match aa with
     | Double _ -> assert false
     | Single (x,y) ->
-      ( MO_class.to_bitlist n_int (Hole.to_mo_class x) ) @
+      Bitlist.to_int64_list
-      ( MO_class.to_bitlist n_int (Particle.to_mo_class y) )
+        ( MO_class.to_bitlist n_int (    Hole.to_mo_class x) )  @
-      |> Bitlist.to_int64_list
+      Bitlist.to_int64_list
+        ( MO_class.to_bitlist n_int (Particle.to_mo_class y) )
   in
   Ezfio.set_bitmasks_n_mask_cas 1;
   Ezfio.ezfio_array_of_list ~rank:3 ~dim:([| (N_int_number.to_int n_int) ; 2; 1|]) ~data:result

ocaml/qp_tunnel.ml

@@ -0,0 +1,470 @@
+open Qputils
+open Qptypes
+
+type ezfio_or_address = EZFIO of string | ADDRESS of string
+type req_or_sub = REQ | SUB 
+
+let localport = 42379
+
+
+let in_time_sum = ref 1.e-9
+and in_size_sum = ref 0.
+
+
+let () =
+  let open Command_line in
+  begin
+    "Creates an ssh tunnel for using slaves on another network. Launch a server on the front-end node of the cluster on which the master process runs. Then start a client ont the front-end node of the distant cluster."
+
+    |> set_footer_doc ;
+
+    [ { short='g' ; long="get-input" ; opt=Optional ;
+        doc="Downloads the EZFIO directory." ;
+        arg=Without_arg; } ;
+
+      { short='v' ; long="verbose" ; opt=Optional ;
+        doc="Prints the transfer speed." ;
+        arg=Without_arg; } ;
+
+      anonymous
+        "(EZFIO_DIR|ADDRESS)"
+         Mandatory
+        "EZFIO directory or address.";
+    ] |> set_specs
+  end;
+
+  let arg =
+    let x = 
+      match Command_line.anon_args () with
+      | [x] -> x
+      | _ -> begin
+               Command_line.help () ;
+               failwith "EZFIO_FILE or ADDRESS is missing"
+             end
+    in
+    if Sys.file_exists x && Sys.is_directory x then
+      EZFIO x
+    else
+      ADDRESS x
+  in
+
+  let verbose =
+    Command_line.get_bool "verbose"
+  in
+
+
+
+  let localhost = 
+     Lazy.force TaskServer.ip_address
+  in
+
+
+  let long_address  =
+    match arg with
+    | ADDRESS x -> x
+    | EZFIO   x -> 
+	let ic = 
+	  Filename.concat  (Qpackage.ezfio_work x)  "qp_run_address"
+          |> open_in
+	in
+        let result = 
+          input_line ic
+          |> String.trim
+        in
+        close_in ic;
+        result
+  in
+  
+  let protocol, address, port =
+    match String.split_on_char ':' long_address with
+    | t :: a :: p :: [] -> t, a, int_of_string p
+    | _ -> failwith @@ 
+           Printf.sprintf "%s : Malformed address" long_address
+  in
+
+
+  let zmq_context = 
+    Zmq.Context.create ()
+  in
+
+
+  (** Check availability of the ports *)
+  let localport =
+    let dummy_socket =
+      Zmq.Socket.create zmq_context Zmq.Socket.rep
+    in
+    let rec try_new_port port_number =
+      try
+        List.iter (fun i ->
+            let address =
+              Printf.sprintf "tcp://%s:%d" localhost (port_number+i)
+            in
+            Zmq.Socket.bind    dummy_socket  address;
+            Zmq.Socket.unbind  dummy_socket  address
+        ) [ 0;1;2;3;4;5;6;7;8;9 ] ;
+        port_number
+      with
+      | Unix.Unix_error _ -> try_new_port (port_number+100)
+    in
+    let result =
+      try_new_port localport
+    in
+    Zmq.Socket.close dummy_socket;
+    result
+  in
+
+
+  let create_socket  sock_type  bind_or_connect  addr =
+    let socket = 
+      Zmq.Socket.create zmq_context sock_type
+    in
+    let () = 
+      try
+        bind_or_connect socket addr
+      with
+      | _ -> failwith @@
+             Printf.sprintf "Unable to establish connection to %s." addr
+    in
+    socket
+  in
+
+
+  (* Handle termination *)
+  let run_status = ref true in
+  let handler =
+    Sys.Signal_handle (fun signum ->
+        run_status := false;
+        Sys.set_signal  signum  Sys.Signal_default
+      )
+  in
+  Sys.set_signal  Sys.sigusr1  handler;
+  Sys.set_signal  Sys.sigint   handler;
+
+
+  let new_thread req_or_sub addr_in addr_out =
+    let socket_in, socket_out =
+      match req_or_sub with
+      | REQ ->
+          create_socket  Zmq.Socket.router  Zmq.Socket.bind     addr_in,
+          create_socket  Zmq.Socket.dealer  Zmq.Socket.connect  addr_out
+      | SUB ->
+          create_socket  Zmq.Socket.sub  Zmq.Socket.connect  addr_in,
+          create_socket  Zmq.Socket.pub  Zmq.Socket.bind     addr_out
+    in
+
+    if req_or_sub = SUB then
+        Zmq.Socket.subscribe  socket_in  "";
+
+
+    (*
+    let action = 
+      if verbose then
+        begin
+          match req_or_sub with
+          | REQ -> (fun () ->
+                    let msg = 
+                      Zmq.Socket.recv_all  socket_in
+                    in
+                    let t0 = Unix.gettimeofday () in
+                    Zmq.Socket.send_all socket_out  msg;
+                    let in_size =
+                      float_of_int ( List.fold_left (fun accu x -> accu + String.length x) 0 msg )
+                      /. 8192. /. 1024.
+                    in
+                    let msg = 
+                      Zmq.Socket.recv_all  socket_out
+                    in
+                    let t1 = Unix.gettimeofday () in
+                    Zmq.Socket.send_all  socket_in msg;
+                    let in_time = t1 -. t0 in
+                    in_time_sum := !in_time_sum +. in_time;
+                    in_size_sum := !in_size_sum +. in_size;
+                    Printf.printf " %16.2f MiB/s  --  %16.2f MiB/s\n%!" (in_size /. in_time) (!in_size_sum /. !in_time_sum);
+                    )
+          | SUB -> (fun () ->
+                    Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out)
+        end
+      else
+        begin
+          match req_or_sub with
+          | REQ -> (fun () ->
+                    Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out;
+                    Zmq.Socket.recv_all  socket_out |> Zmq.Socket.send_all  socket_in )
+          | SUB -> (fun () ->
+                    Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out)
+        end
+    in
+    *)
+
+   let action_in = 
+      match req_or_sub with
+        | REQ -> (fun () -> Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out)
+        | SUB -> (fun () -> Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out)
+    in
+
+    let action_out = 
+      match req_or_sub with
+        | REQ -> (fun () -> Zmq.Socket.recv_all  socket_out |> Zmq.Socket.send_all  socket_in )
+        | SUB -> (fun () -> () )
+    in
+
+    let pollitem =
+      Zmq.Poll.mask_of
+        [| (socket_in,  Zmq.Poll.In) ; (socket_out, Zmq.Poll.In) |]
+    in
+
+
+    while !run_status do
+
+        let polling =
+          Zmq.Poll.poll  ~timeout:1000  pollitem
+        in
+
+        match polling with
+          | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |]     -> ( action_out () ; action_in () )
+          | [| _                ; Some Zmq.Poll.In |]     ->   action_out () 
+          | [| Some Zmq.Poll.In ; _                |]     ->   action_in  () 
+          | _                    -> ()
+    done;
+
+    Zmq.Socket.close  socket_in;
+    Zmq.Socket.close  socket_out;
+  in
+
+
+
+  let ocaml_thread =
+    let addr_out =
+      Printf.sprintf "tcp:%s:%d"  address  port
+    in
+
+    let addr_in =
+      Printf.sprintf "tcp://*:%d"  localport
+    in
+
+    let f () = 
+      new_thread  REQ  addr_in  addr_out
+    in
+
+    (Thread.create f) ()
+  in
+  Printf.printf "Connect to:\ntcp://%s:%d\n%!" localhost localport;
+
+
+  let fortran_thread =
+    let addr_out =
+      Printf.sprintf "tcp:%s:%d" address (port+2)
+    in
+
+    let addr_in =
+      Printf.sprintf "tcp://*:%d" (localport+2)
+    in
+
+    let f () = 
+      new_thread  REQ  addr_in  addr_out
+    in
+    (Thread.create f) ()
+  in
+
+
+  let pub_thread =
+    let addr_in  =
+      Printf.sprintf "tcp:%s:%d" address (port+1)
+    in
+
+    let addr_out =
+      Printf.sprintf "tcp://*:%d" (localport+1)
+    in
+
+    let f () = 
+      new_thread  SUB  addr_in  addr_out
+    in
+    (Thread.create f) ()
+  in
+
+
+
+  let input_thread =
+    let f () = 
+      let addr_out =
+        match arg with
+        | EZFIO _ -> None
+        | ADDRESS _ -> Some (
+            Printf.sprintf "tcp:%s:%d"  address  (port+9) )
+      in
+
+      let addr_in =
+        Printf.sprintf "tcp://*:%d"  (localport+9)
+      in
+
+      let socket_in = 
+        create_socket  Zmq.Socket.rep  Zmq.Socket.bind     addr_in
+      in
+
+      let socket_out =
+        match addr_out with 
+        | Some addr_out -> Some (
+            create_socket  Zmq.Socket.req  Zmq.Socket.connect  addr_out)
+        | None -> None
+      in
+
+      let temp_file = 
+        Filename.temp_file "qp_tunnel" ".tar.gz"
+      in
+
+      let get_ezfio_filename () = 
+        match arg with
+        | EZFIO x -> x
+        | ADDRESS _ ->
+          begin
+            match socket_out with
+            | None -> assert false
+            | Some socket_out -> (
+                Zmq.Socket.send socket_out "get_ezfio_filename" ;
+                Zmq.Socket.recv socket_out
+)
+          end
+      in
+
+      let get_input () = 
+        match arg with
+        | EZFIO x -> 
+          begin
+            Printf.sprintf "tar --exclude=\"*.gz.*\" -zcf %s %s" temp_file x
+            |> Sys.command |> ignore;
+            let fd =
+              Unix.openfile  temp_file  [Unix.O_RDONLY]  0o640
+            in
+            let len =
+              Unix.lseek  fd  0  Unix.SEEK_END
+            in
+            ignore @@ Unix.lseek  fd  0  Unix.SEEK_SET ;
+            let bstr =
+                Unix.map_file  fd  Bigarray.char 
+                  Bigarray.c_layout false [| len |]
+                |> Bigarray.array1_of_genarray
+            in
+            let result = 
+              String.init len (fun i -> bstr.{i}) ;
+            in
+            Unix.close fd;
+            Sys.remove temp_file;
+            result
+          end
+        | ADDRESS _ ->
+          begin
+            match socket_out with
+            | None -> assert false
+            | Some socket_out -> (
+                Zmq.Socket.send socket_out "get_input" ;
+                Zmq.Socket.recv socket_out
+)
+          end
+      in
+
+      let () = 
+        match socket_out with
+        | None -> ()
+        | Some socket_out ->
+            Zmq.Socket.send  socket_out  "test";
+            Printf.printf "Communication [ %s ]\n%!" (Zmq.Socket.recv  socket_out);
+      in
+
+      (* Download input if asked *)
+      if Command_line.get_bool "get-input" then
+        begin
+          match arg with
+          | EZFIO _ -> ()
+          | ADDRESS _ ->
+            begin
+              Printf.printf "Getting input... %!";
+              let ezfio_filename = 
+                get_ezfio_filename ()
+              in
+              Printf.printf "%s%!" ezfio_filename;
+              let oc =
+                open_out temp_file
+              in
+              get_input ()
+              |> output_string oc;
+              close_out oc;
+              Printf.sprintf "tar -zxf %s" temp_file
+              |> Sys.command |> ignore ;
+              let oc =
+                Filename.concat  (Qpackage.ezfio_work ezfio_filename)  "qp_run_address"
+                |> open_out 
+              in
+              Printf.fprintf oc "tcp://%s:%d\n"  localhost  localport;
+              close_out oc;
+              Printf.printf " ...done\n%!"
+            end
+        end;
+
+      (* Main loop *)
+      let pollitem =
+        Zmq.Poll.mask_of [| (socket_in, Zmq.Poll.In) |]
+      in
+
+      let action () =
+        match Zmq.Socket.recv socket_in with
+        | "get_input" -> get_input ()
+                        |> Zmq.Socket.send socket_in 
+        | "get_ezfio_filename" -> get_ezfio_filename ()
+                        |> Zmq.Socket.send socket_in 
+        | "test" -> Zmq.Socket.send socket_in "OK"
+        | x -> Printf.sprintf "Message '%s' not understood" x
+               |> Zmq.Socket.send socket_in
+      in
+
+      Printf.printf "
+On remote hosts, create ssh tunnel using:
+ ssh -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d %s &
+Or from this host connect to clients using:
+ ssh -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d <host> &
+%!" 
+                (port  ) localhost (localport  )
+                (port+1) localhost (localport+1)
+                (port+2) localhost (localport+2)
+                (port+9) localhost (localport+9)
+                (Unix.gethostname ())
+                (port  ) (localport  )
+                (port+1) (localport+1)
+                (port+2) (localport+2)
+                (port+9) (localport+9);
+      Printf.printf "Ready\n%!";
+      while !run_status do
+
+          let polling =
+            Zmq.Poll.poll  ~timeout:1000  pollitem
+          in
+
+          match polling.(0) with
+            | Some Zmq.Poll.In     -> action ()
+            | None                 -> ()
+            | Some Zmq.Poll.In_out 
+            | Some Zmq.Poll.Out    -> ()
+
+      done;
+
+      let () = 
+        match socket_out with
+        | Some socket_out -> Zmq.Socket.close  socket_out
+        | None -> ()
+      in
+      Zmq.Socket.close  socket_in
+    in
+
+    (Thread.create f) ()
+  in
+
+  (* Termination *)
+  Thread.join input_thread;
+  Thread.join fortran_thread;
+  Thread.join pub_thread;
+  Thread.join ocaml_thread;
+  Zmq.Context.terminate zmq_context;
+  Printf.printf "qp_tunnel exited properly.\n"
+  
+
+
+

ocaml/qptypes_generator.ml

@@ -58,7 +58,7 @@ let input_data = "
 
 * Det_number_max : int
   assert (x > 0) ;
-  if (x > 50_00_000_000) then
+  if (x > 50_000_000_000) then
     warning \"More than 50 billion determinants\";
 
 * States_number : int

scripts/ezfio_interface/qp_edit_template

@@ -79,7 +79,7 @@ let get s =
          | Ao_basis ->
            f Ao_basis.(read, to_rst)
          | Determinants_by_hand ->
-           f Determinants_by_hand.(read_maybe, to_rst)
+           f Determinants_by_hand.(read ~full:false, to_rst)
 {section_to_rst}
       end
     with

src/ao_one_e_ints/README.rst

@@ -6,9 +6,8 @@ All the one-electron integrals in the |AO| basis are here.
 
 The most important providers for usual quantum-chemistry calculation are:
 
-* `ao_kinetic_integral` which are the kinetic operator integrals on the |AO| basis (see :file:`kin_ao_ints.irp.f`)
+* `ao_kinetic_integrals` which are the kinetic operator integrals on the |AO| basis 
-* `ao_nucl_elec_integral` which are the nuclear-elctron operator integrals on the |AO| basis (see :file:`pot_ao_ints.irp.f`)
+* `ao_integrals_n_e` which are the nuclear-elctron operator integrals on the |AO| basis
-* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis (see :file:`ao_mono_ints.irp.f`)
+* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis
 
 
-Note that you can find other interesting integrals related to the position operator in :file:`spread_dipole_ao.irp.f`.

src/bitmask/bitmask_cas_routines.irp.f

@@ -2,6 +2,29 @@ use bitmasks
 integer function number_of_holes(key_in)
  BEGIN_DOC
  ! Function that returns the number of holes in the inact space
+ !
+!   popcnt(
+!      xor(
+!        iand(
+!          reunion_of_core_inact_bitmask(1,1),
+!          xor(
+!            key_in(1,1),
+!            iand(
+!              key_in(1,1),
+!              cas_bitmask(1,1,1))
+!          )
+!        ),
+!        reunion_of_core_inact_bitmask(1,1)) )
+!
+! (key_in && cas_bitmask)
+! +---------------------+
+!    electrons in cas     xor key_in
+! +---------------------------------+
+!        electrons outside of cas     && reunion_of_core_inact_bitmask 
+! +------------------------------------------------------------------+ 
+!            electrons in the core/inact space     xor reunion_of_core_inact_bitmask
+! +---------------------------------------------------------------------------------+ 
+!              holes
  END_DOC
  implicit none
  integer(bit_kind), intent(in) :: key_in(N_int,2)
@@ -78,22 +101,6 @@ integer function number_of_holes(key_in)
    + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) )&
    + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,1), xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1)))), reunion_of_core_inact_bitmask(7,1)) )&
    + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,2), xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1)))), reunion_of_core_inact_bitmask(7,2)) )
- else if(N_int == 8)then
-    number_of_holes = number_of_holes  &
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(8,1), xor(key_in(8,1),iand(key_in(8,1),cas_bitmask(8,1,1)))), reunion_of_core_inact_bitmask(8,1)) )&
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(8,2), xor(key_in(8,2),iand(key_in(8,2),cas_bitmask(8,2,1)))), reunion_of_core_inact_bitmask(8,2)) )
  else
    do i = 1, N_int
     number_of_holes = number_of_holes             &
@@ -108,7 +115,7 @@ integer function number_of_holes(key_in)
             )                                     &
           )                                       &
         ), reunion_of_core_inact_bitmask(i,1)) )  &
-   + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) )
+   + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,2), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,2,1)))), reunion_of_core_inact_bitmask(i,2)) )
    enddo
  endif
 end

src/cipsi/run_pt2_slave.irp.f

@@ -141,8 +141,8 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
     b%cur=0
 
 !    ! Try to adjust n_tasks around nproc/2 seconds per job
-!    n_tasks = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0)))
+    n_tasks = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0)))
-    n_tasks = 1
+!    n_tasks = 1
   end do
 
   integer, external :: disconnect_from_taskserver

src/cipsi/run_selection_slave.irp.f

@@ -85,7 +85,7 @@ subroutine run_selection_slave(thread,iproc,energy)
       if(ctask > 0) then
         call sort_selection_buffer(buf)
 !        call merge_selection_buffers(buf,buf2)
-print *, task_id(1), pt2(1), buf%cur, ctask
+!print *, task_id(1), pt2(1), buf%cur, ctask
         call push_selection_results(zmq_socket_push, pt2, variance, norm, buf, task_id(1), ctask)
 !        buf%mini = buf2%mini
         pt2(:) = 0d0

src/determinants/EZFIO.cfg

@@ -66,6 +66,12 @@ interface: ezfio
 doc: Number of determinants in the current wave function
 type: integer
 
+[n_det_qp_edit]
+interface: ezfio
+doc: Number of determinants to print in qp_edit
+type: integer
+interface: ezfio
+
 [psi_coef]
 interface: ezfio
 doc: Coefficients of the wave function
@@ -78,6 +84,18 @@ doc: Determinants of the variational space
 type: integer*8
 size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det)
 
+[psi_coef_qp_edit]
+interface: ezfio
+doc: Coefficients of the wave function
+type: double precision 
+size: (determinants.n_det_qp_edit,determinants.n_states)
+
+[psi_det_qp_edit]
+interface: ezfio
+doc: Determinants of the variational space
+type: integer*8
+size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det_qp_edit)
+
 [expected_s2]
 interface: ezfio
 doc: Expected value of |S^2|

src/determinants/create_excitations.irp.f

@@ -12,6 +12,7 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok)
   integer(bit_kind), intent(inout) :: key_in(N_int,2)
   integer, intent(out)           :: i_ok
   integer                        :: k,j,i
+  integer(bit_kind)              :: mask
   use bitmasks
   ASSERT (i_hole > 0 )
   ASSERT (i_particle <= mo_num)
@@ -19,31 +20,66 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok)
   ! hole
   k = shiftr(i_hole-1,bit_kind_shift)+1
   j = i_hole-shiftl(k-1,bit_kind_shift)-1
+  mask = ibset(0_bit_kind,j)
 ! check whether position j is occupied
-  if (ibits(key_in(k,ispin),j,1).eq.1) then
+  if (iand(key_in(k,ispin),mask) /= 0_bit_kind) then
    key_in(k,ispin) = ibclr(key_in(k,ispin),j)
   else 
    i_ok= -1
+   return
   end if
 
   ! particle
   k = shiftr(i_particle-1,bit_kind_shift)+1
   j = i_particle-shiftl(k-1,bit_kind_shift)-1
-  key_in(k,ispin) = ibset(key_in(k,ispin),j)
+  mask = ibset(0_bit_kind,j)
+  if (iand(key_in(k,ispin),mask) == 0_bit_kind) then
+   key_in(k,ispin) = ibset(key_in(k,ispin),j)
+  else 
+   i_ok= -1
+   return
+  end if
 
-  integer                        :: n_elec_tmp
+!  integer                        :: n_elec_tmp
-  n_elec_tmp = 0
+!  n_elec_tmp = 0
-  do i = 1, N_int
+!  do i = 1, N_int
-    n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
+!    n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
-  enddo
+!  enddo
-  if(n_elec_tmp .ne. elec_num)then
+!  if(n_elec_tmp .ne. elec_num)then
-    !print*, n_elec_tmp,elec_num
+!    print*, n_elec_tmp,elec_num
-    !call debug_det(key_in,N_int)
+!    call debug_det(key_in,N_int)
-    i_ok = -1
+!    stop -1
-  endif
+!  endif
 end
 
 
+subroutine build_singly_excited_wavefunction(i_hole,i_particle,ispin,det_out,coef_out)
+  implicit none
+  BEGIN_DOC
+  ! Applies the single excitation operator : a^{dager}_(i_particle) a_(i_hole) of
+  ! spin = ispin to the current wave function (psi_det, psi_coef)
+  END_DOC
+  integer, intent(in)            :: i_hole,i_particle,ispin
+  integer(bit_kind), intent(out) :: det_out(N_int,2,N_det)
+  double precision, intent(out)  :: coef_out(N_det,N_states)
+
+  integer :: k
+  integer :: i_ok
+  double precision :: phase
+  do k=1,N_det
+    coef_out(k,:) = psi_coef(k,:)
+    det_out(:,:,k) = psi_det(:,:,k)
+    call do_single_excitation(det_out(1,1,k),i_hole,i_particle,ispin,i_ok)
+    if (i_ok == 1) then
+      call get_phase(psi_det(1,1,k), det_out(1,1,k),phase,N_int)
+      coef_out(k,:) = phase * coef_out(k,:)
+    else
+      coef_out(k,:) = 0.d0
+      det_out(:,:,k) = psi_det(:,:,k)
+    endif
+  enddo
+end
+
 logical function is_spin_flip_possible(key_in,i_flip,ispin)
   implicit none
   BEGIN_DOC

src/determinants/density_matrix.irp.f

@@ -269,7 +269,6 @@ subroutine save_natural_mos
    END_DOC
    call set_natural_mos
    call save_mos
-
 end
 
 

src/determinants/determinants.irp.f

@@ -44,6 +44,16 @@ BEGIN_PROVIDER [ integer, N_det ]
   ASSERT (N_det > 0)
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ integer, N_det_qp_edit  ]
+  implicit none
+  BEGIN_DOC
+! Number of determinants to print in qp_edit
+  END_DOC
+
+  N_det_qp_edit = min(N_det,10000)
+END_PROVIDER
+
 BEGIN_PROVIDER [integer, max_degree_exc]
   implicit none
   integer                        :: i,degree
@@ -476,7 +486,7 @@ subroutine save_wavefunction
   endif
   if (mpi_master) then
     call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
-    endif
+  endif
 end
 
 
@@ -504,12 +514,16 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
   integer*8, allocatable         :: psi_det_save(:,:,:)
   double precision, allocatable  :: psi_coef_save(:,:)
 
-  integer                        :: i,j,k
+  double precision               :: accu_norm
+  integer                        :: i,j,k, ndet_qp_edit
 
   if (mpi_master) then
+    ndet_qp_edit = min(ndet,N_det_qp_edit)
+
     call ezfio_set_determinants_N_int(N_int)
     call ezfio_set_determinants_bit_kind(bit_kind)
     call ezfio_set_determinants_N_det(ndet)
+    call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit)
     call ezfio_set_determinants_n_states(nstates)
     call ezfio_set_determinants_mo_label(mo_label)
 
@@ -522,10 +536,10 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
       enddo
     enddo
     call ezfio_set_determinants_psi_det(psi_det_save)
+    call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
     deallocate (psi_det_save)
 
     allocate (psi_coef_save(ndet,nstates))
-    double precision               :: accu_norm
     do k=1,nstates
       do i=1,ndet
         psi_coef_save(i,k) = psicoef(i,k)
@@ -535,6 +549,18 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
 
     call ezfio_set_determinants_psi_coef(psi_coef_save)
     deallocate (psi_coef_save)
+
+    allocate (psi_coef_save(ndet_qp_edit,nstates))
+    do k=1,nstates
+      do i=1,ndet_qp_edit
+        psi_coef_save(i,k) = psicoef(i,k)
+      enddo
+      call normalize(psi_coef_save(1,k),ndet_qp_edit)
+    enddo
+
+    call ezfio_set_determinants_psi_coef_qp_edit(psi_coef_save)
+    deallocate (psi_coef_save)
+
     call write_int(6,ndet,'Saved determinants')
   endif
 end
@@ -559,54 +585,80 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde
   integer                        :: N_int2
   equivalence (det_8, det_bk)
 
-  integer                        :: i,k
+  integer                        :: i,j,k, ndet_qp_edit
 
-  call ezfio_set_determinants_N_int(N_int)
+  if (mpi_master) then
-  call ezfio_set_determinants_bit_kind(bit_kind)
+    ndet_qp_edit = min(ndetsave,N_det_qp_edit)
-  call ezfio_set_determinants_N_det(ndetsave)
+    call ezfio_set_determinants_N_int(N_int)
-  call ezfio_set_determinants_n_states(nstates)
+    call ezfio_set_determinants_bit_kind(bit_kind)
-  call ezfio_set_determinants_mo_label(mo_label)
+    call ezfio_set_determinants_N_det(ndetsave)
+    call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit)
+    call ezfio_set_determinants_n_states(nstates)
+    call ezfio_set_determinants_mo_label(mo_label)
 
-  N_int2 = (N_int*bit_kind)/8
+    N_int2 = (N_int*bit_kind)/8
-  allocate (psi_det_save(N_int2,2,ndetsave))
+    allocate (psi_det_save(N_int2,2,ndetsave))
-  do i=1,ndetsave
-    do k=1,N_int
-      det_bk(k) = psidet(k,1,index_det_save(i))
-    enddo
-    do k=1,N_int2
-      psi_det_save(k,1,i) = det_8(k)
-    enddo
-    do k=1,N_int
-      det_bk(k) = psidet(k,2,index_det_save(i))
-    enddo
-    do k=1,N_int2
-      psi_det_save(k,2,i) = det_8(k)
-    enddo
-  enddo
-  call ezfio_set_determinants_psi_det(psi_det_save)
-  deallocate (psi_det_save)
-
-  allocate (psi_coef_save(ndetsave,nstates))
-  double precision               :: accu_norm(nstates)
-  accu_norm = 0.d0
-  do k=1,nstates
     do i=1,ndetsave
-      accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k)
+      do k=1,N_int
-      psi_coef_save(i,k) = psicoef(index_det_save(i),k)
+        det_bk(k) = psidet(k,1,index_det_save(i))
+      enddo
+      do k=1,N_int2
+        psi_det_save(k,1,i) = det_8(k)
+      enddo
+      do k=1,N_int
+        det_bk(k) = psidet(k,2,index_det_save(i))
+      enddo
+      do k=1,N_int2
+        psi_det_save(k,2,i) = det_8(k)
+      enddo
     enddo
-  enddo
+    call ezfio_set_determinants_psi_det(psi_det_save)
-  do k = 1, nstates
+    call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
-    accu_norm(k) = 1.d0/dsqrt(accu_norm(k))
+    deallocate (psi_det_save)
-  enddo
-  do k=1,nstates
-    do i=1,ndetsave
-      psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k)
-    enddo
-  enddo
 
-  call ezfio_set_determinants_psi_coef(psi_coef_save)
+    allocate (psi_coef_save(ndetsave,nstates))
-  call write_int(6,ndet,'Saved determinants')
+    double precision               :: accu_norm(nstates)
-  deallocate (psi_coef_save)
+    accu_norm = 0.d0
+    do k=1,nstates
+      do i=1,ndetsave
+        accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k)
+        psi_coef_save(i,k) = psicoef(index_det_save(i),k)
+      enddo
+    enddo
+    do k = 1, nstates
+      accu_norm(k) = 1.d0/dsqrt(accu_norm(k))
+    enddo
+    do k=1,nstates
+      do i=1,ndetsave
+        psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k)
+      enddo
+    enddo
+
+    call ezfio_set_determinants_psi_coef(psi_coef_save)
+    deallocate (psi_coef_save)
+
+    allocate (psi_coef_save(ndet_qp_edit,nstates))
+    accu_norm = 0.d0
+    do k=1,nstates
+      do i=1,ndet_qp_edit
+        accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k)
+        psi_coef_save(i,k) = psicoef(index_det_save(i),k)
+      enddo
+    enddo
+    do k = 1, nstates
+      accu_norm(k) = 1.d0/dsqrt(accu_norm(k))
+    enddo
+    do k=1,nstates
+      do i=1,ndet_qp_edit
+        psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k)
+      enddo
+    enddo
+
+    call ezfio_set_determinants_psi_coef(psi_coef_save)
+    deallocate (psi_coef_save)
+
+    call write_int(6,ndet,'Saved determinants')
+  endif
 end
 
 

src/determinants/occ_pattern.irp.f

@@ -53,7 +53,17 @@ subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint)
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Generate all possible determinants for a give occ_pattern
+  ! Generate all possible determinants for a given occ_pattern
+  ! 
+  ! Input :
+  !    o   : occupation pattern : (doubly occupied, singly occupied)
+  !    sze : Number of produced determinants, computed by `occ_pattern_to_dets_size`
+  !    n_alpha : Number of $\alpha$ electrons
+  !    Nint    : N_int
+  !
+  ! Output:
+  !    d : determinants 
+  !
   END_DOC
   integer          ,intent(in)   :: Nint
   integer          ,intent(in)   :: n_alpha        ! Number of alpha electrons

src/dft_utils_in_r/kin_dens.irp.f

@@ -0,0 +1,15 @@
+BEGIN_PROVIDER [double precision, kinetic_density_generalized,  (n_points_final_grid)]
+ implicit none
+ integer :: i,j,m,i_point
+ kinetic_density_generalized = 0.d0
+ do i_point = 1, n_points_final_grid
+  do i = 1, mo_num
+   do j = 1, mo_num 
+    do m = 1, 3
+     kinetic_density_generalized(i_point) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i_point) * mos_grad_in_r_array_tranp(m,i,i_point) * one_e_dm_mo_for_dft(j,i,1) 
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 

src/dft_utils_in_r/mo_in_r.irp.f

@@ -22,6 +22,7 @@
 
 
  BEGIN_PROVIDER[double precision, mos_grad_in_r_array,(mo_num,n_points_final_grid,3)]
+&BEGIN_PROVIDER[double precision, mos_grad_in_r_array_tranp,(3,mo_num,n_points_final_grid)]
  implicit none
  BEGIN_DOC
  ! mos_grad_in_r_array(i,j,k)          = value of the kth component of the gradient of ith mo on the jth grid point
@@ -35,6 +36,35 @@
  do m=1,3
   call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_grad_in_r_array(1,1,m),mo_num)
  enddo
+ integer  :: i,j
+ do i = 1, n_points_final_grid
+  do j = 1, mo_num
+   do m = 1, 3
+     mos_grad_in_r_array_tranp(m,j,i) = mos_grad_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, alpha_dens_kin_in_r, (n_points_final_grid)]
+&BEGIN_PROVIDER [double precision, beta_dens_kin_in_r, (n_points_final_grid)]
+ implicit none
+ integer  :: i,m,j
+ alpha_dens_kin_in_r = 0.d0
+ beta_dens_kin_in_r = 0.d0
+ do i = 1, n_points_final_grid
+  do j = 1, elec_alpha_num
+   do m = 1, 3
+    alpha_dens_kin_in_r(i) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i)**2.d0 
+   enddo
+  enddo
+  do j = 1, elec_beta_num
+   do m = 1, 3
+    beta_dens_kin_in_r(i)  += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i)**2.d0
+   enddo
+  enddo
+ enddo
+ 
  END_PROVIDER 
 
  BEGIN_PROVIDER[double precision, mos_lapl_in_r_array,(mo_num,n_points_final_grid,3)]

src/dft_utils_one_e/ec_lyp.irp.f

@@ -35,11 +35,8 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
  include 'constants.include.F'
 
 ! Input variables
-
  double precision, intent(in) :: rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2
-
 ! Local variables
-
  double precision :: a,b,c,d,c_f,omega,delta
  double precision :: rho_13,rho_inv_13,rho_83,rho_113,rho_inv_113,denom
  double precision :: thr,huge_num,rho_inv
@@ -47,8 +44,6 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
  double precision :: tmp1,tmp2,tmp3,tmp4
  double precision :: big1,big2,big3
 
-! Output variables
-
 
 ! Constants of the LYP correlation functional
 
@@ -57,15 +52,25 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
  c = 0.2533d0
  d = 0.349d0
 
- thr = 1d-10
+ ec_lyp_88 = 0.d0
+
+ thr = 1d-15
  huge_num = 1.d0/thr
+ if(dabs(rho_a).lt.thr)then
+  return 
+ endif
+
+ if(dabs(rho_b).lt.thr)then
+  return 
+ endif
 
  if(rho.lt.0.d0)then
   print*,'pb !! rho.lt.0.d0'
   stop
  endif
- rho_13  = rho**(1d0/3d0)
+
- rho_113 = rho**(11d0/3d0)
+ rho_13  = rho**(1.d0/3.d0)
+ rho_113 = rho**(11.d0/3.d0)
 
  if(dabs(rho_13) < thr) then
   rho_inv_13 = huge_num
@@ -76,13 +81,13 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
  if (dabs(rho_113) < thr) then
   rho_inv_113 = huge_num
  else
-  rho_inv_113 = 1d0/rho_113
+  rho_inv_113 = 1.d0/rho_113
  endif
 
  if (dabs(rho) < thr) then
   rho_inv = huge_num
  else
-  rho_inv = 1d0/rho
+  rho_inv = 1.d0/rho
  endif
 
 ! Useful quantities to predefine
@@ -90,21 +95,21 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
  denom = 1d0/(1d0 + d*rho_inv_13)
  omega = rho_inv_113*exp(-c*rho_inv_13)*denom
  delta = c*rho_inv_13 + d*rho_inv_13*denom
- c_f   = 0.3d0*(3d0*pi*pi)**(2d0/3d0)
+ c_f   = 0.3d0*(3.d0*pi*pi)**(2.d0/3.d0)
 
  rho_2   = rho  *rho
  rho_a_2 = rho_a*rho_a
  rho_b_2 = rho_b*rho_b
 
- cst_2_113 = 2d0**(11d0/3d0)
+ cst_2_113 = 2.d0**(11.d0/3.d0)
- cst_8_3   = 8d0/3d0
+ cst_8_3   = 8.d0/3.d0
 
  ! first term in the equation (2) of Preuss CPL, 1989
 
- big1 = 4d0*denom*rho_a*rho_b*rho_inv
+ big1 = 4.d0*denom*rho_a*rho_b*rho_inv
 
  tmp1 = cst_2_113*c_f*(rho_a**cst_8_3 + rho_b**cst_8_3)
- tmp2 = (47d0/18d0 - 7d0/18d0*delta)*grad_rho_2
+ tmp2 = (47.d0/18.d0 - 7.d0/18.d0*delta)*grad_rho_2
  tmp3 = - (5d0/2d0 - 1.d0/18d0*delta)*(grad_rho_a_2 + grad_rho_b_2)
  tmp4 = - (delta - 11d0)/9d0*(rho_a*rho_inv*grad_rho_a_2 + rho_b*rho_inv*grad_rho_b_2)
  big2 = rho_a*rho_b*(tmp1 + tmp2 + tmp3 + tmp4)
@@ -114,7 +119,6 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr
  tmp3 = grad_rho_a_2*(2d0/3d0*rho_2 - rho_b_2)
  big3 = tmp1 + tmp2 + tmp3
 
-
  ec_lyp_88 = -a*big1 -a*b*omega*big2 -a*b*omega*big3
 
 end

src/dft_utils_one_e/ec_scan.irp.f

@@ -0,0 +1,100 @@
+double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2)
+ include 'constants.include.F'
+ implicit none
+ double precision, intent(in) :: rho_a,rho_b,tau,grad_rho_2
+ double precision :: cst_13,cst_23,cst_43,cst_53,rho_inv,cst_18,cst_3pi2
+ double precision :: thr,nup,ndo,xi,s,spin_d,drho,drho2,rho,inv_1alph,e_c_lsda1,h0
+ double precision :: rs,t_w,t_unif,ds_xi,alpha,fc_alpha,step_f,cst_1alph,beta_inf
+ double precision :: c_1c,c_2c,d_c,e_c_ldsa1,h1,phi,t,beta_rs,gama,a,w_1,g_at2,phi_3,e_c_1
+ double precision :: b_1c,b_2c,b_3c,dx_xi,gc_xi,e_c_lsda0,w_0,g_inf,cx_xi,x_inf,f0,e_c_0
+ thr = 1.d-12
+ nup = max(rho_a,thr)
+ ndo = max(rho_b,thr)
+ rho = nup + ndo 
+ ec_scan = 0.d0
+ if((rho).lt.thr)return
+ ! constants ...
+ rho_inv  = 1.d0/rho
+ cst_13   = 1.d0/3.d0
+ cst_23   = 2.d0 * cst_13
+ cst_43   = 4.d0 * cst_13
+ cst_53   = 5.d0 * cst_13
+ cst_18   = 1.d0/8.d0
+ cst_3pi2 = 3.d0 * pi*pi
+ drho2    = max(grad_rho_2,thr)
+ drho     = dsqrt(drho2)
+ if((nup-ndo).gt.0.d0)then
+  spin_d  = max(nup-ndo,thr) 
+ else
+  spin_d  = min(nup-ndo,-thr) 
+ endif
+ c_1c     = 0.64d0
+ c_2c     = 1.5d0
+ d_c      = 0.7d0
+ b_1c     = 0.0285764d0
+ b_2c     = 0.0889d0
+ b_3c     = 0.125541d0
+ gama     = 0.031091d0
+ ! correlation energy lsda1
+ call ec_only_lda_sr(0.d0,nup,ndo,e_c_lsda1)
+          
+ xi       = spin_d/rho
+ rs       = (cst_43 * pi * rho)**(-cst_13)
+ s        = drho/( 2.d0 * cst_3pi2**(cst_13) * rho**cst_43  )
+ t_w      = drho2 * cst_18 * rho_inv
+ ds_xi    = 0.5d0 * ( (1.d0+xi)**cst_53 + (1.d0 - xi)**cst_53)
+ t_unif   = 0.3d0 * (cst_3pi2)**cst_23 * rho**cst_53*ds_xi
+ t_unif   = max(t_unif,thr)
+ alpha    = (tau - t_w)/t_unif
+ cst_1alph= 1.d0 - alpha
+ if(cst_1alph.gt.0.d0)then
+  cst_1alph= max(cst_1alph,thr)
+ else
+  cst_1alph= min(cst_1alph,-thr)
+ endif
+ inv_1alph= 1.d0/cst_1alph
+ phi      = 0.5d0 * ( (1.d0+xi)**cst_23 + (1.d0 - xi)**cst_23)
+ phi_3    = phi*phi*phi
+ t        = (cst_3pi2/16.d0)**cst_13 * s / (phi * rs**0.5d0)
+ w_1      = dexp(-e_c_lsda1/(gama * phi_3)) - 1.d0
+ a        = beta_rs(rs) /(gama * w_1)
+ g_at2    = 1.d0/(1.d0 + 4.d0 * a*t*t)**0.25d0
+ h1       = gama * phi_3 * dlog(1.d0 + w_1 * (1.d0 - g_at2))
+ ! interpolation function 
+ fc_alpha = dexp(-c_1c * alpha * inv_1alph) * step_f(cst_1alph) - d_c * dexp(c_2c * inv_1alph) * step_f(-cst_1alph) 
+ ! first part of the correlation energy 
+ e_c_1    = e_c_lsda1 + h1
+ 
+ dx_xi    =  0.5d0 * ( (1.d0+xi)**cst_43 + (1.d0 - xi)**cst_43)
+ gc_xi    = (1.d0 - 2.3631d0 * (dx_xi - 1.d0) ) * (1.d0 - xi**12.d0)
+ e_c_lsda0= - b_1c / (1.d0 + b_2c * rs**0.5d0 + b_3c * rs)
+ w_0      = dexp(-e_c_lsda0/b_1c) - 1.d0
+ beta_inf = 0.066725d0 * 0.1d0 / 0.1778d0
+ cx_xi    = -3.d0/(4.d0*pi) * (9.d0 * pi/4.d0)**cst_13 * dx_xi
+
+ x_inf   = 0.128026d0
+ f0       = -0.9d0
+ g_inf    = 1.d0/(1.d0 + 4.d0 * x_inf * s*s)**0.25d0
+ 
+ h0       = b_1c * dlog(1.d0 + w_0 * (1.d0 - g_inf)) 
+ e_c_0    = (e_c_lsda0 + h0) * gc_xi
+
+ ec_scan = e_c_1 + fc_alpha * (e_c_0 - e_c_1)
+end
+
+double precision function step_f(x)
+ implicit none
+ double precision, intent(in) :: x
+ if(x.lt.0.d0)then
+  step_f = 0.d0
+ else
+  step_f = 1.d0 
+ endif
+end             
+
+double precision function beta_rs(rs)
+ implicit none
+ double precision, intent(in) ::rs
+ beta_rs = 0.066725d0 * (1.d0 + 0.1d0 * rs)/(1.d0 + 0.1778d0 * rs)
+
+end

src/fci/40.fci.bats

@@ -36,24 +36,24 @@ function run_stoch() {
   [[ -n $TRAVIS ]] && skip
   qp set_file f2.ezfio
   qp set_frozen_core
-  run_stoch -199.30496 1.e-4
+  run_stoch -199.30486   1.e-4
 }
 
 @test "NH3" { # 10.6657s
   qp set_file nh3.ezfio
   qp set_mo_class --core="[1-4]" --act="[5-72]"
-  run -56.244753429144986  1.e-5
+  run -56.244753429144986  1.e-4
 }
 
 @test "DHNO" { # 11.4721s
   qp set_file dhno.ezfio
   qp set_mo_class --core="[1-7]" --act="[8-64]" 
-  run -130.459020029816 1.e-5
+  run -130.459020029816 1.e-4
 }
 
 @test "HCO" { # 12.2868s
   qp set_file hco.ezfio
-  run -113.297494345682 2.e-05 
+  run -113.297494345682 1.e-4
 }
 
 @test "H2O2" { # 12.9214s
@@ -65,82 +65,82 @@ function run_stoch() {
 @test "HBO" { # 13.3144s
   [[ -n $TRAVIS ]] && skip
   qp set_file hbo.ezfio
-  run -100.214185815312 1.e-5
+  run -100.212829869715 1.e-4
 }
 
 @test "H2O" { # 11.3727s
   [[ -n $TRAVIS ]] && skip
   qp set_file h2o.ezfio
-  run -76.2359268957699 2.e-5
+  run -76.2359268957699 1.e-4
 }
 
 @test "ClO" { # 13.3755s
   [[ -n $TRAVIS ]] && skip
   qp set_file clo.ezfio
-  run -534.546005867797 5.e-5
+  run -534.545881614967 1.e-4
 }
 
 @test "SO" { # 13.4952s
   [[ -n $TRAVIS ]] && skip
   qp set_file so.ezfio
-  run -26.0124797722154 1.e-5
+  run -26.0158153138924 1.e-4
 }
 
 @test "H2S" { # 13.6745s
   [[ -n $TRAVIS ]] && skip
   qp set_file h2s.ezfio
-  run -398.859480581924 1.e-5
+  run -398.859168655255 1.e-4
 }
 
 @test "OH" { # 13.865s 
   [[ -n $TRAVIS ]] && skip
   qp set_file oh.ezfio
-  run -75.6119887538831 1.e-05 
+  run -75.6120779012574 1.e-4 
 }
 
 @test "SiH2_3B1" { # 13.938ss
   [[ -n $TRAVIS ]] && skip
   qp set_file sih2_3b1.ezfio
-  run -290.017539006762 1.e-5
+  run -290.017539006762 1.e-4
 }
 
 @test "H3COH" { # 14.7299s
   [[ -n $TRAVIS ]] && skip
   qp set_file h3coh.ezfio
-  run -115.205054063687 1.e-5
+  run -115.205941463667 1.e-4
 }
 
 @test "SiH3" { # 15.99s
   [[ -n $TRAVIS ]] && skip
   qp set_file sih3.ezfio
-  run -5.57269434557089 2.e-05 
+  run -5.57241217753818 1.e-4 
 }
 
 @test "CH4" { # 16.1612s
   [[ -n $TRAVIS ]] && skip
   qp set_file ch4.ezfio
   qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]"
-  run -40.2409059687324 2.e-5
+  run -40.2409678239136 1.e-4
 }
 
 @test "ClF" { # 16.8864s
   [[ -n $TRAVIS ]] && skip
   qp set_file clf.ezfio
-  run -559.170406471496 1.e-5
+  run -559.170272077166 1.e-4
 }
 
 @test "SO2" { # 17.5645s
   [[ -n $TRAVIS ]] && skip
   qp set_file so2.ezfio
   qp set_mo_class --core="[1-8]" --act="[9-87]" 
-  run -41.5746738713298 5.e-5
+  run -41.5746738713298 1.e-4
 }
 
 @test "C2H2" { # 17.6827s
   [[ -n $TRAVIS ]] && skip
   qp set_file c2h2.ezfio
   qp set_mo_class --act="[1-30]" --del="[31-36]"
-  run -12.3670840202635 2.e-5
+  run -12.3656179738175 1.e-4
 }
 
 @test "N2" { # 18.0198s
@@ -154,14 +154,14 @@ function run_stoch() {
   [[ -n $TRAVIS ]] && skip
   qp set_file n2h4.ezfio
   qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]"
-  run -111.367234092521 2.e-5
+  run -111.367332681559 1.e-4
 }
 
 @test "CO2" { # 21.1748s
   [[ -n $TRAVIS ]] && skip
   qp set_file co2.ezfio
   qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]"
-  run -187.969676381867 1.e-5
+  run -187.968599504402 1.e-4
 }
 
 
@@ -169,13 +169,13 @@ function run_stoch() {
   [[ -n $TRAVIS ]] && skip
   qp set_file cu_nh3_4_2plus.ezfio
   qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
-  run -1862.98610987882  1.e-05
+  run -1862.98614665139  1.e-04
 }
 
 @test "HCN" { # 20.3273s
   [[ -n $TRAVIS ]] && skip
   qp set_file hcn.ezfio
   qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]"
-  run -93.0799328685679 2.e-5
+  run -93.0728641601823 1.e-4
 }
 

src/nuclei/EZFIO.cfg

@@ -21,7 +21,7 @@ type: double precision
 size: (nuclei.nucl_num,3)
 interface: ezfio
 
-[disk_access_nuclear_repulsion]
+[io_nuclear_repulsion]
 doc: Read/Write Nuclear Repulsion from/to disk [ Write | Read | None ]
 type: Disk_access
 interface: ezfio,provider,ocaml

src/nuclei/nuclei.irp.f

@@ -143,7 +143,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
    END_DOC
 
    PROVIDE mpi_master nucl_coord nucl_charge nucl_num
-   if (disk_access_nuclear_repulsion.EQ.'Read') then
+   if (io_nuclear_repulsion == 'Read') then
      logical                        :: has
 
      if (mpi_master) then
@@ -194,7 +194,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
    call write_time(6)
    call write_double(6,nuclear_repulsion,'Nuclear repulsion energy')
 
-   if (disk_access_nuclear_repulsion.EQ.'Write') then
+   if (io_nuclear_repulsion == 'Write') then
      if (mpi_master) then
        call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
      endif

src/tools/print_ci_vectors.irp.f

@@ -24,6 +24,7 @@ subroutine routine
  implicit none
  integer :: i,k
  integer :: degree
+ call print_energy_components
  do i = 1, N_det
   print *, 'Determinant ', i
   call debug_det(psi_det(1,1,i),N_int)