Merge http://github.com/scemama/quantum_package

.travis.yml

@@ -4,19 +4,22 @@
 #  - sudo apt-get install gfortran liblapack-dev gcc
 #  - sudo apt-get install graphviz
 
+os: linux
+
 dist: trusty
 
 sudo: false
 
+compiler: gfortran
+
 addons:
   apt:
     packages:
     - gfortran
     - gcc
-    - liblapack-dev
+#    - liblapack-dev
+#    - libblas-dev
     - graphviz
-#    - zlib1g-dev
-#    - libgmp3-dev
 
 cache:
   directories:
@@ -29,6 +32,7 @@ python:
 script: 
   - ./configure  ./config/travis.cfg
   - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD_ZMQ mrcepa0 All_singles
+  - source ./quantum_package.rc ; cd ~ ; install_lapack.sh ; cd $QP_ROOT
   - 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

@@ -1,3 +1,10 @@
+## IMPORTANT
+
+If you have problems upgrading to the current version, consider re-installing everything from scratch including the OCaml compiler.
+To do this, you will have to remove the `quantum_package` directory **and** the `$HOME/.opam` directory as well.
+
+
+
 ![QP](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/qp.png)     
 [![Build Status](https://travis-ci.org/LCPQ/quantum_package.svg?branch=master)](https://travis-ci.org/LCPQ/quantum_package)
 [![Gitter](https://badges.gitter.im/Join%20Chat.svg)](https://gitter.im/LCPQ/quantum_package?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)

bin/qp_gaspi_run

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+QP_ROOT=$( cd $(dirname "${BASH_SOURCE}")/..  ; pwd -P )
+source $HOME/.bashrc
+source $QP_ROOT/quantum_package.rc
+exec $QP_ROOT/ocaml/qp_run $@

config/gfortran.cfg

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

config/gfortran_debug.cfg

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

config/travis.cfg

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

configure

@@ -43,6 +43,7 @@ except KeyError:
 
 QP_ROOT_BIN = join(QP_ROOT, "bin")
 QP_ROOT_LIB = join(QP_ROOT, "lib")
+QP_ROOT_LIB64 = join(QP_ROOT, "lib64")
 QP_ROOT_INSTALL = join(QP_ROOT, "install")
 
 os.environ["PATH"] = os.environ["PATH"] + ":" + QP_ROOT_BIN
@@ -65,6 +66,7 @@ d_dependency = {
     "python": [],
     "ninja": ["g++", "python"],
     "make": [],
+    "gpi2": ["g++", "make"],
     "p_graphviz": ["python"],
     "bats": []
 }
@@ -140,6 +142,11 @@ f77zmq = Info(
     description=' F77-ZeroMQ',
     default_path=join(QP_ROOT_LIB, "libf77zmq.a") )
 
+gpi2 = Info( 
+    url='https://github.com/cc-hpc-itwm/GPI-2/archive/v1.3.0.tar.gz',
+    description=' GPI-2',
+    default_path=join(QP_ROOT_LIB64, "libGPI2.a") )
+
 p_graphviz = Info(
     url='https://github.com/xflr6/graphviz/archive/master.tar.gz',
     description=' Python library for graphviz',
@@ -154,7 +161,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", "gpi2"]: 
     exec ("d_info['{0}']={0}".format(m))
 
 
@@ -481,8 +488,8 @@ def create_ninja_and_rc(l_installed):
         'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")),
         'export PYTHONPATH="${QP_EZFIO}/Python":"${QP_PYTHON}":"${PYTHONPATH}"',
         'export PATH="${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml:"${PATH}"',
-        'export LD_LIBRARY_PATH="${QP_ROOT}"/lib:"${LD_LIBRARY_PATH}"',
+        'export LD_LIBRARY_PATH="${QP_ROOT}"/lib:"${QP_ROOT}"/lib64:"${LD_LIBRARY_PATH}"',
-        'export LIBRARY_PATH="${QP_ROOT}"/lib:"${LIBRARY_PATH}"',
+        'export LIBRARY_PATH="${QP_ROOT}"/lib:"${QP_ROOT}"/lib64:"${LIBRARY_PATH}"',
         'export C_INCLUDE_PATH="${C_INCLUDE_PATH}":"${QP_ROOT}"/include',
         '',
         'source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh',

install/scripts/install_gpi2.sh

@@ -0,0 +1,15 @@
+#!/bin/bash -x
+
+TARGET=gpi2
+#GPI_OPTIONS=--with-infiniband 
+GPI_OPTIONS=--with-ethernet
+
+function _install()
+{
+  cd _build/gpi2
+  ./install.sh -p $QP_ROOT $GPI_OPTIONS
+  cp src/GASPI.f90 $QP_ROOT/plugins/GPI2/
+  return 0
+}
+
+source scripts/build.sh

install/scripts/install_lapack.sh

@@ -0,0 +1,10 @@
+#!/bin/bash -x
+
+git clone https://github.com/Reference-LAPACK/lapack-release.git
+cd lapack-release
+cp make.inc.example make.inc
+make -j 8 
+mv librefblas.a liblapack.a libtmglib.a $QP_ROOT/lib
+
+
+

install/scripts/install_ocaml.sh

@@ -5,8 +5,7 @@ QP_ROOT=$PWD
 cd -
 
 # Normal installation
-PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ"
+PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ ppx_sexp_conv ppx_deriving"
-#ppx_sexp_conv
 
 # Needed for ZeroMQ
 export C_INCLUDE_PATH="${QP_ROOT}"/include:"${C_INCLUDE_PATH}"
@@ -64,7 +63,7 @@ fi
 cd Downloads || exit 1
 chmod +x  ocaml.sh || exit 1
 
-echo N | ./ocaml.sh ${QP_ROOT}/bin/ 4.02.1 || exit 1
+echo N | ./ocaml.sh ${QP_ROOT}/bin/ 4.04.2 || exit 1
 
 ${QP_ROOT}/bin/opam config setup -a -q || exit 1
 
@@ -76,7 +75,6 @@ source ${HOME}/.opam/opam-init/init.sh > /dev/null 2> /dev/null || true
 
 NCPUs=$(cat /proc/cpuinfo  | grep -i  MHz | wc -l)
 ${QP_ROOT}/bin/opam install -j ${NCPUs} stdint.0.4.2 -y -q || exit 1
-${QP_ROOT}/bin/opam install -j ${NCPUs} ZMQ -y -q || exit 1
 ${QP_ROOT}/bin/opam install -j ${NCPUs} ${PACKAGES} -y -q || exit 1
 
 rm -f ../_build/ocaml.log

ocaml/Address.ml

@@ -1,4 +1,3 @@
-open Core.Std
 
 module Tcp : sig
   type t 
@@ -8,7 +7,7 @@ module Tcp : sig
 end = struct
   type t = string
   let of_string x =
-    if not (String.is_prefix ~prefix:"tcp://" x) then
+    if not (String_ext.is_prefix ~prefix:"tcp://" x) then
       invalid_arg "Address Invalid"
     ;
     x
@@ -26,7 +25,7 @@ module Ipc : sig
 end = struct
   type t = string
   let of_string x =
-    assert (String.is_prefix ~prefix:"ipc://" x);
+    assert (String_ext.is_prefix ~prefix:"ipc://" x);
     x
   let create name =
     Printf.sprintf "ipc://%s" name
@@ -41,7 +40,7 @@ module Inproc : sig
 end = struct
   type t = string
   let of_string x =
-    assert (String.is_prefix ~prefix:"inproc://" x);
+    assert (String_ext.is_prefix ~prefix:"inproc://" x);
     x
   let create name =
     Printf.sprintf "inproc://%s" name

ocaml/Atom.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 exception AtomError of string
 
@@ -6,7 +6,7 @@ type t =
 { element : Element.t ;
   charge  : Charge.t ;
   coord   : Point3d.t ;
-} with sexp
+} [@@deriving sexp]
 
 (** Read xyz coordinates of the atom *)
 let of_string ~units s =

ocaml/Basis.ml

@@ -1,7 +1,7 @@
-open Core.Std
+open Sexplib.Std
 open Qptypes
 
-type t = (Gto.t * Nucl_number.t) list with sexp
+type t = (Gto.t * Nucl_number.t) list [@@deriving sexp]
 
 (** Read all the basis functions of an element *)
 let read in_channel at_number =
@@ -16,7 +16,7 @@ let read in_channel at_number =
 
 (** Find an element in the basis set file *)
 let find in_channel element =
-  In_channel.seek in_channel 0L;
+  seek_in in_channel 0;
   let element_read = ref Element.X in
   while !element_read <> element
   do
@@ -56,13 +56,13 @@ let to_string_general ~fmt ~atom_sep ?ele_array b =
     do_work ((Gto.to_string ~fmt g)::accu) n tail
   in
   do_work [new_nucleus 1] 1 b
-  |> String.concat ~sep:"\n"
+  |> String.concat "\n"
 
 let to_string_gamess ?ele_array =
     to_string_general ?ele_array ~fmt:Gto.Gamess ~atom_sep:""
 
 let to_string_gaussian ?ele_array b =
-  String.concat ~sep:"\n"
+  String.concat "\n"
   [ to_string_general ?ele_array ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ]
 
 let to_string ?(fmt=Gto.Gamess) =

ocaml/Bit.ml

@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core;;
 
 (*
 Type for bits 
@@ -11,7 +11,7 @@ Zero | One
 type t      = 
 | One 
 | Zero 
-with sexp
+[@@deriving sexp]
 
 let to_string = function
   | Zero -> "0"

ocaml/Bit.mli

@@ -1,4 +1,4 @@
-type t = One | Zero with sexp
+type t = One | Zero [@@deriving sexp]
 
 (** String conversions for printing *)
 val to_string : t -> string

ocaml/Bitlist.ml

@@ -1,5 +1,5 @@
 open Qptypes
-open Core.Std
+open Core
 
 (*
 Type for bits strings

ocaml/Charge.ml

@@ -1,6 +1,6 @@
-open Core.Std
+open Core
 
-type t = float with sexp
+type t = float [@@deriving sexp]
 
 let of_float x = x
 let of_int   i = Float.of_int i

ocaml/Charge.mli

@@ -1,4 +1,4 @@
-type t = float with sexp
+type t = float [@@deriving sexp]
 
 (** Float conversion functions *)
 val to_float : t -> float

ocaml/Determinant.ml

@@ -1,7 +1,7 @@
-open Core.Std;;
+open Qptypes
-open Qptypes;;
+open Sexplib.Std
 
-type t = int64 array with sexp
+type t = int64 array [@@deriving sexp]
 
 let to_int64_array (x:t) = (x:int64 array)
 
@@ -9,8 +9,8 @@ let to_int64_array (x:t) = (x:int64 array)
 let to_alpha_beta x = 
   let x = to_int64_array x in
   let n_int = (Array.length x)/2 in
-  ( Array.init n_int ~f:(fun i -> x.(i)) ,
+  ( Array.init n_int (fun i -> x.(i)) ,
-    Array.init n_int ~f:(fun i -> x.(i+n_int)) )
+    Array.init n_int (fun i -> x.(i+n_int)) )
 
 
 let to_bitlist_couple x =
@@ -28,12 +28,14 @@ let bitlist_to_string ~mo_tot_num x =
   let len = 
     MO_number.to_int mo_tot_num
   in
-  List.map x ~f:(function
+  let s = 
-      | Bit.Zero -> "-"
+    List.map (function
-      | Bit.One  -> "+" 
+        | Bit.Zero -> "-"
-  )
+        | Bit.One  -> "+" 
-  |> String.concat
+    ) x
-  |> String.sub ~pos:0 ~len 
+    |> String.concat ""
+  in
+  String.sub s 0 len 
 
 
 
@@ -77,6 +79,6 @@ let to_string ~mo_tot_num x =
   let (xa,xb) = to_bitlist_couple x in
   [ "  " ; bitlist_to_string ~mo_tot_num xa ; "\n" ;
     "  " ; bitlist_to_string ~mo_tot_num xb ]
-  |> String.concat 
+  |> String.concat ""
 
 

ocaml/Determinant.mli

@@ -5,7 +5,7 @@
   * where each int64 is a list of 64 MOs. When the bit is set
   * to 1, the MO is occupied.
   *)
-type t = int64 array with sexp
+type t = int64 array [@@deriving sexp]
 
 (** Transform to an int64 array *)
 val to_int64_array : t -> int64 array

ocaml/Element.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 exception ElementError of string
@@ -11,7 +11,7 @@ type t =
 |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr
 |Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe
                            |Pt          
-with sexp
+[@@deriving sexp]
 
 let of_string x = 
   match (String.capitalize (String.lowercase x)) with

ocaml/Element.mli

@@ -8,7 +8,7 @@ type t =
 |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr
 |Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe
                            |Pt
-with sexp
+[@@deriving sexp]
 
 (** String conversion functions *)
 val of_string : string -> t

ocaml/Excitation.ml

@@ -1,14 +1,14 @@
-open Core.Std;;
+open Core
-open Qptypes;;
+open Qptypes
 
 module Hole = struct
-    type t = MO_class.t with sexp
+    type t = MO_class.t [@@deriving sexp]
     let of_mo_class x = x
     let to_mo_class x = x
 end
 
 module Particle = struct
-    type t = MO_class.t with sexp
+    type t = MO_class.t [@@deriving sexp]
     let of_mo_class x = x
     let to_mo_class x = x
 end
@@ -16,7 +16,7 @@ end
 type t =
 | Single of Hole.t*Particle.t
 | Double of Hole.t*Particle.t*Hole.t*Particle.t
-with sexp;;
+[@@deriving sexp]
 
 let create_single ~hole ~particle =
   MO_class.(
@@ -29,7 +29,7 @@ let create_single ~hole ~particle =
   | (          _, Inactive _ ) -> failwith "Particles can not be in virtual MOs"
   | (h, p) -> Single ( (Hole.of_mo_class h), (Particle.of_mo_class p) )
   )
-;;
+
 
 let double_of_singles s1 s2 =
   let (h1,p1) = match s1 with
@@ -40,14 +40,14 @@ let double_of_singles s1 s2 =
   | _ -> assert false
   in
   Double (h1,p1,h2,p2)
-;;
+
 
 let create_double ~hole1 ~particle1 ~hole2 ~particle2 =
   let s1 = create_single ~hole:hole1 ~particle:particle1
   and s2 =  create_single ~hole:hole2 ~particle:particle2
   in
   double_of_singles s1 s2
-;;
+
 
 let to_string = function
   | Single (h,p) ->
@@ -68,5 +68,5 @@ let to_string = function
         (MO_class.to_string (Particle.to_mo_class p2));
         "]"]
       |> String.concat ~sep:" "
-;;
+
 

ocaml/Excitation.mli

@@ -18,7 +18,7 @@ module Particle :
 type t =
   | Single of Hole.t * Particle.t
   | Double of Hole.t * Particle.t * Hole.t * Particle.t
-with sexp
+[@@deriving sexp]
 
 val create_single : hole:MO_class.t -> particle:MO_class.t -> t
 

ocaml/GaussianPrimitive.ml

@@ -1,10 +1,10 @@
 open Qptypes
-open Core.Std
+open Core
 
 type t =
 { sym  : Symmetry.t ;
   expo : AO_expo.t ;
-} with sexp
+} [@@deriving sexp]
 
 let to_string p = 
   let { sym = s ; expo = e } = p in

ocaml/Generic_input_of_rst.ml

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

ocaml/Gto.ml

@@ -1,5 +1,5 @@
-open Core.Std
 open Qptypes
+open Sexplib.Std
 
 exception GTO_Read_Failure of string
 exception End_Of_Basis
@@ -11,11 +11,11 @@ type fmt =
 type t =
 { sym  : Symmetry.t ;
   lc   : ((GaussianPrimitive.t * AO_coef.t) list)
-} with sexp
+} [@@deriving sexp]
 
 
 let of_prim_coef_list pc =
-  let (p,c) = List.hd_exn pc in
+  let (p,c) = List.hd pc in
   let sym = p.GaussianPrimitive.sym in
   let rec check = function
   | [] -> `OK
@@ -37,12 +37,12 @@ let of_prim_coef_list pc =
 let read_one in_channel =
   (* Fetch number of lines to read on first line *)
   let buffer = input_line in_channel in
-  if ( (String.strip buffer) = "" ) then
+  if ( (String_ext.strip buffer) = "" ) then
      raise End_Of_Basis;
   let sym_str = String.sub buffer 0 2 in
   let   n_str = String.sub buffer 2 ((String.length buffer)-2) in
-  let sym = Symmetry.of_string (String.strip sym_str) in
+  let sym = Symmetry.of_string (String_ext.strip sym_str) in
-  let n = Int.of_string (String.strip n_str) in
+  let n = int_of_string (String_ext.strip n_str) in
   (* Read all the primitives *)
   let rec read_lines result = function
   | 0 -> result
@@ -50,18 +50,19 @@ let read_one in_channel =
     begin
       let line_buffer = input_line in_channel in
       let buffer = line_buffer 
-      |> String.split ~on:' ' 
+      |> String_ext.split ~on:' ' 
-      |> List.filter ~f:(fun x -> x <> "")
+      |> List.filter (fun x -> x <> "")
       in
       match buffer with
       | [ j ; expo ; coef ] ->
         begin
-          let coef = String.tr ~target:'D' ~replacement:'e' coef
+          let coef =
+            Str.global_replace (Str.regexp "D") "e" coef
           in
           let p =
             GaussianPrimitive.of_sym_expo sym 
-              (AO_expo.of_float (Float.of_string expo) )
+              (AO_expo.of_float (float_of_string expo) )
-          and c = AO_coef.of_float (Float.of_string coef) in
+          and c = AO_coef.of_float (float_of_string coef) in
           read_lines ( (p,c)::result) (i-1)
         end
       | _ -> raise (GTO_Read_Failure line_buffer)
@@ -89,7 +90,7 @@ let to_string_gamess { sym = sym ; lc = lc } =
     do_work (result::accu) (i+1) tail
   in
   (do_work [result] 1 lc)
-  |> String.concat ~sep:"\n"
+  |> String.concat "\n"
 
 
 (** Write the GTO in Gaussian format *)
@@ -109,7 +110,7 @@ let to_string_gaussian { sym = sym ; lc = lc } =
     do_work (result::accu) (i+1) tail
   in
   (do_work [result] 1 lc)
-  |> String.concat ~sep:"\n"
+  |> String.concat "\n"
 
 
 (** Transform the gto to a string *)

ocaml/Gto.mli

@@ -7,7 +7,7 @@ type fmt =
 type t =
   { sym : Symmetry.t ;
     lc  : (GaussianPrimitive.t * Qptypes.AO_coef.t) list;
-  } with sexp
+  } [@@deriving sexp]
 
 (** Create from a list of GaussianPrimitive.t * Qptypes.AO_coef.t *)
 val of_prim_coef_list :

ocaml/Input.ml

@@ -1,6 +1,6 @@
 open Qputils;;
 open Qptypes;;
-open Core.Std;;
+open Core;;
 
 include Input_ao_basis;;
 include Input_bitmasks;;

ocaml/Input_ao_basis.ml

@@ -1,6 +1,6 @@
 open Qptypes;;
 open Qputils;;
-open Core.Std;;
+open Core;;
 
 module Ao_basis : sig
   type t = 
@@ -13,7 +13,7 @@ module Ao_basis : sig
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
-    } with sexp
+    } [@@deriving sexp]
   ;;
   val read : unit -> t option
   val to_string : t -> string
@@ -32,7 +32,7 @@ end = struct
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
-    } with sexp
+    } [@@deriving sexp]
   ;;
 
   let get_default = Qpackage.get_ezfio_default "ao_basis";;

ocaml/Input_bi_integrals.ml

@@ -1,6 +1,6 @@
 open Qptypes;;
 open Qputils;;
-open Core.Std;;
+open Core;;
 
 module Bielec_integrals : sig
   type t = 
@@ -11,7 +11,7 @@ module Bielec_integrals : sig
       threshold_ao       : Threshold.t;
       threshold_mo       : Threshold.t;
       direct             : bool;
-    } with sexp
+    } [@@deriving sexp]
   ;;
   val read  : unit -> t option
   val write : t -> unit
@@ -27,7 +27,7 @@ end = struct
       threshold_ao       : Threshold.t;
       threshold_mo       : Threshold.t;
       direct             : bool;
-    } with sexp
+    } [@@deriving sexp]
   ;;
 
   let get_default = Qpackage.get_ezfio_default "bielec_integrals";;

ocaml/Input_bitmasks.ml

@@ -1,6 +1,6 @@
 open Qptypes;;
 open Qputils;;
-open Core.Std;;
+open Core;;
 
 module Bitmasks : sig
   type t = 
@@ -10,7 +10,7 @@ module Bitmasks : sig
       generators         : int64 array;
       n_mask_cas         : Bitmask_number.t;
       cas                : int64 array;
-    } with sexp
+    } [@@deriving sexp]
   ;;
   val read : unit -> t option 
   val to_string : t -> string
@@ -22,7 +22,7 @@ end = struct
       generators         : int64 array;
       n_mask_cas         : Bitmask_number.t;
       cas                : int64 array;
-    } with sexp
+    } [@@deriving sexp]
   ;;
 
   let get_default = Qpackage.get_ezfio_default "bitmasks";;

ocaml/Input_determinants_by_hand.ml

@@ -1,6 +1,6 @@
 open Qptypes;;
 open Qputils;;
-open Core.Std;;
+open Core;;
 
 module Determinants_by_hand : sig
   type t = 
@@ -11,7 +11,7 @@ module Determinants_by_hand : sig
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
-    } with sexp
+    } [@@deriving sexp]
   val read  : unit -> t
   val read_maybe  : unit -> t option
   val write : t -> unit
@@ -30,7 +30,7 @@ end = struct
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
-    } with sexp
+    } [@@deriving sexp]
   ;;
 
   let get_default = Qpackage.get_ezfio_default "determinants";;

ocaml/Input_electrons.ml

@@ -1,12 +1,12 @@
 open Qptypes;;
 open Qputils;;
-open Core.Std;;
+open Core;;
 
 module Electrons : sig
   type t = 
     { elec_alpha_num     : Elec_alpha_number.t;
       elec_beta_num      : Elec_beta_number.t;
-    } with sexp
+    } [@@deriving sexp]
   ;;
   val read  : unit -> t option
   val write : t -> unit
@@ -18,7 +18,7 @@ end = struct
   type t = 
     { elec_alpha_num     : Elec_alpha_number.t;
       elec_beta_num      : Elec_beta_number.t;
-    } with sexp
+    } [@@deriving sexp]
   ;;
 
   let get_default = Qpackage.get_ezfio_default "electrons";;

ocaml/Input_mo_basis.ml

@@ -1,23 +1,29 @@
 open Qptypes
 open Qputils
-open Core.Std
+open Core
 
-type t_mo = 
-    { mo_tot_num      : MO_number.t ;
-      mo_label        : MO_label.t;
-      mo_class        : MO_class.t array;
-      mo_occ          : MO_occ.t array;
-      mo_coef         : (MO_coef.t array) array;
-      ao_md5          : MD5.t;
-    } with sexp
 
 module Mo_basis : sig
-  type t = t_mo
+  type t = 
+      { mo_tot_num      : MO_number.t ;
+        mo_label        : MO_label.t;
+        mo_class        : MO_class.t array;
+        mo_occ          : MO_occ.t array;
+        mo_coef         : (MO_coef.t array) array;
+        ao_md5          : MD5.t;
+      } [@@deriving sexp]
   val read : unit -> t option
   val to_string : t -> string
   val to_rst : t -> Rst_string.t
 end = struct
-  type t = t_mo
+  type t = 
+      { mo_tot_num      : MO_number.t ;
+        mo_label        : MO_label.t;
+        mo_class        : MO_class.t array;
+        mo_occ          : MO_occ.t array;
+        mo_coef         : (MO_coef.t array) array;
+        ao_md5          : MD5.t;
+      } [@@deriving sexp]
   let get_default = Qpackage.get_ezfio_default "mo_basis"
 
   let read_mo_label () = 

ocaml/Input_nuclei_by_hand.ml

@@ -1,6 +1,6 @@
 open Qptypes;;
 open Qputils;;
-open Core.Std;;
+open Core;;
 
 module Nuclei_by_hand : sig
   type t = 
@@ -8,7 +8,7 @@ module Nuclei_by_hand : sig
       nucl_label      : Element.t array;
       nucl_charge     : Charge.t array;
       nucl_coord      : Point3d.t array;
-    } with sexp
+    } [@@deriving sexp]
   ;;
   val read  : unit -> t option
   val write : t -> unit 
@@ -22,7 +22,7 @@ end = struct
       nucl_label      : Element.t array;
       nucl_charge     : Charge.t array;
       nucl_coord      : Point3d.t array;
-    } with sexp
+    } [@@deriving sexp]
   ;;
 
   let get_default = Qpackage.get_ezfio_default "nuclei";;

ocaml/Io_ext.ml

@@ -0,0 +1,24 @@
+let input_lines filename =
+  let in_channel =
+    open_in filename
+  in
+  let rec aux accu = 
+    try
+      let newline = 
+        input_line in_channel
+      in
+      aux (newline::accu)
+    with End_of_file -> accu
+  in
+  let result = 
+    List.rev (aux [])
+  in
+  close_in in_channel;
+  result
+
+
+
+let read_all filename = 
+  input_lines filename
+  |> String.concat "\n"
+

ocaml/Long_basis.ml

@@ -1,7 +1,7 @@
-open Core.Std;;
+open Qptypes
-open Qptypes;;
+open Sexplib.Std
 
-type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list with sexp
+type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp]
 
 let of_basis b =
   let rec do_work accu = function
@@ -10,14 +10,14 @@ let of_basis b =
         begin
           let new_accu = 
             Symmetry.Xyz.of_symmetry g.Gto.sym 
-            |> List.rev_map ~f:(fun x-> (x,g,n)) 
+            |> List.rev_map (fun x-> (x,g,n)) 
           in
           do_work (new_accu@accu) tail
         end
   in
   do_work [] b
   |> List.rev
-;;
+
 
 let to_basis b =
   let rec do_work accu = function
@@ -25,7 +25,7 @@ let to_basis b =
   | (s,g,n)::tail -> 
     let first_sym = 
       Symmetry.Xyz.of_symmetry g.Gto.sym
-      |> List.hd_exn
+      |> List.hd
     in
     let new_accu = 
       if ( s = first_sym ) then
@@ -36,19 +36,19 @@ let to_basis b =
     do_work new_accu tail
   in
   do_work [] b 
-;;
+
 
 let to_string b =
-  let middle = List.map ~f:(fun (x,y,z) -> 
+  let middle = List.map (fun (x,y,z) -> 
-     "( "^((Int.to_string (Nucl_number.to_int z)))^", "^
+     "( "^((string_of_int (Nucl_number.to_int z)))^", "^
      (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y)
      ^" )"
   ) b
-  |> String.concat ~sep:",\n"
+  |> String.concat ",\n"
   in "("^middle^")"
-;;
 
-include To_md5;;
+
+include To_md5
 let to_md5 = to_md5 sexp_of_t
-;;
+
 

ocaml/Long_basis.mli

@@ -5,7 +5,7 @@ open Qptypes;;
   * all the D orbitals are converted to xx, xy, xz, yy, yx
   * etc
 *)
-type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list with sexp
+type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp]
 
 (** Transform a basis to a long basis *)
 val of_basis :

ocaml/MO_class.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes 
 
 type t = 
@@ -7,7 +7,7 @@ type t =
   | Active   of MO_number.t list
   | Virtual  of MO_number.t list
   | Deleted  of MO_number.t list
-with sexp
+[@@deriving sexp]
 
 
 let to_string  x = 

ocaml/MO_class.mli

@@ -4,7 +4,7 @@ type t =
   | Active   of Qptypes.MO_number.t list
   | Virtual  of Qptypes.MO_number.t list
   | Deleted  of Qptypes.MO_number.t list
-with sexp
+[@@deriving sexp]
 
 
 (** Create different excitation classes *)

ocaml/MO_label.ml

@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core;;
 
 type t =
 | Guess
@@ -7,7 +7,7 @@ type t =
 | Localized
 | Orthonormalized
 | None
-with sexp
+[@@deriving sexp]
 ;;
 
 let to_string = function

ocaml/MO_label.mli

@@ -5,7 +5,7 @@ type t =
   | Localized
   | Orthonormalized
   | None
-with sexp
+[@@deriving sexp]
 
 (** String representation *)
 val to_string : t -> string

ocaml/Makefile

@@ -11,8 +11,8 @@ endif
 
 LIBS=
 PKGS=
-OCAMLCFLAGS="-g -warn-error A"
+OCAMLCFLAGS="-g"
-OCAMLBUILD=ocamlbuild -j 0 -syntax camlp4o -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) 
+OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) 
 MLLFILES=$(wildcard *.mll) 
 MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml 
 MLIFILES=$(wildcard *.mli) git

ocaml/Message.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 (** New job : Request to create a new multi-tasked job *)

ocaml/Molecule.ml

@@ -1,14 +1,14 @@
-open Core.Std ;;
+open Qptypes 
-open Qptypes ;;
+open Sexplib.Std
 
-exception MultiplicityError of string;;
+exception MultiplicityError of string
-exception XYZError ;;
+exception XYZError 
 
 type t = {
   nuclei     : Atom.t list ;
   elec_alpha : Elec_alpha_number.t ;
   elec_beta  : Elec_beta_number.t ;
-} with sexp
+} [@@deriving sexp]
 
 let get_charge { nuclei  ; elec_alpha ; elec_beta } =
   let result =
@@ -19,7 +19,7 @@ let get_charge { nuclei  ; elec_alpha ; elec_beta } =
   | a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest
   | [] -> 0.
   in
-  Charge.of_float (nucl_charge nuclei  -. (Float.of_int result))
+  Charge.of_float (nucl_charge nuclei  -. (float_of_int result))
 
 
 let get_multiplicity m = 
@@ -59,9 +59,10 @@ let name m =
   | a::rest ->
       begin
         let e = a.Atom.element in
-        match (List.Assoc.find accu e) with
+        try
-        | None   -> build_list (List.Assoc.add accu e 1) rest
+          let i = List.assoc e accu in
-        | Some i -> build_list (List.Assoc.add accu e (i+1)) rest
+          build_list ( (e,i+1)::(List.remove_assoc e accu) ) rest
+        with Not_found -> build_list ( (e,1)::accu ) rest
       end
   | [] -> accu
   in
@@ -83,7 +84,7 @@ let name m =
   let result =
      build_list [] nuclei |> build_name [c ; ", " ; mult]
   in
-  String.concat (result)
+  String.concat "" result
 
 
 let to_string_general ~f m = 
@@ -95,8 +96,8 @@ let to_string_general ~f m =
   let title =
      name m
   in
-  [ Int.to_string n ; title ] @  (List.map ~f nuclei)
+  [ string_of_int n ; title ] @  (List.map f nuclei)
-  |> String.concat ~sep:"\n"
+  |> String.concat "\n"
 
 let to_string =
   to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x)
@@ -109,9 +110,9 @@ let of_xyz_string
     ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
     ?(units=Units.Angstrom)
     s =
-  let l = String.split s ~on:'\n'
+  let l = String_ext.split s ~on:'\n'
-       |> List.filter ~f:(fun x -> x <> "")
+       |> List.filter (fun x -> x <> "")
-       |> List.map ~f:(fun x -> Atom.of_string units x)
+       |> List.map (fun x -> Atom.of_string units x)
   in
   let ne = ( get_charge { 
         nuclei=l ;
@@ -145,25 +146,28 @@ let of_xyz_file
     ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
     ?(units=Units.Angstrom)
     filename =
-  let (x,buffer) = In_channel.read_all filename 
+  let lines = 
-  |> String.lsplit2_exn ~on:'\n'
+    match Io_ext.input_lines filename with
+    | natoms :: title :: rest -> 
+        begin
+          try
+            if (int_of_string @@ String_ext.strip natoms) <= 0 then
+              raise XYZError
+          with
+          | _ -> raise XYZError
+        end;
+        String.concat "\n" rest
+    | _ -> raise XYZError
   in
-  let result =
+  of_xyz_string ~charge:charge ~multiplicity:multiplicity 
-    try
+    ~units:units lines
-      (int_of_string @@ String.strip x) > 0
-    with
-    | Failure "int_of_string"  -> false
-  in
-  if not result then raise XYZError;
-  let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in
-  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
+  Io_ext.read_all filename
   |> Zmatrix.of_string
   |> Zmatrix.to_xyz_string
   |> of_xyz_string ~charge ~multiplicity ~units
@@ -182,14 +186,14 @@ let of_file
 let distance_matrix molecule =
   let coord =
     molecule.nuclei
-    |> List.map ~f:(fun x -> x.Atom.coord)
+    |> List.map (fun x -> x.Atom.coord)
     |> Array.of_list
   in
   let n = 
     Array.length coord
   in
   let result =
-    Array.make_matrix ~dimx:n  ~dimy:n 0.
+    Array.make_matrix n n 0.
   in
   for i = 0 to (n-1)
   do
@@ -203,6 +207,7 @@ let distance_matrix molecule =
   
   
 
+open Core ;;
 include To_md5
 let to_md5 = to_md5 sexp_of_t
 

ocaml/Molecule.mli

@@ -4,7 +4,7 @@ type t = {
   nuclei     : Atom.t list;
   elec_alpha : Qptypes.Elec_alpha_number.t;
   elec_beta  : Qptypes.Elec_beta_number.t;
-} with sexp
+} [@@deriving sexp]
 
 (** Returns the charge of the molecule *)
 val get_charge       : t -> Charge.t

ocaml/Multiplicity.ml

@@ -1,7 +1,7 @@
-open Core.Std;;
+open Core;;
 open Qptypes ;;
 
-type t = Strictly_positive_int.t with sexp
+type t = Strictly_positive_int.t [@@deriving sexp]
 
 let of_int = Strictly_positive_int.of_int ;;
 let to_int = Strictly_positive_int.to_int ;;

ocaml/Multiplicity.mli

@@ -1,4 +1,4 @@
-type t = Qptypes.Strictly_positive_int.t with sexp
+type t = Qptypes.Strictly_positive_int.t [@@deriving sexp]
 
 (** Conversion from int *)
 val of_int : int -> t

ocaml/Point3d.ml

@@ -1,11 +1,11 @@
-open Core.Std;;
+open Core;;
 open Qptypes;;
 
 type t = {
   x : float ;
   y : float ;
   z : float ;
-} with sexp
+} [@@deriving sexp]
 
 let of_tuple ~units (x,y,z) =
   let f = match units with 

ocaml/Point3d.mli

@@ -2,7 +2,7 @@ type t =
 { x : float;
   y : float;
   z : float;
-} with sexp
+} [@@deriving sexp]
 
 (** Create from a tuple of floats *)
 val of_tuple  : units:Units.units -> float*float*float -> t

ocaml/Primitive.mli

@@ -1,7 +1,7 @@
 type t =
 { sym : Symmetry.t;
   expo : Qptypes.AO_expo.t;
-} with sexp
+} [@@deriving sexp]
 
 (** Conversion to string for printing *)
 val to_string : t -> string

ocaml/Progress_bar.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 type t =
 { 
@@ -53,13 +53,13 @@ let display_tty bar =
     in
     let stop_time =
       let x = 
-        Time.Span.to_float running_time
+        Time.Span.to_sec running_time
       in
       if (percent > 0.) then
         x *. 100. /. percent -. x
-        |> Time.Span.of_float
+        |> Time.Span.of_sec
       else
-        Time.Span.of_float 0.
+        Time.Span.of_sec 0.
     in
     Printf.eprintf "%s : [%s] %4.1f%%  | %10s, ~%10s left\r%!"
       bar.title
@@ -67,7 +67,7 @@ let display_tty bar =
       percent
       (Time.Span.to_string running_time)
       (stop_time |> Time.Span.to_string );
-    { bar with dirty = false ; next = Time.add now (Time.Span.of_float 0.1) }
+    { bar with dirty = false ; next = Time.add now (Time.Span.of_sec 0.1) }
   
 
 let display_file bar =
@@ -80,19 +80,19 @@ let display_file bar =
     in
     let stop_time =
       let x = 
-        Time.Span.to_float running_time
+        Time.Span.to_sec running_time
       in
       if (percent > 0.) then
         x *. 100. /. percent -. x
-        |> Time.Span.of_float
+        |> Time.Span.of_sec
       else
-        Time.Span.of_float 0.
+        Time.Span.of_sec 0.
     in
     Printf.eprintf "%5.2f %%  in  %20s, ~%20s left\n%!"
       percent
       (Time.Span.to_string running_time)
       (Time.Span.to_string stop_time);
-    { bar with dirty = false ; next = Time.add (Time.now ()) (Time.Span.of_float 2.) }
+    { bar with dirty = false ; next = Time.add (Time.now ()) (Time.Span.of_sec 2.) }
   
 
 

ocaml/Pseudo.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 
@@ -7,7 +7,7 @@ module GaussianPrimitive_local : sig
   type t = { 
     expo    : AO_expo.t ;
     r_power : R_power.t ;
-  } with sexp
+  } [@@deriving sexp]
 
   val of_expo_r_power : AO_expo.t -> R_power.t -> t
   val to_string : t -> string
@@ -17,7 +17,7 @@ end = struct
   type t = { 
     expo    : AO_expo.t ;
     r_power : R_power.t ;
-  } with sexp
+  } [@@deriving sexp]
 
   let of_expo_r_power dz n = 
      { expo = dz ; r_power = n }
@@ -35,7 +35,7 @@ module GaussianPrimitive_non_local : sig
     expo    : AO_expo.t ;
     r_power : R_power.t ;
     proj    : Positive_int.t
-  } with sexp
+  } [@@deriving sexp]
 
   val of_proj_expo_r_power : Positive_int.t -> AO_expo.t -> R_power.t -> t
   val to_string : t -> string
@@ -46,7 +46,7 @@ end = struct
     expo    : AO_expo.t ;
     r_power : R_power.t ;
     proj    : Positive_int.t
-  } with sexp
+  } [@@deriving sexp]
 
   let of_proj_expo_r_power p dz n = 
      { expo = dz ; r_power = n ; proj = p }
@@ -66,7 +66,7 @@ type t = {
   n_elec    : Positive_int.t ;
   local     : (GaussianPrimitive_local.t * AO_coef.t ) list ;
   non_local : (GaussianPrimitive_non_local.t * AO_coef.t ) list 
-} with sexp
+} [@@deriving sexp]
 
 let empty e =
   { element   = e;

ocaml/Qpackage.ml

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

ocaml/Qputils.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Sexplib
 
 (*
 let rec transpose = function
@@ -14,12 +14,12 @@ let rec transpose = function
 
 let input_to_sexp s =
     let result = 
-      String.split_lines s 
+      String_ext.split ~on:'\n' s 
-      |> List.filter ~f:(fun x->
+      |> List.filter (fun x-> (String_ext.strip x) <> "")
-         (String.strip x) <> "")
+      |> List.map (fun x-> "("^ 
-      |> List.map ~f:(fun x->
+        (Str.global_replace (Str.regexp "=") " " x)
-         "("^(String.tr '=' ' ' x)^")")
+        ^")")
-      |> String.concat 
+      |> String.concat ""
     in
     print_endline ("("^result^")");
     "("^result^")"
@@ -29,10 +29,10 @@ let rmdir dirname =
   let rec remove_one dir =
     Sys.chdir dir;
     Sys.readdir "."
-    |> Array.iter ~f:(fun x ->
+    |> Array.iter (fun x ->
-      match (Sys.is_directory x, Sys.is_file x) with
+      match (Sys.is_directory x, Sys.file_exists x) with
-      | (`Yes, _) -> remove_one x
+      | (true, _) -> remove_one x
-      | (_, `Yes) -> Sys.remove x
+      | (_, true) -> Sys.remove x
       | _ -> failwith ("Unable to remove file "^x^".")
     );
     Sys.chdir "..";

ocaml/README_qp_edit.rst

@@ -19,7 +19,7 @@ of the block.
          r_y    : Y_type.t
         ...
          last_r : bool
-      } with sexp
+      } [@@deriving sexp]
     ;;
     val read  : unit -> t
     val write : t -> unit
@@ -31,7 +31,7 @@ of the block.
          r_y    : Y_type.t
         ...
          last_r : bool
-      } with sexp
+      } [@@deriving sexp]
     ;;
 
     let get_default = Qpackage.get_ezfio_default "new_keyword";;

ocaml/Range.ml

@@ -1,4 +1,4 @@
-open Core.Std;;
+open Sexplib.Std
 
 (* A range is a string of the type:
  * 
@@ -12,14 +12,14 @@ open Core.Std;;
 *)
 
 
-type t = int list with sexp
+type t = int list [@@deriving sexp] 
 
 let expand_range r =
-  match String.lsplit2 ~on:'-' r with
+  match String_ext.lsplit2 ~on:'-' r with
   | Some (s, f) ->
       begin
-        let start = Int.of_string s
+        let start = int_of_string s
-        and finish =  Int.of_string f
+        and finish =  int_of_string f
         in
         assert (start <= finish) ;
         let rec do_work = function
@@ -31,9 +31,9 @@ let expand_range r =
       begin
         match r with
           | "" -> []
-          | _  -> [Int.of_string r] 
+          | _  -> [int_of_string r] 
       end
-;;
+
 
 let of_string s =
   match s.[0] with
@@ -43,36 +43,37 @@ let of_string s =
     assert (s.[0] = '[') ;
     assert (s.[(String.length s)-1] = ']') ;
     let s = String.sub s 1 ((String.length s) - 2) in
-    let l = String.split ~on:',' s in
+    let l = String_ext.split ~on:',' s in
-    let l = List.map ~f:expand_range l in
+    let l = List.map expand_range l in
-  List.concat l |> List.dedup ~compare:Int.compare |> List.sort ~cmp:Int.compare
+  List.concat l
-;;
+  |> List.sort_uniq compare 
+
 
 let to_string l =
   let rec do_work buf symbol = function
     | [] -> buf
     | a::([] as t) -> 
-          do_work (buf^symbol^(Int.to_string a)) "" t
+          do_work (buf^symbol^(string_of_int a)) "" t
     | a::(b::q as t) ->
         if (b-a = 1) then
           do_work buf "-" t
         else
-          do_work (buf^symbol^(Int.to_string a)^","^(Int.to_string b)) "" t
+          do_work (buf^symbol^(string_of_int a)^","^(string_of_int b)) "" t
   in
   let result = 
     match l with
     | [] ->
         "[]"
     | h::t  ->
-        do_work ("["^(Int.to_string h)) "" l in
+        do_work ("["^(string_of_int h)) "" l in
   (String.sub result 0 ((String.length result)))^"]"
-;;
+
 
 let test_module () =
   let s = "[72-107,36-53,126-131]" in
   let l = of_string s in
   print_string s ; print_newline () ;
-  List.iter ~f:(fun x -> Printf.printf "%d, " x) l ; print_newline () ;
+  List.iter (fun x -> Printf.printf "%d, " x) l ; print_newline () ;
-  to_string l |> print_string ;  print_newline () ;
+  to_string l |> print_string ;  print_newline ();
-;;
+
 

ocaml/Range.mli

@@ -1,4 +1,4 @@
-type t = int list with sexp
+type t = int list [@@deriving sexp] 
 
 (** A range is a sorted list of ints in an interval.
     It is created using a string :

ocaml/String_ext.ml

@@ -0,0 +1,142 @@
+include String
+
+(** Split a string on a given character *)
+let split ?(on=' ') str = 
+  split_on_char on str
+(*
+    let rec do_work ?(accu=[]) ?(left="") = function
+    | "" -> List.rev (left::accu)
+    | s  ->
+        let new_s =
+            (length s) - 1
+            |> sub s 1 
+        in
+        if (s.[0] = on) then
+            let new_accu = 
+                left :: accu
+            in
+            do_work ~accu:new_accu new_s
+        else
+            let new_left =
+                concat "" [ left ; make 1 s.[0] ]
+            in
+            do_work ~accu ~left:new_left new_s
+    in
+    do_work str
+*)
+  
+
+(** Strip blanks on the left of a string *)
+let ltrim s =
+    let rec do_work s l =
+        match s.[0] with
+        | '\n'
+        | ' ' -> do_work (sub s 1 (l-1)) (l-1)
+        | _ -> s
+    in
+    let l =
+        length s
+    in
+    if (l > 0) then
+        do_work s l
+    else
+        s
+
+(** Strip blanks on the right of a string *)
+let rtrim s =
+    let rec do_work s l =
+        let newl =
+            l-1
+        in
+        match s.[newl] with
+        | '\n'
+        | ' ' -> do_work (sub s 0 (newl)) (newl)
+        | _ -> s
+    in
+    let l = 
+        length s
+    in
+    if (l > 0) then
+        do_work s l
+    else
+        s
+
+
+(** Strip blanks on the right and left of a string *)
+let strip = String.trim
+
+
+(** Split a string in two pieces when a character is found the 1st time from the left *)
+let lsplit2_exn ?(on=' ') s =
+    let length = 
+        String.length s
+    in
+    let rec do_work i = 
+        if (i = length) then
+           begin
+              raise Not_found
+           end
+        else if (s.[i] = on) then
+           ( String.sub s 0 i, 
+             String.sub s (i+1) (length-i-1) )
+        else
+           do_work (i+1)
+    in
+    do_work 0
+
+
+(** Split a string in two pieces when a character is found the 1st time from the right *)
+let rsplit2_exn ?(on=' ') s =
+    let length = 
+        String.length s
+    in
+    let rec do_work i = 
+        if (i = -1) then
+           begin
+              raise Not_found
+           end
+        else if (s.[i] = on) then
+           ( String.sub s 0 i, 
+             String.sub s (i+1) (length-i-1) )
+        else
+           do_work (i-1)
+    in
+    do_work length
+
+
+let lsplit2 ?(on=' ') s =
+  try
+    Some (lsplit2_exn ~on s) 
+  with
+  | Not_found -> None
+
+
+let rsplit2 ?(on=' ') s =
+  try
+    Some (rsplit2_exn ~on s) 
+  with
+  | Not_found -> None
+
+
+let to_list s =
+  Array.init (String.length s) (fun i -> s.[i])
+  |> Array.to_list
+
+
+let fold ~init ~f s =
+    to_list s
+    |> List.fold_left f init 
+
+  
+let is_prefix ~prefix s =
+  let len =
+    String.length prefix
+  in
+  if len > String.length s then
+    false
+  else
+    prefix = String.sub s 0 len 
+
+
+let of_char c =
+  String.make 1 c

ocaml/Symmetry.ml

@@ -1,7 +1,7 @@
 open Qptypes
-open Core.Std
+open Sexplib.Std
 
-type t = S|P|D|F|G|H|I|J|K|L with sexp
+type t = S|P|D|F|G|H|I|J|K|L [@@deriving sexp]
 
 let to_string = function
   | S -> "S"
@@ -77,7 +77,7 @@ type st = t
 module Xyz = struct
   type t = { x: Positive_int.t ;
              y: Positive_int.t ;
-             z: Positive_int.t } with sexp
+             z: Positive_int.t } [@@deriving sexp]
   type state_type = Null | X | Y | Z
 
   (** Builds an XYZ triplet from a string.
@@ -86,7 +86,7 @@ module Xyz = struct
     let flush state accu number =
       let n = 
         if (number = "") then 1
-        else (Int.of_string number) 
+        else (int_of_string number) 
       in
       match state with
       | X -> { x= Positive_int.(of_int ( (to_int accu.x) +n));
@@ -111,10 +111,9 @@ module Xyz = struct
     | 'Z'::rest | 'z'::rest -> 
         let new_accu = flush state accu number in 
         do_work Z new_accu "" rest
-    | c::rest -> do_work state accu (number^(String.of_char c)) rest
+    | c::rest -> do_work state accu (number^(String_ext.of_char c)) rest
     in
-    String.to_list_rev s 
+    String_ext.to_list s 
-    |> List.rev 
     |> do_work Null 
      { x=Positive_int.of_int 0 ; 
        y=Positive_int.of_int 0 ;

ocaml/Symmetry.mli

@@ -1,4 +1,4 @@
-type t = S | P | D | F | G | H | I | J | K | L with sexp
+type t = S | P | D | F | G | H | I | J | K | L [@@deriving sexp]
 
 (** Creatio from strings *)
 val to_string : t -> string
@@ -16,7 +16,7 @@ module Xyz :
       x : Qptypes.Positive_int.t;
       y : Qptypes.Positive_int.t;
       z : Qptypes.Positive_int.t;
-    } with sexp
+    } [@@deriving sexp]
 
     (** The string format contains the powers of x,y and z in a
         format like "x2z3" *)

ocaml/TaskServer.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 
@@ -63,7 +63,7 @@ let bind_socket ~socket_type ~socket ~port =
         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) )
+      | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_sec 1. ; loop (i-1) )
       | other_exception -> raise other_exception
   in loop 60
 

ocaml/To_md5.ml

@@ -1,5 +1,5 @@
-open Core.Std;;
+open Qptypes
-open Qptypes;;
+open Sexplib
 
 let to_md5 sexp_of_t t =
   sexp_of_t t

ocaml/_tags

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

ocaml/create_git_sha1.sh

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

ocaml/qp_basis_clean.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let filenames =
   let dir_name = Qpackage.root^"/data/basis/"

ocaml/qp_create_ezfio_from_xyz.ml

@@ -1,6 +1,6 @@
 open Qputils
 open Qptypes
-open Core.Std
+open Core
 
 let spec =
   let open Command.Spec in

ocaml/qp_create_guess.ml

@@ -1,6 +1,6 @@
 open Qputils
 open Qptypes
-open Core.Std
+open Core
 
 let run ~multiplicity ezfio_file =
   if (not (Sys.file_exists_exn ezfio_file)) then

ocaml/qp_find_pi_space.ml

@@ -0,0 +1,109 @@
+open Core
+open Qputils
+open Qptypes
+
+let run ?(sym="None") ezfio_filename =
+  Ezfio.set_file ezfio_filename ;
+
+  let aos = 
+    match Input.Ao_basis.read () with
+    | Some aos -> aos
+    | None -> failwith "Unable to read AOs"
+  and mos = 
+    match Input.Mo_basis.read () with
+    | Some mos -> mos
+    | None -> failwith "Unable to read MOs"
+  in
+  let rec find_power symmetry accu = function
+    | -1 -> accu
+    | i  -> let new_accu = 
+              if aos.Input.Ao_basis.ao_power.(i) = symmetry then (i::accu) 
+              else accu
+           in
+           find_power symmetry new_accu (i-1)
+  and n =
+    (AO_number.to_int aos.Input.Ao_basis.ao_num) 
+  and m =
+    (MO_number.to_int mos.Input.Mo_basis.mo_tot_num) 
+  and one  = Positive_int.of_int 1
+  and zero =  Positive_int.of_int 0
+  in
+
+  (* Indices of the px, py and pz AOs *)
+  let x_indices =
+    find_power Symmetry.Xyz.{x=one  ; y=zero ; z=zero} [] (n-1)
+  and y_indices =
+    find_power Symmetry.Xyz.{x=zero ; y=one  ; z=zero} [] (n-1)
+  and z_indices =
+    find_power Symmetry.Xyz.{x=zero ; y=zero ; z=one } [] (n-1)
+  in
+
+  (* Compute the relative weight of each MO on the px, py, pz spaces *)
+  let compute_weight mo_i list_aos =
+    let num =
+       List.fold_left ~f:(fun s i -> s +. (MO_coef.to_float @@ mos.Input.Mo_basis.mo_coef.(mo_i).(i)) ** 2.) ~init:0. list_aos
+    and denom =
+       Array.fold ~f:(fun s x -> s +. (MO_coef.to_float x) ** 2.) ~init:0. mos.Input.Mo_basis.mo_coef.(mo_i) 
+    in
+    num /. denom
+  in
+  let result = 
+    Array.init ~f:(fun mo_i -> 
+      (mo_i+1,
+      compute_weight mo_i x_indices,
+      compute_weight mo_i y_indices,
+      compute_weight mo_i z_indices) ) m
+    |> Array.to_list
+  in
+
+  let pi, sigma =
+    let rec aux test_xyz (accu_pi: int list) (accu_sigma: int list) = function
+    | [] -> (List.rev accu_pi, List.rev accu_sigma)
+    | (mo_i,x,y,z)::rest ->
+        if test_xyz (x,y,z) then
+           aux test_xyz (mo_i::accu_pi) accu_sigma rest
+         else
+           aux test_xyz accu_pi (mo_i::accu_sigma) rest
+    in
+    match sym with
+    | "x" | "X" -> aux (fun (x,y,z) -> (x>y && x>z)) [] [] result
+    | "y" | "Y" -> aux (fun (x,y,z) -> (y>x && y>z)) [] [] result
+    | "z" | "Z" -> aux (fun (x,y,z) -> (z>x && z>y)) [] [] result
+    | _         -> ([],[])
+  in
+
+  match sym with
+  | "x" | "X"  | "y" | "Y"  | "z" | "Z" -> 
+      begin
+        Printf.printf "Pi: [";
+        List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) pi;
+        Printf.printf "\b]\n\nSigma: [";
+        List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) sigma;
+        Printf.printf "\b]\n"
+      end
+  | _ -> List.iter ~f:(fun (mo_i,x,y,z) -> Printf.printf "%d: (%f,%f,%f)\n" mo_i x y z) result
+  
+
+
+
+
+let spec =
+  let open Command.Spec in
+  empty
+  +> flag "sym" (optional string) ~doc:"{x,y,z} Axis perpendicular to the plane"
+  +> anon ("ezfio_filename" %: string)
+
+
+let command =
+    Command.basic
+    ~summary: "Quantum Package command"
+    ~readme:(fun () ->
+     "Find all the pi molecular orbitals to create a pi space.
+     ")
+    spec
+    (fun sym ezfio_filename () -> run ?sym ezfio_filename )
+
+
+let () =
+    Command.run command
+

ocaml/qp_overlap_of_wf.ml

@@ -1,6 +1,12 @@
+(** 
+ * Computes the overlap <Psi_0 | Psi_1> where both Psi_0 and Psi_1 are truncated in the set
+ * of common determinants and normalized
+ *)
+
 open Input_determinants_by_hand
 open Qptypes
 
+
 let () =
   let ezfio, ezfio' =
     try
@@ -42,16 +48,16 @@ let () =
   let overlap wf wf' =
     let result, norm, norm' = 
       Hashtbl.fold (fun k c (accu,norm,norm') ->
-        let c' =
+        let (c',c) =
-          try  Hashtbl.find wf' k 
+          try  (Hashtbl.find wf' k, c)
-          with Not_found -> 0.
+          with Not_found -> (0.,0.)
         in
         (accu +. c *. c' , 
         norm +. c *. c  , 
         norm'+. c'*. c' ) 
       ) wf (0.,0.,0.)
     in 
-    result /. (norm *. norm')
+    result /. (sqrt (norm *. norm'))
   in
 
   let wf, wf' = 
@@ -62,5 +68,6 @@ let () =
   let o = 
     overlap wf wf'
   in
-  print_float (abs_float o)
+  print_float (abs_float o) ;
+  print_newline ()
   

ocaml/qp_print.ml

@@ -1,6 +1,6 @@
 open Qputils;;
 open Qptypes;;
-open Core.Std;;
+open Core;;
 
 type input_action =
   | Basis

ocaml/qp_print_basis.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 let basis () =

ocaml/qp_run.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qputils
 
 (* Environment variables :
@@ -132,7 +132,7 @@ let run slave exe ezfio_file =
     Sys.remove qp_run_address_filename;
 
   let duration = Time.diff (Time.now()) time_start 
-  |> Core.Span.to_string in
+  |> Time.Span.to_string in
   Printf.printf "Wall time : %s\n\n" duration;
   if (exit_code <> 0) then
     exit exit_code

ocaml/qp_set_mo_class.ml

@@ -1,6 +1,6 @@
 open Qputils
 open Qptypes
-open Core.Std
+open Core
 
 (*
  * Command-line arguments
@@ -210,7 +210,7 @@ let get () =
   in
 
   let mo_tot_num =
-    MO_number.to_int data.Input_mo_basis.mo_tot_num
+    MO_number.to_int data.Input.Mo_basis.mo_tot_num
   in
 
   let n_int =
@@ -244,7 +244,7 @@ let get () =
     | (MO_class.Deleted  _) :: rest ->
         work ~del:(Printf.sprintf "%s,%d" del   i) ~inact ~act  ~virt ~core (i+1) rest
   in
-  work 1 (Array.to_list data.Input_mo_basis.mo_class)
+  work 1 (Array.to_list data.Input.Mo_basis.mo_class)
 
 
 

ocaml/qptypes_generator.ml

@@ -1,4 +1,10 @@
-open Core.Std;;
+let global_replace x = 
+  x
+  |> Str.global_replace (Str.regexp "Float.to_string") "string_of_float"
+  |> Str.global_replace (Str.regexp "Float.of_string") "float_of_string"
+  |> Str.global_replace (Str.regexp "Int.to_string") "string_of_int"
+  |> Str.global_replace (Str.regexp "Int.of_string") "int_of_string"
+  |> Str.global_replace (Str.regexp "String.\\(to\\|of\\)_string") ""
 
 let input_data = "
 * Positive_float : float  
@@ -118,8 +124,12 @@ let input_data = "
 
 * MD5 : string
   assert ((String.length x) = 32);
-  assert (String.fold x ~init:true ~f:(fun accu x ->
+  assert (
-     accu && (x < 'g')));
+    let a =
+      Array.init (String.length x) (fun i -> x.[i])
+    in
+    Array.fold_left (fun accu x -> accu && (x < 'g')) true a
+    );
 
 * Rst_string : string
 
@@ -127,7 +137,7 @@ let input_data = "
   assert (x <> \"\") ;
 
 "
-;;
+
 
 let input_ezfio = "
 * MO_number : int
@@ -156,25 +166,25 @@ let input_ezfio = "
   More than 10 billion of determinants
 
 "
-;;
+
 
 let untouched = "
 module MO_guess : sig
-  type t with sexp
+  type t [@@deriving sexp]
   val to_string : t -> string
   val of_string : string -> t
 end = struct
   type t = 
   | Huckel
   | HCore
-  with sexp
+  [@@deriving sexp]
 
   let to_string = function
   | Huckel -> \"Huckel\"
   | HCore  -> \"HCore\"
 
   let of_string  s = 
-    match (String.lowercase s) with
+    match (String.lowercase_ascii s) with
     | \"huckel\" -> Huckel
     | \"hcore\"  -> HCore
     | _ -> raise (Invalid_argument (\"Wrong Guess type : \"^s))
@@ -182,7 +192,7 @@ end = struct
 end
 
 module Disk_access : sig
-  type t with sexp
+  type t [@@deriving sexp]
   val to_string : t -> string
   val of_string : string -> t
 end = struct
@@ -190,14 +200,14 @@ end = struct
   | Read
   | Write
   | None
-  with sexp
+  [@@deriving sexp]
 
   let to_string = function
   | Read   -> \"Read\"
   | Write  -> \"Write\"
   | None   -> \"None\"
   let of_string  s = 
-    match (String.lowercase s) with
+    match (String.lowercase_ascii s) with
     | \"read\"  -> Read
     | \"write\" -> Write
     | \"none\"  -> None
@@ -206,62 +216,63 @@ end = struct
 end
 "
 
-;;
+
 
 let template = format_of_string "
 module %s : sig
-  type t with sexp
+  type t [@@deriving sexp]
   val to_%s : t -> %s
   val of_%s : %s %s -> t
   val to_string : t -> string
 end = struct
-  type t = %s with sexp
+  type t = %s [@@deriving sexp]
   let to_%s x = x
   let of_%s %s x = ( %s x )
   let to_string x = %s.to_string x
 end
 
 "
-;;
+
 
 let parse_input input=
-  print_string "open Core.Std;;\nlet warning = print_string;;\n" ;
+  print_string "open Sexplib.Std\nlet warning = print_string\n" ;
   let rec parse result = function
     | [] -> result
     | ( "" , ""   )::tail -> parse result tail
     | ( t  , text )::tail -> 
         let name,typ,params,params_val = 
-          match String.split ~on:':' t with
+          match String_ext.split ~on:':' t with
           | [name;typ] -> (name,typ,"","")
           | name::typ::params::params_val -> (name,typ,params,
-            (String.concat params_val ~sep:":") )
+            (String.concat ":" params_val) )
           | _ -> assert false
         in
-        let typ  = String.strip typ
+        let typ  = String_ext.strip typ
-        and name = String.strip name in
+        and name = String_ext.strip name in
         let typ_cap = String.capitalize typ in
         let newstring = Printf.sprintf template name typ typ typ params_val typ typ 
-          typ typ params ( String.strip text ) typ_cap
+          typ typ params ( String_ext.strip text ) typ_cap
         in
         List.rev (parse (newstring::result) tail )
   in
-     String.split ~on:'*' input
+     String_ext.split ~on:'*' input
-  |> List.map ~f:(String.lsplit2_exn ~on:'\n') 
+  |> List.map (String_ext.lsplit2_exn ~on:'\n') 
   |> parse []
-  |> String.concat 
+  |> String.concat  ""
+  |> global_replace
   |> print_string
-;;
+
 
 
 let ezfio_template = format_of_string "
 module %s : sig
-  type t with sexp
+  type t [@@deriving sexp]
   val to_%s : t -> %s
   val get_max : unit -> %s
   val of_%s : ?min:%s -> ?max:%s -> %s -> t
   val to_string : t -> string
 end = struct
-  type t = %s with sexp
+  type t = %s [@@deriving sexp]
   let to_string x = %s.to_string x
   let get_max () =
     if (Ezfio.has_%s ()) then
@@ -287,24 +298,24 @@ end = struct
     end
 end
 "
-;;
+
 
 let parse_input_ezfio input=
   let parse s = 
     match (
-      String.split s ~on:'\n'
+      String_ext.split s ~on:'\n'
-      |> List.filter ~f:(fun x -> (String.strip x) <> "")
+      |> List.filter (fun x -> (String_ext.strip x) <> "")
     ) with
     | [] -> ""
     | a :: b :: c :: d :: [] ->
       begin
-        let (name,typ) = String.lsplit2_exn ~on:':' a
+        let (name,typ) = String_ext.lsplit2_exn ~on:':' a
         and ezfio_func = b
-        and (min, max) = String.lsplit2_exn ~on:':' c
+        and (min, max) = String_ext.lsplit2_exn ~on:':' c
         and msg = d
         in 
         let (name, typ, ezfio_func, min, max, msg) = 
-        match (List.map [ name ; typ ; ezfio_func ; min ; max ; msg ] ~f:String.strip) with
+        match List.map String_ext.strip [ name ; typ ; ezfio_func ; min ; max ; msg ] with
         | [ name ; typ ; ezfio_func ; min ; max ; msg ] -> (name, typ, ezfio_func, min, max, msg)
         | _ -> assert false
         in
@@ -314,16 +325,17 @@ let parse_input_ezfio input=
       end
     | _ -> failwith "Error in input_ezfio"
   in
-     String.split ~on:'*' input
+     String_ext.split ~on:'*' input
-  |> List.map ~f:parse
+  |> List.map parse
-  |> String.concat 
+  |> String.concat ""
+  |> global_replace
   |> print_string
-;;
+
 
 let () = 
   parse_input input_data ;
   parse_input_ezfio input_ezfio;
-  print_endline untouched;
+  print_endline untouched
 
 
 

ocaml/test_basis.ml

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

ocaml/test_gto.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 let test_prim () = 

ocaml/test_message.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let () =
   Message.of_string "new_job ao_integrals tcp://127.0.0.1 inproc://ao_ints:12345"

ocaml/test_molecule.ml

@@ -1,4 +1,4 @@
-open Core.Std ;;
+open Core ;;
 open Qptypes ;;
 
 let test_molecule () =

ocaml/test_pseudo.ml

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

ocaml/test_queuing_system.ml

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

ocaml/test_symmetry.ml

@@ -1,15 +1,13 @@
-open Core.Std
 open Qputils
 open Qptypes
 open Symmetry
 
 let () =
   "SPDFGHIJKL"
-  |> String.to_list_rev
+  |> String_ext.to_list
-  |> List.rev
+  |> List.map of_char
-  |> List.map ~f:of_char
+  |> List.map Xyz.of_symmetry
-  |> List.map ~f:Xyz.of_symmetry
+  |> List.iter (fun x -> List.iter (fun y -> Xyz.to_string y |> print_endline) x ;
-  |> List.iter ~f:(fun x -> List.iter x ~f:(fun y -> Xyz.to_string y |> print_endline) ;
      print_newline ();)
 
 

ocaml/test_task_server.ml

@@ -1,5 +1,3 @@
-open Core
-
 let () =
   TaskServer.run 12345
 

plugins/CAS_SD_ZMQ/selection.irp.f

@@ -1,28 +1,6 @@
 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
-  integer :: ii
-  ii = l-mo_integrals_cache_min
-  ii = ior(ii, k-mo_integrals_cache_min)
-  ii = ior(ii, j-mo_integrals_cache_min)
-  ii = ior(ii, i-mo_integrals_cache_min)
-  if (iand(ii, -64) /= 0) then
-    integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
-  else
-    ii = l-mo_integrals_cache_min
-    ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
-    ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
-    ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
-    integral8 = mo_integrals_cache(ii)
-  endif
-end function
-
-
 BEGIN_PROVIDER [ integer(1), psi_phasemask, (N_int*bit_kind_size, 2, N_det)]
   use bitmasks
   implicit none
@@ -287,7 +265,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   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
+  double precision, external :: get_phase_bi, mo_bielec_integral
   
   integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
   integer, parameter :: turn2(2) = (/2,1/) 
@@ -300,7 +278,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
       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 = mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2, p1, h1, h2)
       hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2)
       vect(:, puti) += hij * coefs
     end do
@@ -313,7 +291,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
       puti = p(j, sp)
       if(bannedOrb(puti)) cycle
       pmob = p(turn2(j), sp)
-      hij = integral8(pfix, pmob, hfix, hmob)
+      hij = mo_bielec_integral(pfix, pmob, hfix, hmob)
       hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix)
       vect(:, puti) += hij * coefs
     end do
@@ -325,7 +303,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
       p2 = p(2,sfix)
       h1 = h(1,sfix)
       h2 = h(2,sfix)
-      hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2))
+      hij = (mo_bielec_integral(p1,p2,h1,h2) - mo_bielec_integral(p2,p1,h1,h2))
       hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2)
       vect(:, puti) += hij * coefs
     end if
@@ -348,7 +326,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   logical :: ok, lbanned(mo_tot_num)
   integer(bit_kind) :: det(N_int, 2)
   double precision :: hij
-  double precision, external :: get_phase_bi, integral8
+  double precision, external :: get_phase_bi, mo_bielec_integral
   
   lbanned = bannedOrb
   sh = 1
@@ -366,13 +344,13 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
     
     do i=1,hole-1
       if(lbanned(i)) cycle
-      hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole))
+      hij = (mo_bielec_integral(p1, p2, i, hole) - mo_bielec_integral(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 = (mo_bielec_integral(p1, p2, hole, i) - mo_bielec_integral(p2, p1, hole, i))
       hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2)
       vect(:,i) += hij * coefs
     end do
@@ -384,7 +362,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
     p2 = p(1, sh)
     do i=1,mo_tot_num
       if(lbanned(i)) cycle
-      hij = integral8(p1, p2, i, hole)
+      hij = mo_bielec_integral(p1, p2, i, hole)
       hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2)
       vect(:,i) += hij * coefs
     end do
@@ -787,7 +765,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   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
+  double precision, external :: get_phase_bi, mo_bielec_integral
   
   integer :: i, j, tip, ma, mi, puti, putj
   integer :: h1, h2, p1, p2, i1, i2
@@ -822,7 +800,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
         if(ma == 1) then
           mat(:, putj, puti) += coefs * hij
         else
@@ -841,7 +819,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         h1 = h(1,1)
         h2 = h(1,2)
         
-        hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
+        hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
         mat(:, puti, putj) += coefs * hij
       end do
       end do
@@ -861,7 +839,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
         mat(:, puti, putj) += coefs * hij
       end do
       end do
@@ -875,7 +853,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = mo_bielec_integral(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
@@ -886,7 +864,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2)
         mat(:, puti, putj) += coefs * hij
       end if
     end if
@@ -905,7 +883,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   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
+  double precision, external :: get_phase_bi, mo_bielec_integral
   
   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
@@ -944,12 +922,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       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)
+        hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(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)
+        hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(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
 
@@ -969,13 +947,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       !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)
+        hij = mo_bielec_integral(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)
+        hij = mo_bielec_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
         tmp_row2(:,puti) += hij * coefs
       end if
     end do
@@ -997,12 +975,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+          hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(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)
+          hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
           tmp_row(:,putj) += hij * coefs
         end do
 
@@ -1020,13 +998,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+          hij = mo_bielec_integral(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)
+          hij = mo_bielec_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
           tmp_row2(:,puti) += hij * coefs
         end if
       end do
@@ -1077,7 +1055,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   
   integer :: i, j, s, h1, h2, p1, p2, puti, putj
   double precision :: hij, phase
-  double precision, external :: get_phase_bi, integral8
+  double precision, external :: get_phase_bi, mo_bielec_integral
   logical :: ok
   
   integer :: bant
@@ -1096,7 +1074,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
           call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
           call i_h_j(gen, det, N_int, hij)
         else
-          hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
+          hij = mo_bielec_integral(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
@@ -1114,7 +1092,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
           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)
+          hij = (mo_bielec_integral(p1, p2, puti, putj) -  mo_bielec_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2)
         end if
         mat(:, puti, putj) += coefs(:) * hij
       end do

plugins/FOBOCI/diag_fock_inactiv_virt.irp.f

@@ -29,7 +29,7 @@ subroutine diag_inactive_virt_and_update_mos
 
 
  label = "Canonical"
- call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
+ call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.)
  soft_touch mo_coef
 
 
@@ -76,7 +76,7 @@ subroutine diag_inactive_virt_new_and_update_mos
 
 
  label = "Canonical"
- call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
+ call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.)
  soft_touch mo_coef
 
 

plugins/FOBOCI/new_approach.irp.f

@@ -655,7 +655,7 @@ subroutine new_approach
  integer :: sign
  sign = -1
  
- call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign)
+ call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign,.true.)
  soft_touch mo_coef
  call save_mos
  

plugins/FOBOCI/routines_foboci.irp.f

@@ -449,7 +449,7 @@ subroutine save_osoci_natural_mos
 
  label = "Natural"
  
- call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
+ call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.)
 !if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Write" )then
 ! disk_access_ao_integrals = "Read"
 ! touch disk_access_ao_integrals
@@ -584,7 +584,7 @@ subroutine set_osoci_natural_mos
  enddo
 
  label = "Natural"
- call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
+ call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.)
  soft_touch mo_coef
  deallocate(tmp,occ)
 

plugins/FourIdx/four_index.irp.f

@@ -0,0 +1,180 @@
+subroutine four_index_transform(map_a,map_c,matrix_B,LDB,            &
+      i_start, j_start, k_start, l_start,                            &
+      i_end  , j_end  , k_end  , l_end  ,                            &
+      a_start, b_start, c_start, d_start,                            &
+      a_end  , b_end  , c_end  , d_end  )
+  implicit none
+  use map_module
+  use mmap_module
+  BEGIN_DOC
+! Performs a four-index transformation of map_a(N^4) into map_c(M^4) using b(NxM)
+! C_{abcd} = \sum_{ijkl} A_{ijkl}.B_{ia}.B_{jb}.B_{kc}.B_{ld}
+! Loops run over *_start->*_end
+  END_DOC
+  type(map_type), intent(in)     :: map_a
+  type(map_type), intent(inout)  :: map_c
+  integer, intent(in)            :: LDB
+  double precision, intent(in)   :: matrix_B(LDB,*)
+  integer, intent(in)            :: i_start, j_start, k_start, l_start
+  integer, intent(in)            :: i_end  , j_end  , k_end  , l_end
+  integer, intent(in)            :: a_start, b_start, c_start, d_start
+  integer, intent(in)            :: a_end  , b_end  , c_end  , d_end
+
+  double precision, allocatable  :: T(:,:,:), U(:,:,:), V(:,:,:)
+  integer                        :: i_max, j_max, k_max, l_max
+  integer                        :: i_min, j_min, k_min, l_min
+  integer                        :: i, j, k, l
+  integer                        :: a, b, c, d
+  double precision, external     :: get_ao_bielec_integral
+  integer(key_kind)              :: idx
+  real(integral_kind)            :: tmp
+  integer(key_kind), allocatable :: key(:)
+  real(integral_kind), allocatable :: value(:)
+
+  ASSERT (k_start == i_start)
+  ASSERT (l_start == j_start)
+  ASSERT (a_start == c_start)
+  ASSERT (b_start == d_start)
+
+  i_min = min(i_start,a_start)
+  i_max = max(i_end  ,a_end  )
+  j_min = min(j_start,b_start)
+  j_max = max(j_end  ,b_end  )
+  k_min = min(k_start,c_start)
+  k_max = max(k_end  ,c_end  )
+  l_min = min(l_start,d_start)
+  l_max = max(l_end  ,d_end  )
+
+  ASSERT (0 < i_max)
+  ASSERT (0 < j_max)
+  ASSERT (0 < k_max)
+  ASSERT (0 < l_max)
+  ASSERT (LDB >= i_max)
+  ASSERT (LDB >= j_max)
+  ASSERT (LDB >= k_max)
+  ASSERT (LDB >= l_max)
+
+  ! Create a temporary memory-mapped file
+  integer                        :: fd
+  type(c_ptr)                    :: c_pointer
+  integer*8, pointer             :: a_array(:,:,:)
+  call mmap(trim(ezfio_filename)//'/work/four_idx',                  &
+      (/ 4_8,int(i_end-i_start+1,8),int(j_end-j_start+1,8),int(k_end-k_start+1,8), int(l_end-l_start+1,8) /), 8, fd, .False., c_pointer)
+  call c_f_pointer(c_pointer, a_array, (/ 4, (i_end-i_start+1)*(j_end-j_start+1)*(k_end-k_start+1), l_end-l_start+1 /))
+
+
+  !$OMP PARALLEL DEFAULT(NONE) SHARED(a_array,c_pointer,fd,          &
+      !$OMP  a_start,a_end,b_start,b_end,c_start,c_end,d_start,d_end,&
+      !$OMP  i_start,i_end,j_start,j_end,k_start,k_end,l_start,l_end,&
+      !$OMP  i_min,i_max,j_min,j_max,k_min,k_max,l_min,l_max,        &
+      !$OMP  map_a,map_c,matrix_B)                                   &
+      !$OMP  PRIVATE(key,value,T,U,V,i,j,k,l,idx,   &
+      !$OMP  a,b,c,d,tmp)
+  allocate( key(i_max*j_max*k_max), value(i_max*j_max*k_max) )
+  allocate( U(a_start:a_end, c_start:c_end, b_start:b_end) )
+
+
+  !$OMP DO SCHEDULE(dynamic,4)
+  do l=l_start,l_end
+    a = 1
+    do j=j_start,j_end
+      do k=k_start,k_end
+        do i=i_start,i_end
+          call bielec_integrals_index(i,j,k,l,idx)
+          call map_get(map_a,idx,tmp)
+          if (tmp /= 0.d0) then
+            a = a+1
+            a_array(1,a,l-l_start+1) = i
+            a_array(2,a,l-l_start+1) = j
+            a_array(3,a,l-l_start+1) = k
+            a_array(4,a,l-l_start+1) = transfer(dble(tmp), 1_8)
+          endif
+        enddo
+      enddo
+    enddo
+    a_array(1,1,l-l_start+1) = a
+    print *,  l
+  enddo
+  !$OMP END DO
+
+  !$OMP DO SCHEDULE(dynamic)
+  do d=d_start,d_end
+    U = 0.d0
+    do l=l_start,l_end
+      if (dabs(matrix_B(l,d)) < 1.d-10) then
+        cycle
+      endif
+      print *,  d, l
+
+      allocate( T(i_start:i_end, k_start:k_end, j_start:j_end), &
+                V(a_start:a_end, k_start:k_end, j_start:j_end) )
+
+      T = 0.d0
+      do a=2,a_array(1,1,l-l_start+1)
+        i = a_array(1,a,l-l_start+1)
+        j = a_array(2,a,l-l_start+1)
+        k = a_array(3,a,l-l_start+1)
+        T(i, k,j) = transfer(a_array(4,a,l-l_start+1), 1.d0)
+      enddo
+
+      call DGEMM('T','N', (a_end-a_start+1),                         &
+          (k_end-k_start+1)*(j_end-j_start+1),                       &
+          (i_end-i_start+1), 1.d0,                                   &
+          matrix_B(i_start,a_start), size(matrix_B,1),               &
+          T(i_start,k_start,j_start), size(T,1),  0.d0,              &
+          V(a_start,k_start,j_start), size(V, 1) )
+
+      deallocate(T)
+      allocate( T(a_start:a_end, k_start:k_end, b_start:d) )
+
+      call DGEMM('N','N', (a_end-a_start+1)*(k_end-k_start+1),       &
+              (b_end-b_start+1),                                     &
+              (j_end-j_start+1), 1.d0,                               &
+              V(a_start,k_start,j_start), size(V,1)*size(V,2),       &
+              matrix_B(j_start,b_start), size(matrix_B,1),0.d0,      &
+              T(a_start,k_start,b_start), size(T,1)*size(T,2) )
+
+      deallocate(V)
+
+      do b=b_start,b_end
+        call DGEMM('N','N', (a_end-a_start+1), (c_end-c_start+1),    &
+            (k_end-k_start+1), matrix_B(l, d),                   &
+            T(a_start,k_start,b), size(T,1),                     &
+            matrix_B(k_start,c_start), size(matrix_B,1), 1.d0,   &
+            U(a_start,c_start,b), size(U,1) )
+      enddo
+
+      deallocate(T)
+
+    enddo
+
+    idx = 0_8
+    do b=b_start,b_end
+      do c=c_start,c_end
+        do a=a_start,a_end
+          if (dabs(U(a,c,b)) < 1.d-15) then
+            cycle
+          endif
+          idx = idx+1_8
+          call bielec_integrals_index(a,b,c,d,key(idx))
+          value(idx) = U(a,c,b)
+        enddo
+      enddo
+    enddo
+
+    !$OMP CRITICAL
+    call map_append(map_c, key, value, idx) 
+    call map_sort(map_c)
+    !$OMP END CRITICAL
+
+
+  enddo
+  !$OMP END DO
+
+  deallocate(key,value)
+  !$OMP END PARALLEL
+
+  call munmap( &
+      (/ 4_8,int(i_end-i_start+1,8),int(j_end-j_start+1,8),int(k_end-k_start+1,8), int(l_end-l_start+1,8) /), 8, fd, c_pointer)
+
+end

plugins/FourIdx/four_index_sym.irp.f

@@ -0,0 +1,277 @@
+subroutine four_index_transform_sym(map_a,map_c,matrix_B,LDB,            &
+      i_start, j_start, k_start, l_start,                            &
+      i_end  , j_end  , k_end  , l_end  ,                            &
+      a_start, b_start, c_start, d_start,                            &
+      a_end  , b_end  , c_end  , d_end  )
+  implicit none
+  use map_module
+  use mmap_module
+  BEGIN_DOC
+! Performs a four-index transformation of map_a(N^4) into map_c(M^4) using b(NxM)
+! C_{abcd} = \sum_{ijkl} A_{ijkl}.B_{ia}.B_{jb}.B_{kc}.B_{ld}
+! Loops run over *_start->*_end
+  END_DOC
+  type(map_type), intent(in)     :: map_a
+  type(map_type), intent(inout)  :: map_c
+  integer, intent(in)            :: LDB
+  double precision, intent(in)   :: matrix_B(LDB,*)
+  integer, intent(in)            :: i_start, j_start, k_start, l_start
+  integer, intent(in)            :: i_end  , j_end  , k_end  , l_end
+  integer, intent(in)            :: a_start, b_start, c_start, d_start
+  integer, intent(in)            :: a_end  , b_end  , c_end  , d_end
+
+  double precision, allocatable  :: T(:,:), U(:,:,:), V(:,:)
+  double precision, allocatable  :: T2d(:,:), V2d(:,:)
+  integer                        :: i_max, j_max, k_max, l_max
+  integer                        :: i_min, j_min, k_min, l_min
+  integer                        :: i, j, k, l, ik, ll
+  integer                        :: a, b, c, d
+  double precision, external     :: get_ao_bielec_integral
+  integer*8                      :: ii
+  integer(key_kind)              :: idx
+  real(integral_kind)            :: tmp
+  integer(key_kind), allocatable :: key(:)
+  real(integral_kind), allocatable :: value(:)
+  integer*8, allocatable         :: l_pointer(:)
+
+  ASSERT (k_start == i_start)
+  ASSERT (l_start == j_start)
+  ASSERT (a_start == c_start)
+  ASSERT (b_start == d_start)
+
+  i_min = min(i_start,a_start)
+  i_max = max(i_end  ,a_end  )
+  j_min = min(j_start,b_start)
+  j_max = max(j_end  ,b_end  )
+  k_min = min(k_start,c_start)
+  k_max = max(k_end  ,c_end  )
+  l_min = min(l_start,d_start)
+  l_max = max(l_end  ,d_end  )
+
+  ASSERT (0 < i_max)
+  ASSERT (0 < j_max)
+  ASSERT (0 < k_max)
+  ASSERT (0 < l_max)
+  ASSERT (LDB >= i_max)
+  ASSERT (LDB >= j_max)
+  ASSERT (LDB >= k_max)
+  ASSERT (LDB >= l_max)
+
+  ! Create a temporary memory-mapped file
+  integer                        :: fd
+  type(c_ptr)                    :: c_pointer
+  integer*8, pointer             :: a_array(:)
+  call mmap(trim(ezfio_filename)//'/work/four_idx',                  &
+      (/ 12_8 * map_a % n_elements /), 8, fd, .False., c_pointer)
+  call c_f_pointer(c_pointer, a_array, (/ 12_8 * map_a % n_elements /))
+
+  allocate(l_pointer(l_start:l_end+1), value((i_max*k_max)) )
+  ii = 1_8
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,ik,idx) 
+  do l=l_start,l_end
+    !$OMP SINGLE
+    l_pointer(l) = ii
+    !$OMP END SINGLE
+    do j=j_start,j_end
+      !$OMP DO SCHEDULE(static,1)
+      do k=k_start,k_end
+        do i=i_start,k
+          ik = (i-i_start+1) + ishft( (k-k_start)*(k-k_start+1), -1 )
+          call bielec_integrals_index(i,j,k,l,idx)
+          call map_get(map_a,idx,value(ik))
+        enddo
+      enddo
+      !$OMP END DO
+
+      !$OMP SINGLE
+      ik=0
+      do k=k_start,k_end
+        do i=i_start,k
+          ik = ik+1
+          tmp=value(ik)
+          if (tmp /= 0.d0) then
+            a_array(ii) = ik
+            ii = ii+1_8
+            a_array(ii) = j
+            ii = ii+1_8
+            a_array(ii) = transfer(dble(tmp), 1_8)
+            ii = ii+1_8
+          endif
+        enddo
+      enddo
+      !$OMP END SINGLE
+    enddo
+  enddo
+  !$OMP SINGLE
+  l_pointer(l_end+1) = ii
+  !$OMP END SINGLE
+  !$OMP END PARALLEL  
+  deallocate(value)
+
+!INPUT DATA
+!open(unit=10,file='INPUT',form='UNFORMATTED')
+!write(10) i_start, j_start, i_end, j_end
+!write(10) a_start, b_start, a_end, b_end
+!write(10) LDB, mo_tot_num
+!write(10) matrix_B(1:LDB,1:mo_tot_num)
+!idx=size(a_array)
+!write(10) idx
+!write(10) a_array
+!write(10) l_pointer
+!close(10)
+!open(unit=10,file='OUTPUT',form='FORMATTED')
+! END INPUT DATA
+
+
+  !$OMP PARALLEL DEFAULT(NONE) SHARED(a_array,c_pointer,fd,          &
+      !$OMP  a_start,a_end,b_start,b_end,c_start,c_end,d_start,d_end,&
+      !$OMP  i_start,i_end,j_start,j_end,k_start,k_end,l_start,l_end,&
+      !$OMP  i_min,i_max,j_min,j_max,k_min,k_max,l_min,l_max,        &
+      !$OMP  map_c,matrix_B,l_pointer)                         &
+      !$OMP  PRIVATE(key,value,T,U,V,i,j,k,l,idx,ik,ll,   &
+      !$OMP  a,b,c,d,tmp,T2d,V2d,ii)
+  allocate( key(i_max*j_max*k_max), value(i_max*j_max*k_max) )
+  allocate( U(a_start:a_end, c_start:c_end, b_start:b_end) )
+
+
+
+  allocate( T2d((i_end-i_start+1)*(k_end-k_start+2)/2, j_start:j_end), &
+            V2d((i_end-i_start+1)*(k_end-k_start+2)/2, b_start:b_end), &
+            V(i_start:i_end, k_start:k_end), &
+            T(k_start:k_end, a_start:a_end))
+
+
+  !$OMP DO SCHEDULE(dynamic)
+  do d=d_start,d_end
+    U = 0.d0
+    do l=l_start,l_end
+      if (dabs(matrix_B(l,d)) < 1.d-10) then
+        cycle
+      endif
+      
+      ii=l_pointer(l)
+      do j=j_start,j_end
+        ik=0
+        do k=k_start,k_end
+          do i=i_start,k
+            ik = ik+1
+            if ( (ik /= a_array(ii)).or.(j /= a_array(ii+1_8))  &
+                 .or.(ii >= l_pointer(l+1)) ) then
+              T2d(ik,j) = 0.d0
+            else
+              T2d(ik,j) = transfer(a_array(ii+2_8), 1.d0)
+              ii=ii+3_8
+            endif
+          enddo
+        enddo
+      enddo
+      call DGEMM('N','N', ishft( (i_end-i_start+1)*(i_end-i_start+2), -1),&
+          (d-b_start+1),                                             &
+          (j_end-j_start+1), 1.d0,                                   &
+          T2d(1,j_start), size(T2d,1),                               &
+          matrix_B(j_start,b_start), size(matrix_B,1),0.d0,          &
+          V2d(1,b_start), size(V2d,1) )
+
+      do b=b_start,d
+        ik = 0
+        do k=k_start,k_end
+          do i=i_start,k
+            ik = ik+1
+            V(i,k) = V2d(ik,b)
+          enddo
+        enddo
+
+!        T = 0.d0
+!        do a=a_start,b
+!          do k=k_start,k_end
+!            do i=i_start,k
+!              T(k,a) = T(k,a) + V(i,k)*matrix_B(i,a)
+!            enddo
+!            do i=k+1,i_end
+!              T(k,a) = T(k,a) + V(k,i)*matrix_B(i,a)
+!            enddo
+!          enddo
+!        enddo
+        call DSYMM('L','U', (k_end-k_start+1), (b-a_start+1),        &
+            1.d0,                                                    &
+            V(i_start,k_start), size(V,1),                           &
+            matrix_B(i_start,a_start), size(matrix_B,1),0.d0,        &
+            T(k_start,a_start), size(T,1) )
+
+!        do c=c_start,b
+!          do a=a_start,c
+!            do k=k_start,k_end
+!              U(a,c,b) = U(a,c,b) + T(k,a)*matrix_B(k,c)*matrix_B(l,d)
+!            enddo
+!          enddo
+!        enddo
+        call DGEMM('T','N', (b-a_start+1), (b-c_start+1),            &
+            (k_end-k_start+1), matrix_B(l, d),                       &
+            T(k_start,a_start), size(T,1),                           &
+            matrix_B(k_start,c_start), size(matrix_B,1), 1.d0,       &
+            U(a_start,c_start,b), size(U,1) )
+!        do c=b+1,c_end
+!          do a=a_start,b
+!            do k=k_start,k_end
+!              U(a,c,b) = U(a,c,b) + T(k,a)*matrix_B(k,c)*matrix_B(l,d)
+!            enddo
+!          enddo
+!        enddo
+        if (b < b_end) then
+          call DGEMM('T','N', (b-a_start+1), (c_end-b),              &
+              (k_end-k_start+1), matrix_B(l, d),                     &
+              T(k_start,a_start), size(T,1),                         &
+              matrix_B(k_start,b+1), size(matrix_B,1), 1.d0,         &
+              U(a_start,b+1,b), size(U,1) )
+        endif
+      enddo
+
+    enddo
+
+    idx = 0_8
+    do b=b_start,d
+      do c=c_start,c_end
+        do a=a_start,min(b,c)
+          if (dabs(U(a,c,b)) < 1.d-15) then
+            cycle
+          endif
+          idx = idx+1_8
+          call bielec_integrals_index(a,b,c,d,key(idx))
+          value(idx) = U(a,c,b)
+        enddo
+      enddo
+    enddo
+
+    !$OMP CRITICAL
+    call map_append(map_c, key, value, idx) 
+    !$OMP END CRITICAL
+
+!WRITE OUTPUT
+! OMP CRITICAL
+!print *,  d
+!do b=b_start,d
+!  do c=c_start,c_end
+!    do a=a_start,min(b,c)
+!      if (dabs(U(a,c,b)) < 1.d-15) then
+!        cycle
+!      endif
+!      write(10,*) d,c,b,a,U(a,c,b)
+!    enddo
+!  enddo
+!enddo
+! OMP END CRITICAL
+!END WRITE OUTPUT
+
+
+  enddo
+  !$OMP END DO
+
+  deallocate(key,value,V,T)
+  !$OMP END PARALLEL
+  call map_sort(map_c)
+
+  call munmap( &
+      (/ 12_8 * map_a % n_elements /), 8, fd, c_pointer)
+  deallocate(l_pointer)
+
+end

plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_full Generators_full ZMQ 
+Perturbation Selectors_full Generators_full ZMQ FourIdx

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -351,12 +351,12 @@ subroutine get_first_tooth(computed, first_teeth)
 end subroutine
 
 
-BEGIN_PROVIDER [ integer, size_tbc ]
+BEGIN_PROVIDER [ integer*8, size_tbc ]
   implicit none
   BEGIN_DOC
 ! Size of the tbc array
   END_DOC
-  size_tbc = (comb_teeth+1)*N_det_generators + fragment_count*fragment_first
+  size_tbc = int((comb_teeth+1),8)*int(N_det_generators,8) + fragment_count*fragment_first
 END_PROVIDER
 
 subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
@@ -409,7 +409,8 @@ end subroutine
 
 subroutine add_comb(comb, computed, tbc, stbc, ct)
   implicit none
-  integer, intent(in) :: stbc, ct
+  integer*8, intent(in) :: stbc
+  integer, intent(in) :: ct
   double precision, intent(in) :: comb
   logical, intent(inout) :: computed(N_det_generators)
   integer, intent(inout) :: tbc(0:stbc)

plugins/Full_CI_ZMQ/run_selection_slave.irp.f

@@ -57,7 +57,6 @@ subroutine run_selection_slave(thread,iproc,energy)
     endif
 
     if(done .or. ctask == size(task_id)) then
-      ASSERT (.not.(buf%N == 0 .and. ctask > 0))
       do i=1, ctask
          call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
       end do

plugins/Full_CI_ZMQ/selection.irp.f

@@ -9,29 +9,6 @@ BEGIN_PROVIDER [ integer, fragment_count ]
 END_PROVIDER
 
 
-double precision function integral8(i,j,k,l)
-  implicit none
-  
-  integer, intent(in) :: i,j,k,l
-  double precision, external :: get_mo_bielec_integral
-  integer :: ii
-  ii = l-mo_integrals_cache_min
-  ii = ior(ii, k-mo_integrals_cache_min)
-  ii = ior(ii, j-mo_integrals_cache_min)
-  ii = ior(ii, i-mo_integrals_cache_min)
-  if (iand(ii, -64) /= 0) then
-    integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
-  else
-    ii = l-mo_integrals_cache_min
-    ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
-    ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
-    ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
-    integral8 = mo_integrals_cache(ii)
-  endif
-end function
-
-
-
 subroutine assert(cond, msg)
   character(*), intent(in) :: msg
   logical, intent(in) :: cond
@@ -135,7 +112,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   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
+  double precision, external :: get_phase_bi, mo_bielec_integral
   
   integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
   integer, parameter :: turn2(2) = (/2,1/) 
@@ -148,7 +125,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
       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 = mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2, p1, h1, h2)
       hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2)
       vect(:, puti) += hij * coefs
     end do
@@ -161,7 +138,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
       puti = p(j, sp)
       if(bannedOrb(puti)) cycle
       pmob = p(turn2(j), sp)
-      hij = integral8(pfix, pmob, hfix, hmob)
+      hij = mo_bielec_integral(pfix, pmob, hfix, hmob)
       hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix)
       vect(:, puti) += hij * coefs
     end do
@@ -173,7 +150,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
       p2 = p(2,sfix)
       h1 = h(1,sfix)
       h2 = h(2,sfix)
-      hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2))
+      hij = (mo_bielec_integral(p1,p2,h1,h2) - mo_bielec_integral(p2,p1,h1,h2))
       hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2)
       vect(:, puti) += hij * coefs
     end if
@@ -198,7 +175,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   logical, allocatable           :: lbanned(:)
   integer(bit_kind)              :: det(N_int, 2)
   double precision               :: hij
-  double precision, external     :: get_phase_bi, integral8
+  double precision, external     :: get_phase_bi, mo_bielec_integral
   
   allocate (lbanned(mo_tot_num))
   lbanned = bannedOrb
@@ -217,13 +194,13 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
     
     do i=1,hole-1
       if(lbanned(i)) cycle
-      hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole))
+      hij = (mo_bielec_integral(p1, p2, i, hole) - mo_bielec_integral(p2, p1, i, hole))
       hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2)
       vect(1:N_states,i) += hij * coefs(1:N_states)
     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 = (mo_bielec_integral(p1, p2, hole, i) - mo_bielec_integral(p2, p1, hole, i))
       hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2)
       vect(1:N_states,i) += hij * coefs(1:N_states)
     end do
@@ -235,7 +212,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
     p2 = p(1, sh)
     do i=1,mo_tot_num
       if(lbanned(i)) cycle
-      hij = integral8(p1, p2, i, hole)
+      hij = mo_bielec_integral(p1, p2, i, hole)
       hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2)
       vect(1:N_states,i) += hij * coefs(1:N_states)
     end do
@@ -442,37 +419,82 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
       fullinteresting(0) = 0
       
       do ii=1,preinteresting(0)
-        i = preinteresting(ii)
+        select case (N_int)
-        mobMask(1,1) = iand(negMask(1,1), preinteresting_det(1,1,ii))
+          case (1)
-        mobMask(1,2) = iand(negMask(1,2), preinteresting_det(1,2,ii))
+            mobMask(1,1) = iand(negMask(1,1), preinteresting_det(1,1,ii))
-        nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+            mobMask(1,2) = iand(negMask(1,2), preinteresting_det(1,2,ii))
-        do j=2,N_int
+            nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
-          mobMask(j,1) = iand(negMask(j,1), preinteresting_det(j,1,ii))
+          case (2)
-          mobMask(j,2) = iand(negMask(j,2), preinteresting_det(j,2,ii))
+            mobMask(1:2,1) = iand(negMask(1:2,1), preinteresting_det(1:2,1,ii))
-          nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+            mobMask(1:2,2) = iand(negMask(1:2,2), preinteresting_det(1:2,2,ii))
-        end do
+            nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) + &
+                 popcnt(mobMask(2, 1)) + popcnt(mobMask(2, 2)) 
+          case (3)
+            mobMask(1:3,1) = iand(negMask(1:3,1), preinteresting_det(1:3,1,ii))
+            mobMask(1:3,2) = iand(negMask(1:3,2), preinteresting_det(1:3,2,ii))
+            nt = 0
+            do j=3,1,-1
+              if (mobMask(j,1) /= 0_bit_kind) then
+                nt = nt+ popcnt(mobMask(j, 1))
+                if (nt > 4) exit
+              endif
+              if (mobMask(j,2) /= 0_bit_kind) then
+                nt = nt+ popcnt(mobMask(j, 2))
+                if (nt > 4) exit
+              endif
+            end do
+          case (4)
+            mobMask(1:4,1) = iand(negMask(1:4,1), preinteresting_det(1:4,1,ii))
+            mobMask(1:4,2) = iand(negMask(1:4,2), preinteresting_det(1:4,2,ii))
+            nt = 0
+            do j=4,1,-1
+              if (mobMask(j,1) /= 0_bit_kind) then
+                nt = nt+ popcnt(mobMask(j, 1))
+                if (nt > 4) exit
+              endif
+              if (mobMask(j,2) /= 0_bit_kind) then
+                nt = nt+ popcnt(mobMask(j, 2))
+                if (nt > 4) exit
+              endif
+            end do
+          case default
+            mobMask(1:N_int,1) = iand(negMask(1:N_int,1), preinteresting_det(1:N_int,1,ii))
+            mobMask(1:N_int,2) = iand(negMask(1:N_int,2), preinteresting_det(1:N_int,2,ii))
+            nt = 0 
+            do j=N_int,1,-1
+              if (mobMask(j,1) /= 0_bit_kind) then
+                nt = nt+ popcnt(mobMask(j, 1))
+                if (nt > 4) exit
+              endif
+              if (mobMask(j,2) /= 0_bit_kind) then
+                nt = nt+ popcnt(mobMask(j, 2))
+                if (nt > 4) exit
+              endif
+            end do
+        end select
         
-         if(nt <= 4) then
+        if(nt <= 4) then
-           interesting(0) += 1
+          i = preinteresting(ii)
-           interesting(interesting(0)) = i
+          interesting(0) += 1
+          interesting(interesting(0)) = i
           minilist(1,1,interesting(0)) = preinteresting_det(1,1,ii)
           minilist(1,2,interesting(0)) = preinteresting_det(1,2,ii)
-           do j=2,N_int
+          do j=2,N_int
             minilist(j,1,interesting(0)) = preinteresting_det(j,1,ii)
             minilist(j,2,interesting(0)) = preinteresting_det(j,2,ii)
-           enddo
+          enddo
-           if(nt <= 2) then
+          if(nt <= 2) then
-             fullinteresting(0) += 1
+            fullinteresting(0) += 1
-             fullinteresting(fullinteresting(0)) = i
+            fullinteresting(fullinteresting(0)) = i
             fullminilist(1,1,fullinteresting(0)) = preinteresting_det(1,1,ii)
             fullminilist(1,2,fullinteresting(0)) = preinteresting_det(1,2,ii)
-             do j=2,N_int
+            do j=2,N_int
               fullminilist(j,1,fullinteresting(0)) = preinteresting_det(j,1,ii)
               fullminilist(j,2,fullinteresting(0)) = preinteresting_det(j,2,ii)
-             enddo
+            enddo
-           end if
+          end if
-         end if
+        end if
-
+        
       end do
       
       do ii=1,prefullinteresting(0)
@@ -481,12 +503,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
         mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
         mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i))
         nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
-        do j=2,N_int
+        if (nt > 2) cycle
+        do j=N_int,2,-1
           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 = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+          if (nt > 2) exit
         end do
-
+        
         if(nt <= 2) then
           fullinteresting(0) += 1
           fullinteresting(fullinteresting(0)) = i
@@ -722,7 +746,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   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
+  double precision, external :: get_phase_bi, mo_bielec_integral
   
   integer :: i, j, tip, ma, mi, puti, putj
   integer :: h1, h2, p1, p2, i1, i2
@@ -757,7 +781,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
         if(ma == 1) then
           mat(:, putj, puti) += coefs * hij
         else
@@ -776,7 +800,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
           if(banned(puti,putj,bant)) cycle
           p1 = p(turn2(i), 1)
           
-          hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
+          hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
           mat(:, puti, putj) += coefs * hij
         end do
       end do
@@ -796,7 +820,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
         mat(:, puti, putj) += coefs * hij
       end do
       end do
@@ -810,7 +834,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = mo_bielec_integral(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
@@ -821,7 +845,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+        hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2)
         mat(:, puti, putj) += coefs * hij
       end if
     end if
@@ -841,7 +865,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   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               :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num)
-  double precision, external     :: get_phase_bi, integral8
+  double precision, external     :: get_phase_bi, mo_bielec_integral
   logical                        :: ok
 
   logical, allocatable           :: lbanned(:,:)
@@ -881,12 +905,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       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)
+        hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(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)
+        hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(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
 
@@ -906,13 +930,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       !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)
+        hij = mo_bielec_integral(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)
+        hij = mo_bielec_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
         tmp_row2(:,puti) += hij * coefs
       end if
     end do
@@ -934,12 +958,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+          hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(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)
+          hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
           tmp_row(:,putj) += hij * coefs
         end do
 
@@ -957,13 +981,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         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)
+          hij = mo_bielec_integral(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)
+          hij = mo_bielec_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
           tmp_row2(:,puti) += hij * coefs
         end if
       end do
@@ -1015,7 +1039,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   
   integer :: i, j, s, h1, h2, p1, p2, puti, putj
   double precision :: hij, phase
-  double precision, external :: get_phase_bi, integral8
+  double precision, external :: get_phase_bi, mo_bielec_integral
   logical :: ok
   
   integer :: bant
@@ -1035,7 +1059,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
           call i_h_j(gen, det, N_int, hij)
         else
           phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
-          hij = integral8(p1, p2, h1, h2) * phase
+          hij = mo_bielec_integral(p1, p2, h1, h2) * phase
         end if
         mat(:, p1, p2) += coefs(:) * hij
       end do
@@ -1052,7 +1076,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
           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)
+          hij = (mo_bielec_integral(p1, p2, puti, putj) -  mo_bielec_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2)
         end if
         mat(:, puti, putj) += coefs(:) * hij
       end do

plugins/Full_CI_ZMQ/selection_slave.irp.f

@@ -19,13 +19,14 @@ 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)
   character*(64) :: states(4)
-  integer :: rc, i
+  integer :: rc, i, ierr
   
   call provide_everything
   

plugins/Full_CI_ZMQ_GPI2/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Full_CI_ZMQ GPI2

plugins/Full_CI_ZMQ_GPI2/README.rst

@@ -1,8 +1,9 @@
-===============
+================
-Full_CI_ZMQ_MPI
+Full_CI_ZMQ_GPI2
-===============
+================
 
-MPI Slave for Full_CI with ZMQ
+GPI2 Slave for Full_CI with ZMQ. There should be one instance of the slave
+per compute node.
 
 Needed Modules
 ==============

plugins/Full_CI_ZMQ_GPI2/selection_davidson_slave_gpi2.irp.f

@@ -14,7 +14,7 @@ end
 
 subroutine provide_everything
   PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
-  PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count MPI_Initialized
+  PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count GASPI_is_Initialized
 end
 
 subroutine run_wf
@@ -51,12 +51,10 @@ subroutine run_wf
       ! ---------
 
       print *,  'Selection'
-      if (is_mpi_master) then
+      if (is_gaspi_master) then
           call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
       endif
-      IRP_IF MIP
+      call broadcast_wf(energy)
-        call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) 
-      IRP_ENDIF
   
       !$OMP PARALLEL PRIVATE(i)
       i = omp_get_thread_num()
@@ -70,7 +68,10 @@ subroutine run_wf
       ! --------
 
       print *,  'Davidson'
-      call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
+      if (is_gaspi_master) then
+        call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
+      endif
+      call broadcast_wf(energy)
       call omp_set_nested(.True.)
       call davidson_slave_tcp(0)
       call omp_set_nested(.False.)
@@ -82,7 +83,10 @@ subroutine run_wf
       ! ---
 
       print *,  'PT2'
-      call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
+      if (is_gaspi_master) then
+        call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
+      endif
+      call broadcast_wf(energy)
   
       logical :: lstop
       lstop = .False.

plugins/Full_CI_ZMQ_MPI/NEEDED_CHILDREN_MODULES

@@ -1 +0,0 @@
-Full_CI_ZMQ MPI

plugins/GPI2/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Determinants