diff --git a/.travis.yml b/.travis.yml
index 262f1147..dd28c132 100644
--- a/.travis.yml
+++ b/.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
diff --git a/README.md b/README.md
index eebb67dd..a11c5713 100644
--- a/README.md
+++ b/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)
diff --git a/bin/qp_gaspi_run b/bin/qp_gaspi_run
new file mode 100755
index 00000000..ec54afdc
--- /dev/null
+++ b/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 $@
diff --git a/config/gfortran.cfg b/config/gfortran.cfg
index c0aa875f..2256ccf4 100644
--- a/config/gfortran.cfg
+++ b/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
 
diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg
index 91a12345..8a7ce9b8 100644
--- a/config/gfortran_debug.cfg
+++ b/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
 
diff --git a/config/travis.cfg b/config/travis.cfg
index c444db09..d2d75867 100644
--- a/config/travis.cfg
+++ b/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
 
diff --git a/configure b/configure
index c26eb7d3..cd76e907 100755
--- a/configure
+++ b/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 LIBRARY_PATH="${QP_ROOT}"/lib:"${LIBRARY_PATH}"',
+        'export LD_LIBRARY_PATH="${QP_ROOT}"/lib:"${QP_ROOT}"/lib64:"${LD_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',
diff --git a/install/scripts/install_gpi2.sh b/install/scripts/install_gpi2.sh
new file mode 100755
index 00000000..751f4ef8
--- /dev/null
+++ b/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 gpi2
+  ./install.sh -p $QP_ROOT $GPI_OPTIONS
+  cp src/GASPI.f90 $QP_ROOT/src/plugins/GPI2/
+  return 0
+}
+
+source scripts/build.sh
diff --git a/install/scripts/install_lapack.sh b/install/scripts/install_lapack.sh
new file mode 100755
index 00000000..25cbefc6
--- /dev/null
+++ b/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
+
+
+
diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh
index 80632dd2..4d356bd3 100755
--- a/install/scripts/install_ocaml.sh
+++ b/install/scripts/install_ocaml.sh
@@ -5,8 +5,7 @@ QP_ROOT=$PWD
 cd -
 
 # Normal installation
-PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ"
-#ppx_sexp_conv
+PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ ppx_sexp_conv ppx_deriving"
 
 # 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
diff --git a/ocaml/Address.ml b/ocaml/Address.ml
index c819a463..a419806c 100644
--- a/ocaml/Address.ml
+++ b/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
diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml
index 72932b1f..5d385e89 100644
--- a/ocaml/Atom.ml
+++ b/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 =
diff --git a/ocaml/Basis.ml b/ocaml/Basis.ml
index 797d53f2..30af6577 100644
--- a/ocaml/Basis.ml
+++ b/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) =
diff --git a/ocaml/Bit.ml b/ocaml/Bit.ml
index 28b8dac9..e5958ba6 100644
--- a/ocaml/Bit.ml
+++ b/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"
diff --git a/ocaml/Bit.mli b/ocaml/Bit.mli
index 6dd5abdd..bc62d13c 100644
--- a/ocaml/Bit.mli
+++ b/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
diff --git a/ocaml/Bitlist.ml b/ocaml/Bitlist.ml
index d7b9fc50..46d9ab9a 100644
--- a/ocaml/Bitlist.ml
+++ b/ocaml/Bitlist.ml
@@ -1,5 +1,5 @@
 open Qptypes
-open Core.Std
+open Core
 
 (*
 Type for bits strings
diff --git a/ocaml/Charge.ml b/ocaml/Charge.ml
index 714a5690..8c4cc429 100644
--- a/ocaml/Charge.ml
+++ b/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
diff --git a/ocaml/Charge.mli b/ocaml/Charge.mli
index 07685531..9ed8e41a 100644
--- a/ocaml/Charge.mli
+++ b/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
diff --git a/ocaml/Determinant.ml b/ocaml/Determinant.ml
index 3791e07e..a6d754b0 100644
--- a/ocaml/Determinant.ml
+++ b/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 ~f:(fun i -> x.(i+n_int)) )
+  ( Array.init n_int (fun i -> x.(i)) ,
+    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
-      | Bit.Zero -> "-"
-      | Bit.One  -> "+" 
-  )
-  |> String.concat
-  |> String.sub ~pos:0 ~len 
+  let s = 
+    List.map (function
+        | Bit.Zero -> "-"
+        | Bit.One  -> "+" 
+    ) x
+    |> 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 ""
 
 
diff --git a/ocaml/Determinant.mli b/ocaml/Determinant.mli
index da9fe02e..78a2f52f 100644
--- a/ocaml/Determinant.mli
+++ b/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
diff --git a/ocaml/Element.ml b/ocaml/Element.ml
index ebfd5e17..fd08b8da 100644
--- a/ocaml/Element.ml
+++ b/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
diff --git a/ocaml/Element.mli b/ocaml/Element.mli
index e78a57ad..fc6c679f 100644
--- a/ocaml/Element.mli
+++ b/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
diff --git a/ocaml/Excitation.ml b/ocaml/Excitation.ml
index d58c3093..58e18b11 100644
--- a/ocaml/Excitation.ml
+++ b/ocaml/Excitation.ml
@@ -1,14 +1,14 @@
-open Core.Std;;
-open Qptypes;;
+open Core
+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:" "
-;;
+
 
diff --git a/ocaml/Excitation.mli b/ocaml/Excitation.mli
index 982cfd0e..ab9e083c 100644
--- a/ocaml/Excitation.mli
+++ b/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
 
diff --git a/ocaml/GaussianPrimitive.ml b/ocaml/GaussianPrimitive.ml
index cf3d7cdb..cdaabd87 100644
--- a/ocaml/GaussianPrimitive.ml
+++ b/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
diff --git a/ocaml/Generic_input_of_rst.ml b/ocaml/Generic_input_of_rst.ml
index 81388824..02aec5ee 100644
--- a/ocaml/Generic_input_of_rst.ml
+++ b/ocaml/Generic_input_of_rst.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core;;
 open Qptypes;;
 
 
diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml
index 2c6efb88..ab265202 100644
--- a/ocaml/Gto.ml
+++ b/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 n = Int.of_string (String.strip n_str) in
+  let sym = Symmetry.of_string (String_ext.strip sym_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:' ' 
-      |> List.filter ~f:(fun x -> x <> "")
+      |> String_ext.split ~on:' ' 
+      |> 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) )
-          and c = AO_coef.of_float (Float.of_string coef) in
+              (AO_expo.of_float (float_of_string expo) )
+          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 *)
diff --git a/ocaml/Gto.mli b/ocaml/Gto.mli
index 93b6c0f3..91534ebe 100644
--- a/ocaml/Gto.mli
+++ b/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 :
diff --git a/ocaml/Input.ml b/ocaml/Input.ml
index b1bb1060..65155f7c 100644
--- a/ocaml/Input.ml
+++ b/ocaml/Input.ml
@@ -1,6 +1,6 @@
 open Qputils;;
 open Qptypes;;
-open Core.Std;;
+open Core;;
 
 include Input_ao_basis;;
 include Input_bitmasks;;
diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml
index 8b0f72a2..348ddd7f 100644
--- a/ocaml/Input_ao_basis.ml
+++ b/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";;
diff --git a/ocaml/Input_bi_integrals.ml b/ocaml/Input_bi_integrals.ml
index c5fc4fe5..1caf38db 100644
--- a/ocaml/Input_bi_integrals.ml
+++ b/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";;
diff --git a/ocaml/Input_bitmasks.ml b/ocaml/Input_bitmasks.ml
index 0469c660..75a22a26 100644
--- a/ocaml/Input_bitmasks.ml
+++ b/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";;
diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml
index 6cc83745..48887ca0 100644
--- a/ocaml/Input_determinants_by_hand.ml
+++ b/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";;
diff --git a/ocaml/Input_electrons.ml b/ocaml/Input_electrons.ml
index 24d0fe00..3779cfd2 100644
--- a/ocaml/Input_electrons.ml
+++ b/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";;
diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml
index df47abfb..c51a7f4a 100644
--- a/ocaml/Input_mo_basis.ml
+++ b/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 () = 
diff --git a/ocaml/Input_nuclei_by_hand.ml b/ocaml/Input_nuclei_by_hand.ml
index f36b6b82..b4b8b0fe 100644
--- a/ocaml/Input_nuclei_by_hand.ml
+++ b/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";;
diff --git a/ocaml/Io_ext.ml b/ocaml/Io_ext.ml
new file mode 100644
index 00000000..489ed1ed
--- /dev/null
+++ b/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"
+
diff --git a/ocaml/Long_basis.ml b/ocaml/Long_basis.ml
index 5c0ea6d6..06ea2ed5 100644
--- a/ocaml/Long_basis.ml
+++ b/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) -> 
-     "( "^((Int.to_string (Nucl_number.to_int z)))^", "^
+  let middle = List.map (fun (x,y,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
-;;
+
 
diff --git a/ocaml/Long_basis.mli b/ocaml/Long_basis.mli
index 7e69ecce..81f00539 100644
--- a/ocaml/Long_basis.mli
+++ b/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 :
diff --git a/ocaml/MO_class.ml b/ocaml/MO_class.ml
index adf1a215..e85e9dda 100644
--- a/ocaml/MO_class.ml
+++ b/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 = 
diff --git a/ocaml/MO_class.mli b/ocaml/MO_class.mli
index 953e1afe..634fa95a 100644
--- a/ocaml/MO_class.mli
+++ b/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 *)
diff --git a/ocaml/MO_label.ml b/ocaml/MO_label.ml
index 72bca28f..4e0b82f6 100644
--- a/ocaml/MO_label.ml
+++ b/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
diff --git a/ocaml/MO_label.mli b/ocaml/MO_label.mli
index d5061095..732bf1f2 100644
--- a/ocaml/MO_label.mli
+++ b/ocaml/MO_label.mli
@@ -5,7 +5,7 @@ type t =
   | Localized
   | Orthonormalized
   | None
-with sexp
+[@@deriving sexp]
 
 (** String representation *)
 val to_string : t -> string
diff --git a/ocaml/Makefile b/ocaml/Makefile
index 8519c973..3534c614 100644
--- a/ocaml/Makefile
+++ b/ocaml/Makefile
@@ -11,8 +11,8 @@ endif
 
 LIBS=
 PKGS=
-OCAMLCFLAGS="-g -warn-error A"
-OCAMLBUILD=ocamlbuild -j 0 -syntax camlp4o -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) 
+OCAMLCFLAGS="-g"
+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
diff --git a/ocaml/Message.ml b/ocaml/Message.ml
index faf5ed69..2ffd1da1 100644
--- a/ocaml/Message.ml
+++ b/ocaml/Message.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 (** New job : Request to create a new multi-tasked job *)
diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml
index eb4e0582..27ac4bd0 100644
--- a/ocaml/Molecule.ml
+++ b/ocaml/Molecule.ml
@@ -1,14 +1,14 @@
-open Core.Std ;;
-open Qptypes ;;
+open Qptypes 
+open Sexplib.Std
 
-exception MultiplicityError of string;;
-exception XYZError ;;
+exception MultiplicityError of string
+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
-        | None   -> build_list (List.Assoc.add accu e 1) rest
-        | Some i -> build_list (List.Assoc.add accu e (i+1)) rest
+        try
+          let i = List.assoc e accu in
+          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.concat ~sep:"\n"
+  [ string_of_int n ; title ] @  (List.map f nuclei)
+  |> 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'
-       |> List.filter ~f:(fun x -> x <> "")
-       |> List.map ~f:(fun x -> Atom.of_string units x)
+  let l = String_ext.split s ~on:'\n'
+       |> List.filter (fun x -> 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 
-  |> String.lsplit2_exn ~on:'\n'
+  let lines = 
+    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 =
-    try
-      (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
+  of_xyz_string ~charge:charge ~multiplicity:multiplicity 
+    ~units:units lines
 
 
 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
 
diff --git a/ocaml/Molecule.mli b/ocaml/Molecule.mli
index f6201b18..fe80b48e 100644
--- a/ocaml/Molecule.mli
+++ b/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
diff --git a/ocaml/Multiplicity.ml b/ocaml/Multiplicity.ml
index e56341c6..7570f483 100644
--- a/ocaml/Multiplicity.ml
+++ b/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 ;;
diff --git a/ocaml/Multiplicity.mli b/ocaml/Multiplicity.mli
index c6f8c6bf..0a28fb11 100644
--- a/ocaml/Multiplicity.mli
+++ b/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
diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml
index 5717ed39..7723ca82 100644
--- a/ocaml/Point3d.ml
+++ b/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 
diff --git a/ocaml/Point3d.mli b/ocaml/Point3d.mli
index 6d7428ec..476461a3 100644
--- a/ocaml/Point3d.mli
+++ b/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
diff --git a/ocaml/Primitive.mli b/ocaml/Primitive.mli
index 77cb633a..f7d8809d 100644
--- a/ocaml/Primitive.mli
+++ b/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
diff --git a/ocaml/Progress_bar.ml b/ocaml/Progress_bar.ml
index 3473ac4b..ac5ae0d7 100644
--- a/ocaml/Progress_bar.ml
+++ b/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.) }
   
 
 
diff --git a/ocaml/Pseudo.ml b/ocaml/Pseudo.ml
index 3791167d..8a59213c 100644
--- a/ocaml/Pseudo.ml
+++ b/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;
diff --git a/ocaml/Qpackage.ml b/ocaml/Qpackage.ml
index 8011b23b..e4759db5 100644
--- a/ocaml/Qpackage.ml
+++ b/ocaml/Qpackage.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core;;
 open Qptypes;;
 open Qputils;;
 
diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml
index 36dc1102..4bde831c 100644
--- a/ocaml/Qputils.ml
+++ b/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 
-      |> List.filter ~f:(fun x->
-         (String.strip x) <> "")
-      |> List.map ~f:(fun x->
-         "("^(String.tr '=' ' ' x)^")")
-      |> String.concat 
+      String_ext.split ~on:'\n' s 
+      |> List.filter (fun x-> (String_ext.strip x) <> "")
+      |> List.map (fun x-> "("^ 
+        (Str.global_replace (Str.regexp "=") " " x)
+        ^")")
+      |> 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 ->
-      match (Sys.is_directory x, Sys.is_file x) with
-      | (`Yes, _) -> remove_one x
-      | (_, `Yes) -> Sys.remove x
+    |> Array.iter (fun x ->
+      match (Sys.is_directory x, Sys.file_exists x) with
+      | (true, _) -> remove_one x
+      | (_, true) -> Sys.remove x
       | _ -> failwith ("Unable to remove file "^x^".")
     );
     Sys.chdir "..";
diff --git a/ocaml/README_qp_edit.rst b/ocaml/README_qp_edit.rst
index dfadea20..8af9cceb 100644
--- a/ocaml/README_qp_edit.rst
+++ b/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";;
diff --git a/ocaml/Range.ml b/ocaml/Range.ml
index 9607c6a3..91fcbcce 100644
--- a/ocaml/Range.ml
+++ b/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
-        and finish =  Int.of_string f
+        let start = int_of_string s
+        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 = List.map ~f:expand_range l in
-  List.concat l |> List.dedup ~compare:Int.compare |> List.sort ~cmp:Int.compare
-;;
+    let l = String_ext.split ~on:',' s in
+    let l = List.map expand_range l in
+  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 () ;
-  to_string l |> print_string ;  print_newline () ;
-;;
+  List.iter (fun x -> Printf.printf "%d, " x) l ; print_newline () ;
+  to_string l |> print_string ;  print_newline ();
+
 
diff --git a/ocaml/Range.mli b/ocaml/Range.mli
index 2d56a0fa..e186ccf9 100644
--- a/ocaml/Range.mli
+++ b/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 :
diff --git a/ocaml/String_ext.ml b/ocaml/String_ext.ml
new file mode 100644
index 00000000..ae8378bf
--- /dev/null
+++ b/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
diff --git a/ocaml/Symmetry.ml b/ocaml/Symmetry.ml
index 8647ae99..7b088b73 100644
--- a/ocaml/Symmetry.ml
+++ b/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 
-    |> List.rev 
+    String_ext.to_list s 
     |> do_work Null 
      { x=Positive_int.of_int 0 ; 
        y=Positive_int.of_int 0 ;
diff --git a/ocaml/Symmetry.mli b/ocaml/Symmetry.mli
index baffeb2e..2ab63003 100644
--- a/ocaml/Symmetry.mli
+++ b/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" *)
diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml
index 1ed403f7..45b438f2 100644
--- a/ocaml/TaskServer.ml
+++ b/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
 
diff --git a/ocaml/To_md5.ml b/ocaml/To_md5.ml
index 33015d75..bc6608f9 100644
--- a/ocaml/To_md5.ml
+++ b/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
diff --git a/ocaml/_tags b/ocaml/_tags
index 0935c0bb..aa36989b 100644
--- a/ocaml/_tags
+++ b/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
diff --git a/ocaml/create_git_sha1.sh b/ocaml/create_git_sha1.sh
index f1fb7fa6..35cbb7d5 100755
--- a/ocaml/create_git_sha1.sh
+++ b/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
-let sha1 = "$SHA1" |> String.strip
-let date = "$DATE" |> String.strip
-let message = "$MESSAGE" |> String.strip
+open Core
+let sha1 = "$SHA1" |> String_ext.strip
+let date = "$DATE" |> String_ext.strip
+let message = "$MESSAGE" |> String_ext.strip
 EOF
 
diff --git a/ocaml/qp_basis_clean.ml b/ocaml/qp_basis_clean.ml
index 32945ce3..2cc6d349 100644
--- a/ocaml/qp_basis_clean.ml
+++ b/ocaml/qp_basis_clean.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let filenames =
   let dir_name = Qpackage.root^"/data/basis/"
diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml
index 787d4a0d..410590a2 100644
--- a/ocaml/qp_create_ezfio_from_xyz.ml
+++ b/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
diff --git a/ocaml/qp_create_guess.ml b/ocaml/qp_create_guess.ml
index bebfdad3..b841c350 100644
--- a/ocaml/qp_create_guess.ml
+++ b/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
diff --git a/ocaml/qp_find_pi_space.ml b/ocaml/qp_find_pi_space.ml
new file mode 100644
index 00000000..0f5f7365
--- /dev/null
+++ b/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
+
diff --git a/ocaml/qp_overlap_of_wf.ml b/ocaml/qp_overlap_of_wf.ml
index 816256fa..deea0556 100644
--- a/ocaml/qp_overlap_of_wf.ml
+++ b/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' =
-          try  Hashtbl.find wf' k 
-          with Not_found -> 0.
+        let (c',c) =
+          try  (Hashtbl.find wf' k, c)
+          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 ()
   
diff --git a/ocaml/qp_print.ml b/ocaml/qp_print.ml
index ec584066..ea52bd7f 100644
--- a/ocaml/qp_print.ml
+++ b/ocaml/qp_print.ml
@@ -1,6 +1,6 @@
 open Qputils;;
 open Qptypes;;
-open Core.Std;;
+open Core;;
 
 type input_action =
   | Basis
diff --git a/ocaml/qp_print_basis.ml b/ocaml/qp_print_basis.ml
index 52e10ce8..477ed366 100644
--- a/ocaml/qp_print_basis.ml
+++ b/ocaml/qp_print_basis.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 let basis () =
diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml
index 5a656d2d..5be26539 100644
--- a/ocaml/qp_run.ml
+++ b/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
diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml
index 7451d87d..e82610d4 100644
--- a/ocaml/qp_set_mo_class.ml
+++ b/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)
 
 
 
diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml
index 06006181..600debf4 100644
--- a/ocaml/qptypes_generator.ml
+++ b/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 ->
-     accu && (x < 'g')));
+  assert (
+    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
-        and name = String.strip name in
+        let typ  = String_ext.strip typ
+        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
-  |> List.map ~f:(String.lsplit2_exn ~on:'\n') 
+     String_ext.split ~on:'*' input
+  |> 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'
-      |> List.filter ~f:(fun x -> (String.strip x) <> "")
+      String_ext.split s ~on:'\n'
+      |> 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
-  |> List.map ~f:parse
-  |> String.concat 
+     String_ext.split ~on:'*' input
+  |> List.map parse
+  |> String.concat ""
+  |> global_replace
   |> print_string
-;;
+
 
 let () = 
   parse_input input_data ;
   parse_input_ezfio input_ezfio;
-  print_endline untouched;
+  print_endline untouched
 
 
 
diff --git a/ocaml/test_basis.ml b/ocaml/test_basis.ml
index f58d30db..9007798c 100644
--- a/ocaml/test_basis.ml
+++ b/ocaml/test_basis.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core;;
 open Qputils;;
 open Qptypes;;
 
diff --git a/ocaml/test_gto.ml b/ocaml/test_gto.ml
index 423df62b..0ad6df56 100644
--- a/ocaml/test_gto.ml
+++ b/ocaml/test_gto.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 let test_prim () = 
diff --git a/ocaml/test_message.ml b/ocaml/test_message.ml
index 2f5592ec..6ccc381e 100644
--- a/ocaml/test_message.ml
+++ b/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"
diff --git a/ocaml/test_molecule.ml b/ocaml/test_molecule.ml
index 5d2935ed..e4b0f616 100644
--- a/ocaml/test_molecule.ml
+++ b/ocaml/test_molecule.ml
@@ -1,4 +1,4 @@
-open Core.Std ;;
+open Core ;;
 open Qptypes ;;
 
 let test_molecule () =
diff --git a/ocaml/test_pseudo.ml b/ocaml/test_pseudo.ml
index 10c8aa0c..3b62680d 100644
--- a/ocaml/test_pseudo.ml
+++ b/ocaml/test_pseudo.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qputils
 open Qptypes
 
diff --git a/ocaml/test_queuing_system.ml b/ocaml/test_queuing_system.ml
index aa2fa280..728528b2 100644
--- a/ocaml/test_queuing_system.ml
+++ b/ocaml/test_queuing_system.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let () =
 
diff --git a/ocaml/test_symmetry.ml b/ocaml/test_symmetry.ml
index bdc28e5c..e8323c14 100644
--- a/ocaml/test_symmetry.ml
+++ b/ocaml/test_symmetry.ml
@@ -1,15 +1,13 @@
-open Core.Std
 open Qputils
 open Qptypes
 open Symmetry
 
 let () =
   "SPDFGHIJKL"
-  |> String.to_list_rev
-  |> List.rev
-  |> List.map ~f:of_char
-  |> List.map ~f:Xyz.of_symmetry
-  |> List.iter ~f:(fun x -> List.iter x ~f:(fun y -> Xyz.to_string y |> print_endline) ;
+  |> String_ext.to_list
+  |> List.map of_char
+  |> List.map Xyz.of_symmetry
+  |> List.iter (fun x -> List.iter (fun y -> Xyz.to_string y |> print_endline) x ;
      print_newline ();)
 
 
diff --git a/ocaml/test_task_server.ml b/ocaml/test_task_server.ml
index e6a6106b..00573a9d 100644
--- a/ocaml/test_task_server.ml
+++ b/ocaml/test_task_server.ml
@@ -1,5 +1,3 @@
-open Core
-
 let () =
   TaskServer.run 12345
 
diff --git a/plugins/CAS_SD_ZMQ/selection.irp.f b/plugins/CAS_SD_ZMQ/selection.irp.f
index eb782cb3..b42116bd 100644
--- a/plugins/CAS_SD_ZMQ/selection.irp.f
+++ b/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
diff --git a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f
index 40d75fc4..7c30e175 100644
--- a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f
+++ b/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
 
 
diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f
index 8e2f2e53..6d359696 100644
--- a/plugins/FOBOCI/new_approach.irp.f
+++ b/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
  
diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f
index 7d194a54..05f07b22 100644
--- a/plugins/FOBOCI/routines_foboci.irp.f
+++ b/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)
 
diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f
index acf19392..f404d069 100644
--- a/plugins/Full_CI_ZMQ/selection.irp.f
+++ b/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
@@ -722,7 +699,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 +734,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 +753,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 +773,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 +787,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 +798,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 +818,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 +858,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 +883,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 +911,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 +934,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 +992,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 +1012,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 +1029,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
diff --git a/plugins/Full_CI_ZMQ/selection_slave.irp.f b/plugins/Full_CI_ZMQ/selection_slave.irp.f
index ba85ca82..ebfcbaf2 100644
--- a/plugins/Full_CI_ZMQ/selection_slave.irp.f
+++ b/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
   
diff --git a/plugins/Full_CI_ZMQ_GPI2/NEEDED_CHILDREN_MODULES b/plugins/Full_CI_ZMQ_GPI2/NEEDED_CHILDREN_MODULES
new file mode 100644
index 00000000..dd79ddf2
--- /dev/null
+++ b/plugins/Full_CI_ZMQ_GPI2/NEEDED_CHILDREN_MODULES
@@ -0,0 +1 @@
+Full_CI_ZMQ GPI2
diff --git a/plugins/Full_CI_ZMQ_MPI/README.rst b/plugins/Full_CI_ZMQ_GPI2/README.rst
similarity index 61%
rename from plugins/Full_CI_ZMQ_MPI/README.rst
rename to plugins/Full_CI_ZMQ_GPI2/README.rst
index bf9dfab4..b3e89d9e 100644
--- a/plugins/Full_CI_ZMQ_MPI/README.rst
+++ b/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
 ==============
diff --git a/plugins/Full_CI_ZMQ_MPI/selection_davidson_slave_mpi.irp.f b/plugins/Full_CI_ZMQ_GPI2/selection_davidson_slave_gpi2.irp.f
similarity index 82%
rename from plugins/Full_CI_ZMQ_MPI/selection_davidson_slave_mpi.irp.f
rename to plugins/Full_CI_ZMQ_GPI2/selection_davidson_slave_gpi2.irp.f
index 6e4bf775..52822e41 100644
--- a/plugins/Full_CI_ZMQ_MPI/selection_davidson_slave_mpi.irp.f
+++ b/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 MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) 
-      IRP_ENDIF
+      call broadcast_wf(energy)
   
       !$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.
diff --git a/plugins/Full_CI_ZMQ_MPI/NEEDED_CHILDREN_MODULES b/plugins/Full_CI_ZMQ_MPI/NEEDED_CHILDREN_MODULES
deleted file mode 100644
index 03d04444..00000000
--- a/plugins/Full_CI_ZMQ_MPI/NEEDED_CHILDREN_MODULES
+++ /dev/null
@@ -1 +0,0 @@
-Full_CI_ZMQ MPI
diff --git a/plugins/GPI2/NEEDED_CHILDREN_MODULES b/plugins/GPI2/NEEDED_CHILDREN_MODULES
new file mode 100644
index 00000000..aae89501
--- /dev/null
+++ b/plugins/GPI2/NEEDED_CHILDREN_MODULES
@@ -0,0 +1 @@
+Determinants
diff --git a/plugins/MPI/README.rst b/plugins/GPI2/README.rst
similarity index 76%
rename from plugins/MPI/README.rst
rename to plugins/GPI2/README.rst
index 7962296b..d6be4958 100644
--- a/plugins/MPI/README.rst
+++ b/plugins/GPI2/README.rst
@@ -1,8 +1,8 @@
-===
-MPI
-===
+=====
+GASPI
+=====
 
-Providers for MPI programs.
+Providers for GASPI programs (with the GPI2 library).
 
 Needed Modules
 ==============
diff --git a/plugins/GPI2/gpi_test.irp.f b/plugins/GPI2/gpi_test.irp.f
new file mode 100644
index 00000000..1fc109e6
--- /dev/null
+++ b/plugins/GPI2/gpi_test.irp.f
@@ -0,0 +1,13 @@
+program test
+  double precision :: energy(N_states)
+  if (is_gaspi_master) then
+    energy = 1.d0
+  else
+    energy = 0.d0
+  endif
+  call broadcast_wf(energy)
+  print *,  'energy (1.d0) :', GASPI_rank, energy(1)
+  print *,  'coef          :', GASPI_rank, psi_coef(1,1)
+  print *,  'det           :', GASPI_rank, psi_det (1,1,1)
+  call gaspi_finalize
+end
diff --git a/plugins/GPI2/utils.irp.f b/plugins/GPI2/utils.irp.f
new file mode 100644
index 00000000..cfb17b75
--- /dev/null
+++ b/plugins/GPI2/utils.irp.f
@@ -0,0 +1,76 @@
+ BEGIN_PROVIDER [ logical, GASPI_is_initialized ]
+&BEGIN_PROVIDER [ logical, has_gaspi ]
+  implicit none
+  BEGIN_DOC
+! This is true when GASPI_Init has been called
+  END_DOC
+
+  has_gaspi = .False.
+  IRP_IF GASPI
+      use GASPI
+      integer(gaspi_return_t) :: res
+      res = gaspi_proc_init(GASPI_BLOCK)
+      if (res /= GASPI_SUCCESS) then
+         print *, res
+         print *, 'GASPI failed to initialize'
+         stop -1
+      endif
+      has_gaspi = .True.
+  IRP_ENDIF
+  GASPI_is_initialized = .True.
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, GASPI_rank ]
+&BEGIN_PROVIDER [ integer, GASPI_size ]
+&BEGIN_PROVIDER [ logical, is_GASPI_master ]
+ implicit none
+ BEGIN_DOC
+! Usual GASPI variables
+ END_DOC
+
+ PROVIDE GASPI_is_initialized
+
+ IRP_IF GASPI
+   use GASPI
+   integer(gaspi_return_t) :: res 
+   integer(gaspi_rank_t)   :: n
+   res = gaspi_proc_num(n)
+   GASPI_size = n
+   if (res  /= GASPI_SUCCESS) then
+     print *, res 
+     print *, 'Unable to get GASPI_size'
+     stop -1
+   endif
+   res = gaspi_proc_rank(n)
+   GASPI_rank = n
+   if (res  /= GASPI_SUCCESS) then
+     print *, res 
+     print *, 'Unable to get GASPI_rank'
+     stop -1
+   endif
+   is_GASPI_master = (GASPI_rank == 0)
+ IRP_ELSE
+   GASPI_rank = 0
+   GASPI_size = 1
+   is_GASPI_master = .True.
+ IRP_ENDIF
+
+
+END_PROVIDER
+
+subroutine gaspi_finalize()
+  implicit none
+  PROVIDE GASPI_is_initialized
+  IRP_IF GASPI
+    use GASPI
+    integer(gaspi_return_t) :: res 
+    res = gaspi_proc_term(GASPI_BLOCK)
+    if (res  /= GASPI_SUCCESS) then
+      print *, res 
+      print *, 'Unable to finalize GASPI'
+      stop -1
+    endif
+  IRP_ENDIF
+end subroutine
+
diff --git a/plugins/Hartree_Fock/DIIS.irp.f b/plugins/Hartree_Fock/DIIS.irp.f
index 614a9173..6ebb5879 100644
--- a/plugins/Hartree_Fock/DIIS.irp.f
+++ b/plugins/Hartree_Fock/DIIS.irp.f
@@ -94,7 +94,7 @@ END_PROVIDER
 
   do i=1,AO_num
     do j=1,AO_num
-      Xt(i,j) = X_Matrix_AO(j,i)
+      Xt(i,j) = S_half_inv(j,i)
     enddo
   enddo
 
@@ -103,7 +103,7 @@ END_PROVIDER
   call dgemm('N','N',AO_num,AO_num,AO_num,     &
        1.d0,                                   &
        Fock_matrix_AO,size(Fock_matrix_AO,1),  &
-       X_Matrix_AO,size(X_Matrix_AO,1),        &
+       S_half_inv,size(S_half_inv,1),        &
        0.d0,                                   &
        eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))       
 
@@ -130,67 +130,10 @@ END_PROVIDER
 
   call dgemm('N','N',AO_num,AO_num,AO_num,     &
        1.d0,                                   &
-       X_matrix_AO,size(X_matrix_AO,1),        &
+       S_half_inv,size(S_half_inv,1),        &
        scratch,size(scratch,1),                &
        0.d0,                                   &
        eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))       
    
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num,AO_num) ]
-
-  BEGIN_DOC
-!   Matrix X = S^{-1/2} obtained by SVD
-  END_DOC
-
-  implicit none
-  
-  integer                         :: num_linear_dependencies
-  integer                         :: LDA, LDC
-  double precision, allocatable   :: U(:,:),Vt(:,:), D(:)
-  integer                         :: info, i, j, k
- 
-  LDA = size(AO_overlap,1)
-  LDC = size(X_matrix_AO,1)
- 
-  allocate(         &
-    U(LDC,AO_num),  &
-    Vt(LDA,AO_num), &
-    D(AO_num))
-
-  call svd(            &
-       AO_overlap,LDA, &
-       U,LDC,          &
-       D,              &
-       Vt,LDA,         &
-       AO_num,AO_num)
-
-  num_linear_dependencies = 0
-  do i=1,AO_num
-    print*,D(i)
-    if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
-      D(i) = 0.d0
-      num_linear_dependencies += 1
-    else
-      ASSERT (D(i) > 0.d0)
-      D(i) = 1.d0/sqrt(D(i))
-    endif
-    do j=1,AO_num
-      X_matrix_AO(j,i) = 0.d0
-    enddo
-  enddo
-  write(*,*) 'linear dependencies',num_linear_dependencies
-!  stop
-
-  do k=1,AO_num
-    if(D(k) /= 0.d0) then
-      do j=1,AO_num
-        do i=1,AO_num
-          X_matrix_AO(i,j) = X_matrix_AO(i,j) + U(i,k)*D(k)*Vt(k,j)
-        enddo
-      enddo
-    endif
-  enddo
-  
-
-END_PROVIDER
diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg
index a4b646e1..7f67d961 100644
--- a/plugins/Hartree_Fock/EZFIO.cfg
+++ b/plugins/Hartree_Fock/EZFIO.cfg
@@ -1,9 +1,3 @@
-[threshold_overlap_ao_eigenvalues]
-type: Threshold
-doc: Threshold on the magnitude of the smallest eigenvalues of the overlap matrix in the AO basis
-interface: ezfio,provider,ocaml
-default: 1.e-6
-
 [max_dim_diis]
 type: integer
 doc: Maximum size of the DIIS extrapolation procedure
@@ -32,7 +26,7 @@ default: 500
 type: Positive_float
 doc: Energy shift on the virtual MOs to improve SCF convergence
 interface: ezfio,provider,ocaml
-default: 0.0
+default: 0.1
 
 [scf_algorithm]
 type: character*(32)
diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f
index 087933c8..0764c83f 100644
--- a/plugins/Hartree_Fock/Fock_matrix.irp.f
+++ b/plugins/Hartree_Fock/Fock_matrix.irp.f
@@ -263,17 +263,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num,mo_tot_num)
    BEGIN_DOC
    ! Fock matrix on the MO basis
    END_DOC
-   double precision, allocatable  :: T(:,:)
-   allocate ( T(ao_num,mo_tot_num) )
-   call dgemm('N','N', ao_num, mo_tot_num, ao_num,                   &
-       1.d0, Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1),      &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0, T, size(T,1))
-   call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num,               &
-       1.d0, mo_coef,size(mo_coef,1),                                &
-       T, size(T,1),                                                 &
-       0.d0, Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1))
-   deallocate(T)
+   call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), &
+                 Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1))
 END_PROVIDER
  
  
@@ -282,17 +273,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num,mo_tot_num)
    BEGIN_DOC
    ! Fock matrix on the MO basis
    END_DOC
-   double precision, allocatable  :: T(:,:)
-   allocate ( T(ao_num,mo_tot_num) )
-   call dgemm('N','N', ao_num, mo_tot_num, ao_num,                   &
-       1.d0, Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1),        &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0, T, size(T,1))
-   call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num,               &
-       1.d0, mo_coef,size(mo_coef,1),                                &
-       T, size(T,1),                                                 &
-       0.d0, Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1))
-   deallocate(T)
+   call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), &
+                 Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1))
 END_PROVIDER
  
 BEGIN_PROVIDER [ double precision, HF_energy ]
@@ -330,97 +312,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ]
      enddo
    enddo
  else
-   double precision, allocatable  :: T(:,:), M(:,:)
-   integer                        :: ierr
-   ! F_ao = S C F_mo C^t S
-   allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr)
-   if (ierr /=0 ) then
-      print *,  irp_here, ' : allocation failed'
-   endif
-
-!  ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
-!  -> M(ao_num,mo_tot_num)
-   call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0,               &
-       ao_overlap, size(ao_overlap,1),                               &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0,                                                         &
-       M, size(M,1))
-
-!  M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num)
-!  -> T(ao_num,mo_tot_num)
-   call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,           &
-       M, size(M,1),                                                 &
-       Fock_matrix_mo, size(Fock_matrix_mo,1),                       &
-       0.d0,                                                         &
-       T, size(T,1))
-
-!  T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
-!  -> M(ao_num,ao_num)
-   call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0,           &
-       T, size(T,1),                                                 &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0,                                                         &
-       M, size(M,1))
-
-!  M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
-!  -> Fock_matrix_ao(ao_num,ao_num)
-   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
-       M, size(M,1),                                                 &
-       ao_overlap, size(ao_overlap,1),                               &
-       0.d0,                                                         &
-       Fock_matrix_ao, size(Fock_matrix_ao,1))
-
-
-   deallocate(T)
+   call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), &
+      Fock_matrix_ao,size(Fock_matrix_ao,1)) 
  endif
 END_PROVIDER
 
-subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
-  implicit none
-  integer, intent(in)            :: LDFMO ! size(FMO,1)
-  integer, intent(in)            :: LDFAO ! size(FAO,1)
-  double precision, intent(in)   :: FMO(LDFMO,*)
-  double precision, intent(out)  :: FAO(LDFAO,*)
-  
-  double precision, allocatable  :: T(:,:), M(:,:)
-  integer                        :: ierr
-  ! F_ao = S C F_mo C^t S
-  allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr)
-  if (ierr /=0 ) then
-     print *,  irp_here, ' : allocation failed'
-  endif
-
-!  ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
-!  -> M(ao_num,mo_tot_num)
-  call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0,                    &
-      ao_overlap, size(ao_overlap,1),                                &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0,                                                          &
-      M, size(M,1))
-
-!  M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num)
-!  -> T(ao_num,mo_tot_num)
-  call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
-      M, size(M,1),                                                  &
-      FMO, size(FMO,1),                                              &
-      0.d0,                                                          &
-      T, size(T,1))
-
-!  T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
-!  -> M(ao_num,ao_num)
-  call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0,            &
-      T, size(T,1),                                                  &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0,                                                          &
-      M, size(M,1))
-
-!  M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
-!  -> Fock_matrix_ao(ao_num,ao_num)
-  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
-      M, size(M,1),                                                  &
-      ao_overlap, size(ao_overlap,1),                                &
-      0.d0,                                                          &
-      FAO, size(FAO,1))
-  deallocate(T,M)
-end
 
diff --git a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
index 9e0d29e5..736d1a97 100644
--- a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
+++ b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
@@ -1,7 +1,7 @@
 BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! S^-1 x Alpha density matrix in the AO basis x S^-1
+   ! S^{-1}.P_alpha.S^{-1}
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
@@ -14,7 +14,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! S^-1 Beta density matrix in the AO basis x S^-1
+   ! S^{-1}.P_beta.S^{-1}
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
@@ -27,7 +27,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! S^-1 Density matrix in the AO basis S^-1
+   ! S^{-1}.P.S^{-1}  where P = C.C^t
    END_DOC
    ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
    if (elec_alpha_num== elec_beta_num) then
diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
index 5d8097d9..241721ae 100644
--- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
+++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
@@ -133,7 +133,7 @@ END_DOC
   write(output_hartree_fock,*)
   
   if(.not.no_oa_or_av_opt)then
-   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
+   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
   endif
 
   call write_double(output_hartree_fock, Energy_SCF, 'Hartree-Fock energy')
diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f
index 5336e8e0..3d71d826 100644
--- a/plugins/Hartree_Fock/SCF.irp.f
+++ b/plugins/Hartree_Fock/SCF.irp.f
@@ -23,7 +23,7 @@ subroutine create_guess
       mo_coef = ao_ortho_lowdin_coef
       TOUCH mo_coef
       mo_label = 'Guess'
-      call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label)
+      call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.)
       SOFT_TOUCH mo_coef mo_label
     else if (mo_guess_type == "Huckel") then
       call huckel_guess
diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f
index 2be14ed3..20a8abd7 100644
--- a/plugins/Hartree_Fock/damping_SCF.irp.f
+++ b/plugins/Hartree_Fock/damping_SCF.irp.f
@@ -119,7 +119,7 @@ subroutine damping_SCF
   write(output_hartree_fock,*)
   
   if(.not.no_oa_or_av_opt)then
-   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
+   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
   endif
 
   call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')
diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f
index 9cb84507..c9e32ad5 100644
--- a/plugins/Hartree_Fock/huckel.irp.f
+++ b/plugins/Hartree_Fock/huckel.irp.f
@@ -7,25 +7,28 @@ subroutine huckel_guess
   double precision               :: accu
   double precision               :: c
   character*(64)                 :: label
-
+  double precision, allocatable  :: A(:,:)
   label = "Guess"
-  call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
-                                     size(mo_mono_elec_integral,1),  &
-                                     size(mo_mono_elec_integral,2),label,1)
-  TOUCH mo_coef
-
   c = 0.5d0 * 1.75d0
 
+  allocate (A(ao_num, ao_num))
+  A = 0.d0
   do j=1,ao_num
     do i=1,ao_num
-      Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + &
-                                                 ao_mono_elec_integral_diag(j))
+      A(i,j) = c * ao_overlap(i,j) * (ao_mono_elec_integral_diag(i) + ao_mono_elec_integral_diag(j))
     enddo
-    Fock_matrix_ao(j,j) = Fock_matrix_ao_alpha(j,j)
+    A(j,j) = ao_mono_elec_integral_diag(j) + ao_bi_elec_integral_alpha(j,j)
   enddo
-  TOUCH Fock_matrix_ao
+
+  Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
+  Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
+
+!  TOUCH mo_coef
+
+  TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta
   mo_coef = eigenvectors_fock_matrix_mo
   SOFT_TOUCH mo_coef
   call save_mos
+  deallocate(A)
 
 end
diff --git a/plugins/MPI/NEEDED_CHILDREN_MODULES b/plugins/MPI/NEEDED_CHILDREN_MODULES
deleted file mode 100644
index 8b137891..00000000
--- a/plugins/MPI/NEEDED_CHILDREN_MODULES
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/plugins/MPI/bcast.irp.f b/plugins/MPI/bcast.irp.f
deleted file mode 100644
index 3236f9f5..00000000
--- a/plugins/MPI/bcast.irp.f
+++ /dev/null
@@ -1,46 +0,0 @@
-subroutine mpi_bcast_psi()
-  use f77_zmq
-  implicit none
-  BEGIN_DOC
-! Put the wave function on the qp_run scheduler
-  END_DOC
-  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
-  integer, intent(in)            :: worker_id
-  integer                        :: ierr
-  character*(256)                :: msg
-
-  IRP_IF MPI
-    call MPI_BCast(N_states,     1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    call MPI_BCast(N_det,        1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    call MPI_BCast(psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-  
-    TOUCH psi_det_size N_det N_states
-
-    call MPI_BCast(psi_det,         N_det, MPI_INTEGER8, 0, MPI_COMM_WORLD, ierr)
-    call MPI_BCast(psi_coef, psi_det_size, MPI_DOUBLE_PRECISION* N_states, 0, MPI_COMM_WORLD, ierr)
-
-
-  rc8 = f77_zmq_recv8(zmq_to_qp_run_socket,psi_coef,int(psi_det_size*N_states*8_8,8),0)
-  if (rc8 /= psi_det_size*N_states*8_8) then
-    print *, '77_zmq_recv8(zmq_to_qp_run_socket,psi_coef,psi_det_size*N_states*8,ZMQ_SNDMORE)'
-    stop 'error'
-  endif
-  TOUCH psi_det psi_coef
-
-  rc = f77_zmq_recv(zmq_to_qp_run_socket,energy,size_energy*8,0)
-  if (rc /= size_energy*8) then
-    print *, '77_zmq_recv(zmq_to_qp_run_socket,energy,size_energy*8,0)'
-    stop 'error'
-  endif
-
-  if (N_det_generators_read > 0) then
-    N_det_generators = N_det_generators_read
-    TOUCH N_det_generators
-  endif
-  if (N_det_selectors_read > 0) then
-    N_det_selectors  = N_det_selectors_read
-    TOUCH N_det_selectors
-  endif
-
-end
-
diff --git a/plugins/MPI/utils.irp.f b/plugins/MPI/utils.irp.f
deleted file mode 100644
index aa4e0742..00000000
--- a/plugins/MPI/utils.irp.f
+++ /dev/null
@@ -1,68 +0,0 @@
- BEGIN_PROVIDER [ logical, MPI_Initialized ]
-&BEGIN_PROVIDER [ logical, has_mpi ]
-  implicit none
-  BEGIN_DOC
-! This is true when MPI_Init has been called
-  END_DOC
-
-  IRP_IF MPI
-      integer :: ierr
-      call MPI_Init(ierr)
-      if (ierr /= 0) then
-         print *, ierr
-         print *, 'MPI failed to initialize'
-         stop -1
-      endif
-  IRP_ENDIF
-  MPI_Initialized = .True.
-END_PROVIDER
-
-
- BEGIN_PROVIDER [ integer, MPI_rank ]
-&BEGIN_PROVIDER [ integer, MPI_size ]
-&BEGIN_PROVIDER [ logical, is_MPI_master ]
- implicit none
- BEGIN_DOC
-! Usual MPI variables
- END_DOC
-
- PROVIDE MPI_Initialized
-
- IRP_IF MPI
-   integer :: ierr
-   call mpi_comm_size(MPI_COMM_WORLD, MPI_size, ierr)
-   if (ierr /= 0) then
-     print *, ierr
-     print *, 'Unable to get MPI_size'
-     stop -1
-   endif
-   call mpi_comm_rank(MPI_COMM_WORLD, MPI_rank, ierr)
-   if (ierr /= 0) then
-     print *, ierr
-     print *, 'Unable to get MPI_rank'
-     stop -1
-   endif
-   is_MPI_master = (MPI_rank == 0)
- IRP_ELSE
-   MPI_rank = 0
-   MPI_size = 1
-   is_MPI_master = .True.
- IRP_ENDIF
-
-
-END_PROVIDER
-
-subroutine qp_mpi_finalize()
-  implicit none
-  PROVIDE MPI_Initialized
-  IRP_IF MPI
-      integer :: ierr
-      call MPI_Finalize(ierr)
-      if (ierr /= 0) then
-         print *, ierr
-         print *, 'Unable to finalize MPI'
-         stop -1
-       endif
-  IRP_ENDIF
-end subroutine
-
diff --git a/plugins/MRPT/jm_mrpt2.irp.f b/plugins/MRPT/jm_mrpt2.irp.f
index 89deef2e..df9a1dbc 100644
--- a/plugins/MRPT/jm_mrpt2.irp.f
+++ b/plugins/MRPT/jm_mrpt2.irp.f
@@ -11,7 +11,7 @@ end
 subroutine routine_3
  implicit none
 !provide fock_virt_total_spin_trace
- provide delta_ij 
+ provide delta_ij_mrpt 
  
  print *,  'N_det    = ', N_det
  print *,  'N_states = ', N_states
diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f
index 83f087bc..34d26127 100644
--- a/plugins/MRPT_Utils/mrpt_utils.irp.f
+++ b/plugins/MRPT_Utils/mrpt_utils.irp.f
@@ -1,5 +1,5 @@
 
- BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
+ BEGIN_PROVIDER [ double precision, delta_ij_mrpt, (N_det,N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, second_order_pt_new, (N_states) ]
 &BEGIN_PROVIDER [ double precision, second_order_pt_new_1h, (N_states) ]
 &BEGIN_PROVIDER [ double precision, second_order_pt_new_1p, (N_states) ]
@@ -19,7 +19,7 @@
  double precision, allocatable :: delta_ij_tmp(:,:,:)
 
 
- delta_ij = 0.d0
+ delta_ij_mrpt = 0.d0
 
  allocate (delta_ij_tmp(N_det,N_det,N_states))
 
@@ -32,7 +32,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
    second_order_pt_new_1h(i_state) = accu(i_state) 
@@ -47,7 +47,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_1p(i_state) = accu(i_state) 
@@ -63,7 +63,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_1h1p(i_state) = accu(i_state) 
@@ -79,7 +79,7 @@
   do i = 1, N_det
    do j = 1, N_det
     accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-    delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+    delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
    enddo
   enddo
   second_order_pt_new_1h1p(i_state) = accu(i_state) 
@@ -95,7 +95,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_2h(i_state) = accu(i_state) 
@@ -110,7 +110,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_2p(i_state) = accu(i_state) 
@@ -126,7 +126,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_1h2p(i_state) = accu(i_state) 
@@ -142,7 +142,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_2h1p(i_state) = accu(i_state) 
@@ -157,7 +157,7 @@
 !do i = 1, N_det
 ! do j = 1, N_det
 !  accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-!  delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+!  delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
 ! enddo
 !enddo
 !second_order_pt_new_2h2p(i_state) = accu(i_state) 
@@ -168,7 +168,7 @@
  call give_2h2p(contrib_2h2p)
  do i_state = 1, N_states
  do i = 1, N_det
-   delta_ij(i,i,i_state) += contrib_2h2p(i_state)
+   delta_ij_mrpt(i,i,i_state) += contrib_2h2p(i_state)
  enddo
  second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) 
  enddo
@@ -179,9 +179,9 @@
  accu = 0.d0
  do i_state = 1, N_states
  do i = 1, N_det
-! write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
+! write(*,'(1000(F16.10,x))')delta_ij_mrpt(i,:,:)
   do j = i_state, N_det
-   accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
+   accu(i_state) += delta_ij_mrpt(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
   enddo
  enddo
  second_order_pt_new(i_state) = accu(i_state) 
@@ -199,7 +199,7 @@ END_PROVIDER
  do i_state = 1, N_states
   do i = 1,N_det
    do j = 1,N_det
-    Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij(j,i,i_state)
+    Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij_mrpt(j,i,i_state)
    enddo
   enddo
  enddo
@@ -214,7 +214,7 @@ END_PROVIDER
   do i = 1,N_det
    do j = i,N_det
     Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) =  H_matrix_all_dets(j,i) & 
-                                            + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) )
+                                            + 0.5d0 * ( delta_ij_mrpt(j,i,i_state) + delta_ij_mrpt(i,j,i_state) )
     Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) =  Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) 
    enddo
   enddo
diff --git a/plugins/QMC/qp_convert_qmcpack_to_ezfio.py b/plugins/QMC/qp_convert_qmcpack_to_ezfio.py
index 5d63f191..8e8cea7b 100755
--- a/plugins/QMC/qp_convert_qmcpack_to_ezfio.py
+++ b/plugins/QMC/qp_convert_qmcpack_to_ezfio.py
@@ -258,7 +258,7 @@ def print_mo_coef(mo_coef_block, l_l_sym):
             i_a = int(l[1]) - 1
             sym = l[2]
 
-            print l_label[i_a], sym, " ".join('{0: 3.8f}'.format(i)
+            print l_label[i_a], sym, " ".join('%20.15e'%i
                                               for i in a[i])
 
         if i_block != nb_block - 1:
diff --git a/plugins/analyze_wf/attachment.irp.f b/plugins/analyze_wf/attachment.irp.f
new file mode 100644
index 00000000..d086aa21
--- /dev/null
+++ b/plugins/analyze_wf/attachment.irp.f
@@ -0,0 +1,64 @@
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_detachment, (mo_tot_num,mo_tot_num,2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_attachment, (mo_tot_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Detachment and attachment density matrices in MO basis
+ END_DOC
+ integer :: i,j, k, istate
+ double precision :: km(mo_tot_num), kp(mo_tot_num)
+
+ one_body_dm_mo_detachment = 0.d0
+ one_body_dm_mo_attachment = 0.d0
+
+ do istate=2,N_states
+
+   km(:) = 0.d0
+   kp(:) = 0.d0
+   do i=1,mo_tot_num
+      if (one_body_dm_mo_diff_eigvalues(i,istate) < 0) then
+        km(i) = -one_body_dm_mo_diff_eigvalues(i,istate)
+      else
+        kp(i) =  one_body_dm_mo_diff_eigvalues(i,istate)
+      endif
+   enddo
+
+   ! Attachment
+  do k=1,mo_tot_num
+     do j=1,mo_tot_num
+       do i=1,mo_tot_num
+         one_body_dm_mo_detachment(i,j,istate) = one_body_dm_mo_detachment(i,j,istate) + &
+           one_body_dm_mo_diff_eigvectors(i,k,istate)*km(k)* &
+           one_body_dm_mo_diff_eigvectors(j,k,istate)
+         one_body_dm_mo_attachment(i,j,istate) = one_body_dm_mo_attachment(i,j,istate) + &
+           one_body_dm_mo_diff_eigvectors(i,k,istate)*kp(k)* &
+           one_body_dm_mo_diff_eigvectors(j,k,istate)
+       enddo
+    enddo
+  enddo
+
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, one_body_dm_ao_detachment, (ao_num,ao_num,2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_ao_attachment, (ao_num,ao_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Detachment and attachment density matrices in AO basis
+ END_DOC
+ integer :: istate
+ do istate=2,N_states
+   call mo_to_ao_no_overlap(                                         &
+       one_body_dm_mo_attachment(1,1,istate),                        &
+       size(one_body_dm_mo_attachment,1),                            &
+       one_body_dm_ao_attachment(1,1,istate),                        &
+       size(one_body_dm_ao_attachment,1) )
+   call mo_to_ao_no_overlap(                                         &
+       one_body_dm_mo_detachment(1,1,istate),                        &
+       size(one_body_dm_mo_detachment,1),                            &
+       one_body_dm_ao_detachment(1,1,istate),                        &
+       size(one_body_dm_ao_detachment,1) )
+ enddo
+
+END_PROVIDER
+
diff --git a/plugins/analyze_wf/dump_nto.irp.f b/plugins/analyze_wf/dump_nto.irp.f
new file mode 100644
index 00000000..8d19c3eb
--- /dev/null
+++ b/plugins/analyze_wf/dump_nto.irp.f
@@ -0,0 +1,19 @@
+program dump_nto
+  implicit none
+  integer :: i,j, istate
+
+  print *,  'Phi_S'
+  do i=2,N_states
+    print *,  i, Phi_S(i)
+  enddo
+
+  do istate=2,N_states
+    do j=1,mo_tot_num
+      print *,  'MO: ', j, 'State:', istate, 'Eig:', one_body_dm_mo_diff_eigvalues(j,istate)
+      do i=1,ao_num
+        print *,  i, transition_natorb(i,j,istate)
+      enddo
+    enddo
+  enddo
+
+end
diff --git a/plugins/analyze_wf/dump_one_body_mos.irp.f b/plugins/analyze_wf/dump_one_body_mos.irp.f
new file mode 100644
index 00000000..7ab841ef
--- /dev/null
+++ b/plugins/analyze_wf/dump_one_body_mos.irp.f
@@ -0,0 +1,38 @@
+program dump_one_body_mos
+  implicit none
+  BEGIN_DOC
+! Output density matrices of all the states
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run()
+end
+
+subroutine run
+  implicit none
+  integer                        :: istate
+  integer, parameter             :: iunit = 66
+  character*(64)                 :: filename, fmt
+  integer                        :: i,j,k
+
+  write(fmt,'(''('',I4.4,''(X,E20.14))'')') mo_tot_num
+  do istate=1,N_states
+    write(filename,'(''state.'',I2.2)') istate
+    open(unit=iunit, form='formatted', file=filename)
+    write(iunit,*) mo_tot_num
+    do j=1,mo_tot_num
+      write(iunit,fmt) one_body_dm_mo_alpha(1:mo_tot_num,j,istate) + one_body_dm_mo_beta(1:mo_tot_num,j,istate)
+    enddo
+  enddo
+
+
+  call run2()
+end
+subroutine run2
+  integer :: i,j, istate
+  print *,  'Phi_S'
+  do i=2,N_states
+    print *,  i, Phi_S(i)
+  enddo
+end
+
diff --git a/plugins/analyze_wf/phi_s.irp.f b/plugins/analyze_wf/phi_s.irp.f
new file mode 100644
index 00000000..12bdb970
--- /dev/null
+++ b/plugins/analyze_wf/phi_s.irp.f
@@ -0,0 +1,128 @@
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvalues,  (mo_tot_num, 2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvectors, (mo_tot_num, mo_tot_num, 2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Eigenvalues and eigenvectors of one_body_dm_mo_diff
+ END_DOC
+ integer                        :: i,j,istate
+ integer                        :: liwork, lwork, n, info
+ integer, allocatable           :: iwork(:)
+ double precision, allocatable  :: work(:)
+ 
+ 
+ one_body_dm_mo_diff_eigvectors(1:mo_tot_num, 1:mo_tot_num, 2:N_states) =&
+     one_body_dm_mo_diff(1:mo_tot_num, 1:mo_tot_num, 2:N_states)
+ 
+ n = mo_tot_num
+ lwork = 1+6*n + 2*n*n
+ liwork = 3 + 5*n
+ 
+ allocate(work(lwork), iwork(liwork))
+ 
+ lwork=-1
+ liwork=-1
+ istate=2
+ 
+ call dsyevd( 'V', 'U', mo_tot_num,                                  &
+     one_body_dm_mo_diff_eigvectors(1,1,istate),                     &
+     size(one_body_dm_mo_diff_eigvectors,1),                         &
+     one_body_dm_mo_diff_eigvalues(1,istate),                        &
+     work, lwork, iwork, liwork, info)
+ 
+ 
+ if (info /= 0) then
+   print *,  irp_here//' DSYEVD failed : ', info
+   stop 1
+ endif
+ lwork = int(work(1))
+ liwork = iwork(1)
+ deallocate(iwork,work)
+ 
+ allocate(work(lwork), iwork(liwork))
+ 
+ do istate=2,N_states
+   call dsyevd( 'V', 'U', mo_tot_num,                                &
+       one_body_dm_mo_diff_eigvectors(1,1,istate),                   &
+       size(one_body_dm_mo_diff_eigvectors,1),                       &
+       one_body_dm_mo_diff_eigvalues(1,istate),                      &
+       work, lwork, iwork, liwork, info)
+
+   if (info /= 0) then
+       print *,  irp_here//' DSYEVD failed : ', info
+       stop 1
+   endif
+
+ enddo
+
+ deallocate(iwork,work)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, transition_natorb, (ao_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Natural transition molecular orbitals
+ END_DOC
+
+ integer :: istate
+
+ do istate=2,N_states
+   call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
+       mo_coef, size(mo_coef,1),                                     &
+       one_body_dm_mo_diff_eigvectors(1,1,istate),                   &
+       size(one_body_dm_mo_diff_eigvectors,1), 0.d0,                 &
+       transition_natorb(1,1,istate), size(transition_natorb,1))
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, phi_s, (2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! 
+ END_DOC
+
+ integer :: i,istate
+ double precision, allocatable :: T(:,:), A(:,:), D(:,:)
+ double precision :: trace, norm
+ allocate(T(ao_num,ao_num), A(ao_num,ao_num), D(ao_num,ao_num))
+
+  do istate=2,N_states
+
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        S_half, size(S_half,1),                                      &
+        one_body_dm_ao_attachment(1,1,istate), size(one_body_dm_ao_attachment,1), 0.d0,&
+        T, size(T,1))
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        T, size(T,1),                                                &
+        S_half, size(S_half,1), 0.d0,                                &
+        A, size(A,1))
+!
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        S_half, size(S_half,1),                                      &
+        one_body_dm_ao_detachment(1,1,istate), size(one_body_dm_ao_detachment,1), 0.d0,&
+        T, size(T,1))
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        T, size(T,1),                                                &
+        S_half, size(S_half,1), 0.d0,                                &
+        D, size(D,1))
+
+    trace = 0.d0
+    do i=1,ao_num
+      trace = trace + A(i,i)
+    enddo
+    norm = 0.d0
+    do i=1,ao_num
+      norm = norm + D(i,i)
+    enddo
+    norm = 0.5d0*(norm + trace)
+
+    trace = 0.d0
+    do i=1,mo_tot_num
+      trace = trace + dsqrt(A(i,i)*D(i,i))
+    enddo
+    phi_s(istate) = trace/norm 
+  enddo
+
+END_PROVIDER
+
diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f
index f2630520..53fa74a1 100644
--- a/plugins/mrcepa0/mrcepa0_general.irp.f
+++ b/plugins/mrcepa0/mrcepa0_general.irp.f
@@ -38,12 +38,12 @@ subroutine run(N_st,energy)
     do while (delta_E > thresh_mrcc)
       iteration += 1
       print *,  '===============================================' 
-      print *,  'MRCEPA0 Iteration', iteration, '/', n_it_mrcc_max
+      print *,  'Iteration', iteration, '/', n_it_mrcc_max
       print *,  '===============================================' 
       print *,  ''
       E_old = sum(ci_energy_dressed(1:N_states))
       do i=1,N_st
-        call write_double(6,ci_energy_dressed(i),"MRCEPA0 energy")
+        call write_double(6,ci_energy_dressed(i),"Energy")
       enddo
       call diagonalize_ci_dressed(lambda)
       E_new = sum(ci_energy_dressed(1:N_states))
@@ -61,7 +61,7 @@ subroutine run(N_st,energy)
         exit
       endif
     enddo
-    call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy")
+    call write_double(6,ci_energy_dressed(1),"Final energy")
   endif
   energy(1:N_st) = ci_energy_dressed(1:N_st)
 end
diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py
index 780a7a91..56d79a4b 100755
--- a/scripts/compilation/qp_create_ninja.py
+++ b/scripts/compilation/qp_create_ninja.py
@@ -96,7 +96,8 @@ def ninja_create_env_variable(pwd_config_file):
         l_string.append(str_)
 
     lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB")
-    str_lib = " ".join([LIB, lib_lapack, EZFIO_LIB, ZMQ_LIB])
+    lib_gpi2 = get_compilation_option(pwd_config_file, "GPI2_LIB")
+    str_lib = " ".join([LIB, lib_lapack, lib_gpi2, EZFIO_LIB, ZMQ_LIB])
     l_string.append("LIB = {0} ".format(str_lib))
 
     l_string.append("")
diff --git a/scripts/compilation/read_compilation_cfg.py b/scripts/compilation/read_compilation_cfg.py
index f3ca01ce..43a64631 100755
--- a/scripts/compilation/read_compilation_cfg.py
+++ b/scripts/compilation/read_compilation_cfg.py
@@ -37,7 +37,8 @@ def get_compilation_option(pwd_cfg, flag_name):
 if __name__ == '__main__':
 
     qpackage_root = os.environ['QP_ROOT']
-    pwd_cfg = os.path.join(qpackage_root, "config/gfortran_example.cfg")
+    pwd_cfg = os.path.join(qpackage_root, "config/ifort_gpi2.cfg")
 
     print get_compilation_option(pwd_cfg, "FC")
     print get_compilation_option(pwd_cfg, "FCFLAGS")
+    print get_compilation_option(pwd_cfg, "GPI2_LIB")
diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py
index 4d61062e..ee44a1e1 100755
--- a/scripts/ezfio_interface/ei_handler.py
+++ b/scripts/ezfio_interface/ei_handler.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 """
-Welcom the ei_handler.
+Welcome to the ei_handler.
 We will create all the ezfio related stuff from a EZFIO.cfg file.
 
 Usage:
@@ -543,7 +543,7 @@ def create_ocaml_input(dict_ezfio_cfg, module_lower):
 
     template += ["open Qptypes;;",
                  "open Qputils;;",
-                 "open Core.Std;;",
+                 "open Core;;",
                  "",
                  "module {0} : sig".format(module_lower.capitalize())]
 
diff --git a/scripts/ezfio_interface/ezfio_generate_ocaml.py b/scripts/ezfio_interface/ezfio_generate_ocaml.py
index f36441b6..244f67a3 100755
--- a/scripts/ezfio_interface/ezfio_generate_ocaml.py
+++ b/scripts/ezfio_interface/ezfio_generate_ocaml.py
@@ -1,10 +1,10 @@
 #!/usr/bin/env python
 """
-This programme generate all the
-ocaml template needed by qp_edit
+This program generates all the
+OCaml templates needed by qp_edit
 
 You can see `ezfio_generate_provider.py`
-for an example of utilisation
+for an example.
 """
 
 import sys
@@ -170,7 +170,7 @@ class EZFIO_ocaml(object):
             else:
                 l_template += ["    {0:<30} : {1};".format(p, t.ocaml)]
 
-        l_template += ["  } with sexp",
+        l_template += ["  } [@@deriving sexp]",
                        ";;"]
 
         # ~#~#~#~#~#~ #
@@ -352,7 +352,7 @@ class EZFIO_ocaml(object):
 
         l_template = ['open Qputils;;',
                       'open Qptypes;;',
-                      'open Core.Std;;',
+                      'open Core;;',
                       '']
 
         for m in self.l_module_lower:
diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py
index 2a9bc087..7f4f30be 100755
--- a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py
+++ b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py
@@ -3,7 +3,7 @@
 convert output of gamess/GAU$$IAN to ezfio
 
 Usage:
-    qp_convert_output_to_ezfio.py <file.out> [--ezfio=<ezfio_directory>]
+    qp_convert_output_to_ezfio.py <file.out> [-o <ezfio_directory>]
 
 Option:
     file.out is the file to check (like gamess.out)
@@ -309,16 +309,14 @@ def write_ezfio(res, filename):
             array_l_max_block.append(l_max_block)
             array_z_remove.append(z_remove)
 
-            x = [[ filter(None,coef_n_zeta.split()) for coef_n_zeta in l.split('\n')] \
-               for l in array_party[1:] ]
-#            x = []
-#            for l in array_party[1:]:
-#              y = []
-#              for coef_n_zeta in l.split('\n'):
-#                z = coef_n_zeta.split()
-#                if z : y.append(z)
-#              x.append(y)
-#            matrix.append(x)
+            x = []
+            for l in array_party[1:]:
+              y = []
+              for coef_n_zeta in l.split('\n'):
+                z = coef_n_zeta.split()
+                if z : y.append(z)
+              x.append(y)
+            matrix.append(x)
         return (matrix, array_l_max_block, array_z_remove)
 
     def get_local_stuff(matrix):
@@ -400,8 +398,8 @@ def write_ezfio(res, filename):
         klocmax, m_coef, m_n, m_zeta = get_local_stuff(matrix)
         ezfio.pseudo_pseudo_klocmax = klocmax
 
-        ezfio.pseudo_pseudo_v_k = zip(*m_coef)
-        ezfio.pseudo_pseudo_n_k = zip(*m_n)
+        ezfio.pseudo_pseudo_v_k  = zip(*m_coef)
+        ezfio.pseudo_pseudo_n_k  = zip(*m_n)
         ezfio.pseudo_pseudo_dz_k = zip(*m_zeta)
 
         # ~#~#~#~#~#~#~#~#~ #
@@ -431,8 +429,8 @@ if __name__ == '__main__':
 
     file_ = get_full_path(arguments['<file.out>'])
 
-    if arguments["--ezfio"]:
-        ezfio_file = get_full_path(arguments["--ezfio"])
+    if arguments["-o"]:
+        ezfio_file = get_full_path(arguments["<ezfio_directory>"])
     else:
         ezfio_file = "{0}.ezfio".format(file_)
 
diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template
index 9d327124..55a35f83 100644
--- a/scripts/ezfio_interface/qp_edit_template
+++ b/scripts/ezfio_interface/qp_edit_template
@@ -4,7 +4,7 @@
 
 open Qputils
 open Qptypes
-open Core.Std
+open Core
 
 (** Interactive editing of the input.
 
diff --git a/src/AO_Basis/EZFIO.cfg b/src/AO_Basis/EZFIO.cfg
index 9e548514..07712a4a 100644
--- a/src/AO_Basis/EZFIO.cfg
+++ b/src/AO_Basis/EZFIO.cfg
@@ -16,7 +16,7 @@ interface: ezfio, provider
 
 [ao_prim_num_max]
 type: integer
-doc: number of primitive maximun
+doc: maximum number of primitives
 default: =maxval(ao_basis.ao_prim_num)
 interface: ezfio
 
@@ -76,3 +76,4 @@ size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 default: false
 
+
diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f
index edf48b25..83110293 100644
--- a/src/AO_Basis/ao_overlap.irp.f
+++ b/src/AO_Basis/ao_overlap.irp.f
@@ -129,3 +129,104 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
   !$OMP END PARALLEL DO
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]
+ implicit none
+ BEGIN_DOC
+! S^-1
+ END_DOC
+ call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1))
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]
+
+  BEGIN_DOC
+!   Matrix X = S^{-1/2} obtained by SVD
+  END_DOC
+
+  implicit none
+
+  integer                         :: num_linear_dependencies
+  integer                         :: LDA, LDC
+  double precision, allocatable   :: U(:,:),Vt(:,:), D(:)
+  integer                         :: info, i, j, k
+  double precision, parameter     :: threshold_overlap_AO_eigenvalues = 1.d-6
+ 
+  LDA = size(AO_overlap,1)
+  LDC = size(S_half_inv,1)
+ 
+  allocate(         &
+    U(LDC,AO_num),  &
+    Vt(LDA,AO_num), &
+    D(AO_num))
+
+  call svd(            &
+       AO_overlap,LDA, &
+       U,LDC,          &
+       D,              &
+       Vt,LDA,         &
+       AO_num,AO_num)
+
+  num_linear_dependencies = 0
+  do i=1,AO_num
+    print*,D(i)
+    if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
+      D(i) = 0.d0
+      num_linear_dependencies += 1
+    else
+      ASSERT (D(i) > 0.d0)
+      D(i) = 1.d0/sqrt(D(i))
+    endif
+    do j=1,AO_num
+      S_half_inv(j,i) = 0.d0
+    enddo
+  enddo
+  write(*,*) 'linear dependencies',num_linear_dependencies
+
+  do k=1,AO_num
+    if(D(k) /= 0.d0) then
+      do j=1,AO_num
+        do i=1,AO_num
+          S_half_inv(i,j) = S_half_inv(i,j) + U(i,k)*D(k)*Vt(k,j)
+        enddo
+      enddo
+    endif
+  enddo
+  
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num)  ]
+ implicit none
+ BEGIN_DOC
+ ! S^{1/2}
+ END_DOC
+
+  integer :: i,j,k
+  double precision, allocatable  :: U(:,:)
+  double precision, allocatable  :: Vt(:,:)
+  double precision, allocatable  :: D(:)
+  
+  allocate(U(ao_num,ao_num),Vt(ao_num,ao_num),D(ao_num))
+
+  call svd(ao_overlap,size(ao_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_num,ao_num)
+
+  do i=1,ao_num
+    D(i) = dsqrt(D(i))
+    do j=1,ao_num
+      S_half(j,i) = 0.d0
+    enddo
+  enddo
+
+  do k=1,ao_num
+      do j=1,ao_num
+        do i=1,ao_num
+          S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j)
+        enddo
+      enddo
+  enddo
+  
+  deallocate(U,Vt,D)
+
+END_PROVIDER
+
diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f
index f0f03fab..871c0ee6 100644
--- a/src/AO_Basis/aos.irp.f
+++ b/src/AO_Basis/aos.irp.f
@@ -9,14 +9,11 @@ BEGIN_PROVIDER [ integer, ao_num_align ]
    ao_num_align = align_double(ao_num)
 END_PROVIDER 
 
- BEGIN_PROVIDER [ integer, ao_prim_num_max ]
-&BEGIN_PROVIDER [ integer, ao_prim_num_max_align ]
+BEGIN_PROVIDER [ integer, ao_prim_num_max ]
  implicit none
  ao_prim_num_max = 0
  PROVIDE ezfio_filename
  call ezfio_get_ao_basis_ao_prim_num_max(ao_prim_num_max)
- integer :: align_double
- ao_prim_num_max_align = align_double(ao_prim_num_max)
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num_align,ao_prim_num_max) ]
diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg
index a68a61a5..f4c5d866 100644
--- a/src/Determinants/EZFIO.cfg
+++ b/src/Determinants/EZFIO.cfg
@@ -2,7 +2,7 @@
 type: Det_number_max
 doc: Max number of determinants in the wave function
 interface: ezfio,provider,ocaml
-default: 10000
+default: 1000000
 
 [N_det_max_property]
 type: Det_number_max
@@ -14,7 +14,7 @@ default: 10000
 type: Det_number_max
 doc: Maximum number of determinants diagonalized by Jacobi
 interface: ezfio,provider,ocaml
-default: 1000
+default: 2000
 
 [N_states]
 type: States_number
@@ -97,7 +97,7 @@ type: double precision
 size:  (determinants.n_det)
 
 [expected_s2]
-interface: ezfio,provider
+interface: ezfio
 doc: Expected value of S^2
 type: double precision
 
diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f
index e0764d96..4fa33daa 100644
--- a/src/Determinants/density_matrix.irp.f
+++ b/src/Determinants/density_matrix.irp.f
@@ -15,6 +15,32 @@
    enddo
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff, (mo_tot_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Difference of the one-body density matrix with respect to the ground state
+ END_DOC
+ integer :: i,j, istate
+
+ do istate=2,N_states
+  do j=1,mo_tot_num
+    do i=1,mo_tot_num
+      one_body_dm_mo_diff(i,j,istate) = &
+        one_body_dm_mo_alpha(i,j,istate) - one_body_dm_mo_alpha(i,j,1) + &
+        one_body_dm_mo_beta (i,j,istate) - one_body_dm_mo_beta (i,j,1) 
+    enddo
+  enddo
+  double precision :: trace
+  trace = 0.d0
+  do i=1,mo_tot_num
+    trace += one_body_dm_mo_diff(i,i,istate)
+  enddo
+  print *, irp_here, trace
+ enddo
+ 
+END_PROVIDER
+
+
  BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ]
  implicit none 
  integer :: i,j,ispin,istate
diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f
index ad955b97..9a1d4ee1 100644
--- a/src/Determinants/determinants.irp.f
+++ b/src/Determinants/determinants.irp.f
@@ -133,115 +133,6 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ]
 END_PROVIDER
 
 
- BEGIN_PROVIDER [ integer(bit_kind), psi_occ_pattern, (N_int,2,psi_det_size) ]
-&BEGIN_PROVIDER [ integer, N_occ_pattern ]
- implicit none
- BEGIN_DOC
-  ! array of the occ_pattern present in the wf
-  ! psi_occ_pattern(:,1,j) = jth occ_pattern of the wave function : represent all the single occupation
-  ! psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation
- END_DOC
- integer :: i,j,k
-
- ! create
- do i = 1, N_det
-  do k = 1, N_int
-   psi_occ_pattern(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i))
-   psi_occ_pattern(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i))
-  enddo
- enddo
-
- ! Sort
- integer, allocatable           :: iorder(:)
- integer*8, allocatable         :: bit_tmp(:)
- integer*8, external            :: occ_pattern_search_key
- integer(bit_kind), allocatable :: tmp_array(:,:,:)
- logical,allocatable            :: duplicate(:)
-
-
- allocate ( iorder(N_det), duplicate(N_det), bit_tmp(N_det), tmp_array(N_int,2,psi_det_size) )
-
- do i=1,N_det
-   iorder(i) = i
-   !$DIR FORCEINLINE
-   bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int)
- enddo
- call i8sort(bit_tmp,iorder,N_det)
- !DIR$ IVDEP
- do i=1,N_det
-  do k=1,N_int
-    tmp_array(k,1,i) = psi_occ_pattern(k,1,iorder(i))
-    tmp_array(k,2,i) = psi_occ_pattern(k,2,iorder(i))
-  enddo
-  duplicate(i) = .False.
- enddo
-
- i=1
- integer (bit_kind) :: occ_pattern_tmp
- do i=1,N_det
-  duplicate(i) = .False.
- enddo
-
- do i=1,N_det-1
-  if (duplicate(i)) then
-    cycle
-  endif
-  j = i+1
-  do while (bit_tmp(j)==bit_tmp(i))
-    if (duplicate(j)) then
-      j+=1
-      cycle
-    endif
-    duplicate(j) = .True.
-    do k=1,N_int
-      if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) &
-      .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then
-         duplicate(j) = .False.
-         exit
-      endif
-    enddo
-    j+=1
-    if (j>N_det) then
-      exit
-    endif
-  enddo
- enddo
-
- N_occ_pattern=0
- do i=1,N_det
-  if (duplicate(i)) then
-    cycle
-  endif
-  N_occ_pattern += 1
-  do k=1,N_int
-    psi_occ_pattern(k,1,N_occ_pattern) = tmp_array(k,1,i)
-    psi_occ_pattern(k,2,N_occ_pattern) = tmp_array(k,2,i)
-  enddo
- enddo
-
- deallocate(iorder,duplicate,bit_tmp,tmp_array)
-! !TODO DEBUG
-! integer :: s
-! do i=1,N_occ_pattern
-!   do j=i+1,N_occ_pattern
-!    s = 0
-!    do k=1,N_int
-!      if((psi_occ_pattern(k,1,j) /= psi_occ_pattern(k,1,i)).or. &
-!         (psi_occ_pattern(k,2,j) /= psi_occ_pattern(k,2,i))) then
-!         s=1
-!         exit
-!      endif
-!    enddo
-!    if ( s == 0 ) then
-!      print *,  'Error : occ ', j, 'already in wf'
-!      call debug_det(psi_occ_pattern(1,1,j),N_int)
-!      stop
-!    endif
-!   enddo
-! enddo
-! !TODO DEBUG
-END_PROVIDER 
-
 
 BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
   implicit none
diff --git a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
index f24e4974..1333753d 100644
--- a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
+++ b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
@@ -142,8 +142,8 @@ subroutine ao_bielec_integrals_in_map_collector
   integer(ZMQ_PTR), external     :: new_zmq_pull_socket
   integer(ZMQ_PTR)               :: zmq_socket_pull
   
-  integer*8                      :: control, accu
-  integer                        :: task_id, more, sze
+  integer*8                      :: control, accu, sze
+  integer                        :: task_id, more
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   zmq_socket_pull = new_zmq_pull_socket()
diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f
index 190615c3..996f8464 100644
--- a/src/Integrals_Bielec/map_integrals.irp.f
+++ b/src/Integrals_Bielec/map_integrals.irp.f
@@ -327,41 +327,50 @@ subroutine insert_into_mo_integrals_map(n_integrals,                 &
   call map_update(mo_integrals_map, buffer_i, buffer_values, n_integrals, thr)
 end
 
- BEGIN_PROVIDER [ integer, mo_integrals_cache_min ]
-&BEGIN_PROVIDER [ integer, mo_integrals_cache_max ]
+ BEGIN_PROVIDER [ integer*4, mo_integrals_cache_min ]
+&BEGIN_PROVIDER [ integer*4, mo_integrals_cache_max ]
+&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_min_8 ]
+&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_max_8 ]
  implicit none
  BEGIN_DOC
  ! Min and max values of the MOs for which the integrals are in the cache
  END_DOC
- mo_integrals_cache_min = max(1,elec_alpha_num - 31)
- mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+63)
+ mo_integrals_cache_min_8 = max(1_8,elec_alpha_num - 63_8)
+ mo_integrals_cache_max_8 = min(int(mo_tot_num,8),mo_integrals_cache_min+127_8)
+ mo_integrals_cache_min   = max(1,elec_alpha_num - 63)
+ mo_integrals_cache_max   = min(mo_tot_num,mo_integrals_cache_min+127)
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:64*64*64*64) ]
+BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*128_8) ]
  implicit none
  BEGIN_DOC
  ! Cache of MO integrals for fast access
  END_DOC
  PROVIDE mo_bielec_integrals_in_map
- integer                        :: i,j,k,l
- integer                        :: ii
+ integer*8                      :: i,j,k,l
+ integer*4                      :: i4,j4,k4,l4
+ integer*8                      :: ii
  integer(key_kind)              :: idx
  real(integral_kind)            :: integral
  FREE ao_integrals_cache
- !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
- do l=mo_integrals_cache_min,mo_integrals_cache_max
-   do k=mo_integrals_cache_min,mo_integrals_cache_max
-     do j=mo_integrals_cache_min,mo_integrals_cache_max
-       do i=mo_integrals_cache_min,mo_integrals_cache_max
+ !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral)
+ do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+   l4 = int(l,4)
+   do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+     k4 = int(k,4)
+     do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+       j4 = int(j,4)
+       do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+         i4 = int(i,4)
          !DIR$ FORCEINLINE
-         call bielec_integrals_index(i,j,k,l,idx)
+         call bielec_integrals_index(i4,j4,k4,l4,idx)
          !DIR$ FORCEINLINE
          call map_get(mo_integrals_map,idx,integral)
-         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)
+         ii = l-mo_integrals_cache_min_8
+         ii = ior( ishft(ii,7), k-mo_integrals_cache_min_8)
+         ii = ior( ishft(ii,7), j-mo_integrals_cache_min_8)
+         ii = ior( ishft(ii,7), i-mo_integrals_cache_min_8)
          mo_integrals_cache(ii) = integral
        enddo
      enddo
@@ -381,6 +390,7 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
   integer, intent(in)            :: i,j,k,l
   integer(key_kind)              :: idx
   integer                        :: ii
+  integer*8                      :: ii_8
   type(map_type), intent(inout)  :: map
   real(integral_kind)            :: tmp
   PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
@@ -388,18 +398,18 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
   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
+  if (iand(ii, -128) /= 0) then
     !DIR$ FORCEINLINE
     call bielec_integrals_index(i,j,k,l,idx)
     !DIR$ FORCEINLINE
     call map_get(map,idx,tmp)
     get_mo_bielec_integral = dble(tmp)
   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)
-    get_mo_bielec_integral = mo_integrals_cache(ii)
+    ii_8 = int(l,8)-mo_integrals_cache_min_8
+    ii_8 = ior( ishft(ii_8,7), int(k,8)-mo_integrals_cache_min_8)
+    ii_8 = ior( ishft(ii_8,7), int(j,8)-mo_integrals_cache_min_8)
+    ii_8 = ior( ishft(ii_8,7), int(i,8)-mo_integrals_cache_min_8)
+    get_mo_bielec_integral = mo_integrals_cache(ii_8)
   endif
 end
 
diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f
index d3e2eef9..11fef327 100644
--- a/src/MOGuess/H_CORE_guess.irp.f
+++ b/src/MOGuess/H_CORE_guess.irp.f
@@ -1,15 +1,9 @@
-program H_CORE_guess
+program H_CORE_guess_prog
   BEGIN_DOC
 ! Produce `H_core` MO orbital 
 ! output:  mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ
   END_DOC
   implicit none
   character*(64)                 :: label
-  label = "Guess"
-  call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
-                                     size(mo_mono_elec_integral,1),  &
-                                     size(mo_mono_elec_integral,2),label,1)
-  print *,  'save mos'
-  call save_mos
-  
+  call hcore_guess 
 end
diff --git a/src/MOGuess/h_core_guess_routine.irp.f b/src/MOGuess/h_core_guess_routine.irp.f
index 23899160..5b4ede91 100644
--- a/src/MOGuess/h_core_guess_routine.irp.f
+++ b/src/MOGuess/h_core_guess_routine.irp.f
@@ -7,8 +7,7 @@ subroutine hcore_guess
   label = "Guess"
   call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
                                      size(mo_mono_elec_integral,1),  &
-                                     size(mo_mono_elec_integral,2),label,1)
-  print *,  'save mos'
+                                     size(mo_mono_elec_integral,2),label,1,.false.)
   call save_mos
   SOFT_TOUCH mo_coef mo_label
 end
diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f
index 8962dd00..4bcd7221 100644
--- a/src/MO_Basis/mos.irp.f
+++ b/src/MO_Basis/mos.irp.f
@@ -35,7 +35,6 @@ END_PROVIDER
   END_DOC
   integer                        :: i, j
   double precision, allocatable  :: buffer(:,:)
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer
   logical                        :: exists
   PROVIDE ezfio_filename
   
@@ -153,16 +152,17 @@ BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ]
   endif
 END_PROVIDER
 
+
 subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
   implicit none
   BEGIN_DOC
   ! Transform A from the AO basis to the MO basis
   !
-  ! C.A_ao.Ct
+  ! Ct.A_ao.C
   END_DOC
   integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, intent(in)   :: A_ao(LDA_ao)
-  double precision, intent(out)  :: A_mo(LDA_mo)
+  double precision, intent(in)   :: A_ao(LDA_ao,ao_num)
+  double precision, intent(out)  :: A_mo(LDA_mo,mo_tot_num)
   double precision, allocatable  :: T(:,:)
   
   allocate ( T(ao_num_align,mo_tot_num) )
@@ -189,54 +189,48 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
   ! (S.C).A_mo.(S.C)t
   END_DOC
   integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, intent(in)   :: A_mo(LDA_mo)
-  double precision, intent(out)  :: A_ao(LDA_ao)
-  double precision, allocatable  :: T(:,:), SC(:,:)
+  double precision, intent(in)   :: A_mo(LDA_mo,mo_tot_num)
+  double precision, intent(out)  :: A_ao(LDA_ao,ao_num)
+  double precision, allocatable  :: T(:,:)
   
-  allocate ( SC(ao_num_align,mo_tot_num) )
   allocate ( T(mo_tot_num_align,ao_num) )
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
-  
-  call dgemm('N','N', ao_num, mo_tot_num, ao_num,                    &
-      1.d0, ao_overlap,size(ao_overlap,1),                           &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0, SC, ao_num_align)
   
   call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
-      1.d0, A_mo,LDA_mo,                                             &
-      SC, size(SC,1),                                      &
-      0.d0, T, mo_tot_num_align)
+      1.d0, A_mo,size(A_mo,1),                                       &
+      S_mo_coef, size(S_mo_coef,1),                                  &
+      0.d0, T, size(T,1))
   
   call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
-      1.d0, SC,size(SC,1),                                 &
-      T, mo_tot_num_align,                                           &
-      0.d0, A_ao, LDA_ao)
+      1.d0, S_mo_coef, size(S_mo_coef,1),                            &
+      T, size(T,1),                                                  &
+      0.d0, A_ao, size(A_ao,1))
   
-  deallocate(T,SC)
+  deallocate(T)
 end
 
 subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
   implicit none
   BEGIN_DOC
   ! Transform A from the MO basis to the S^-1 AO basis
+  ! Useful for density matrix
   END_DOC
   integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, intent(in)   :: A_mo(LDA_mo)
-  double precision, intent(out)  :: A_ao(LDA_ao)
+  double precision, intent(in)   :: A_mo(LDA_mo,mo_tot_num)
+  double precision, intent(out)  :: A_ao(LDA_ao,ao_num)
   double precision, allocatable  :: T(:,:)
   
   allocate ( T(mo_tot_num_align,ao_num) )
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
   
   call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
-      1.d0, A_mo,LDA_mo,                                             &
+      1.d0, A_mo,size(A_mo,1),                                       &
       mo_coef, size(mo_coef,1),                                      &
-      0.d0, T, mo_tot_num_align)
+      0.d0, T, size(T,1))
   
   call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
       1.d0, mo_coef,size(mo_coef,1),                                 &
-      T, mo_tot_num_align,                                           &
-      0.d0, A_ao, LDA_ao)
+      T, size(T,1),                                                  &
+      0.d0, A_ao, size(A_ao,1))
   
   deallocate(T)
 end
@@ -288,18 +282,17 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
   double precision, allocatable  :: T(:,:)
   
   allocate ( T(ao_num_align,ao_num) )
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
   
-  call dgemm('T','N', ao_num, ao_num, ao_num,                    &
+  call dgemm('T','N', ao_num, ao_num, ao_num,                        &
       1.d0,                                                          &
-      ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),  &
-      A_ao,LDA_ao,                                                   &
-      0.d0, T, ao_num_align)
+      ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),&
+      A_ao,size(A_ao,1),                                             &
+      0.d0, T, size(T,1))
   
-  call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0,          &
-      T, size(T,1),                                               &
+  call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0,                  &
+      T, size(T,1),                                                  &
       ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),&
-      0.d0, A, LDA)
+      0.d0, A, size(A,1))
   
   deallocate(T)
 end
diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f
index 750e3420..4221fce4 100644
--- a/src/MO_Basis/utils.irp.f
+++ b/src/MO_Basis/utils.irp.f
@@ -44,11 +44,12 @@ subroutine save_mos_truncated(n)
   
 end
 
-subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
+subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign,output)
   implicit none
   integer,intent(in)             :: n,m, sign
   character*(64), intent(in)     :: label
   double precision, intent(in)   :: matrix(n,m)
+  logical, intent(in)            :: output
   
   integer :: i,j
   double precision, allocatable  :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:)
@@ -76,22 +77,26 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
   mo_coef_new = mo_coef
   
   call lapack_diag(eigvalues,R,A,n,m)
-  write (output_mo_basis,'(A)')  'MOs are now **'//trim(label)//'**'
-  write (output_mo_basis,'(A)') ''
-  write (output_mo_basis,'(A)')  'Eigenvalues'
-  write (output_mo_basis,'(A)') '-----------'
-  write (output_mo_basis,'(A)')  ''
-  write (output_mo_basis,'(A)') '======== ================'
+  if (output) then
+    write (output_mo_basis,'(A)')  'MOs are now **'//trim(label)//'**'
+    write (output_mo_basis,'(A)') ''
+    write (output_mo_basis,'(A)')  'Eigenvalues'
+    write (output_mo_basis,'(A)') '-----------'
+    write (output_mo_basis,'(A)')  ''
+    write (output_mo_basis,'(A)') '======== ================'
+  endif
   if (sign == -1) then
     do i=1,m
       eigvalues(i) = -eigvalues(i)
     enddo
   endif
-  do i=1,m
-    write (output_mo_basis,'(I8,1X,F16.10)')  i,eigvalues(i)
-  enddo
-  write (output_mo_basis,'(A)') '======== ================'
-  write (output_mo_basis,'(A)')  ''
+  if (output) then
+    do i=1,m
+      write (output_mo_basis,'(I8,1X,F16.10)')  i,eigvalues(i)
+    enddo
+    write (output_mo_basis,'(A)') '======== ================'
+    write (output_mo_basis,'(A)')  ''
+  endif
   
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
   deallocate(A,mo_coef_new,R,eigvalues)
diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f
index 1a3d1d36..7e81756c 100644
--- a/src/Utils/LinearAlgebra.irp.f
+++ b/src/Utils/LinearAlgebra.irp.f
@@ -69,7 +69,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
   double precision, allocatable  :: U(:,:)
   double precision, allocatable  :: Vt(:,:)
   double precision, allocatable  :: D(:)
-  double precision, allocatable  :: S_half(:,:)
+  double precision, allocatable  :: S(:,:)
   !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D
   integer                        :: info, i, j
   
@@ -77,7 +77,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
     return
   endif
 
-  allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n))
+  allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n))
 
   call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
 
@@ -103,13 +103,13 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
 
 
   !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP SHARED(S_half,U,D,Vt,n,C,m) &
+  !$OMP SHARED(S,U,D,Vt,n,C,m) &
   !$OMP PRIVATE(i,j)
 
   !$OMP DO
   do j=1,n
     do i=1,n
-      S_half(i,j) = U(i,j)*D(j)
+      S(i,j) = U(i,j)*D(j)
     enddo
     do i=1,n
       U(i,j) = C(i,j)
@@ -119,8 +119,8 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
   
   !$OMP END PARALLEL
 
-  call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1))
-  deallocate (U, Vt, D, S_half)
+  call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1))
+  deallocate (U, Vt, D, S)
   
 end
 
@@ -210,19 +210,19 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   double precision, allocatable  :: U(:,:)
   double precision, allocatable  :: Vt(:,:)
   double precision, allocatable  :: D(:)
-  double precision, allocatable  :: S_half(:,:)
+  double precision, allocatable  :: S(:,:)
   integer                        :: info, i, j, k
   
   if (n < 2) then
     return
   endif
 
-  allocate(U(ldc,n),Vt(lda,n),S_half(lda,n),D(n))
+  allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n))
 
   call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
 
   !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP SHARED(S_half,U,D,Vt,n,C,m) &
+  !$OMP SHARED(S,U,D,Vt,n,C,m) &
   !$OMP PRIVATE(i,j,k)
 
   !$OMP DO
@@ -233,7 +233,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
       D(i) = 1.d0/dsqrt(D(i))
     endif
     do j=1,n
-      S_half(j,i) = 0.d0
+      S(j,i) = 0.d0
     enddo
   enddo
   !$OMP END DO
@@ -243,7 +243,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
       !$OMP DO
       do j=1,n
         do i=1,n
-          S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j)
+          S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j)
         enddo
       enddo
       !$OMP END DO NOWAIT
@@ -261,9 +261,9 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   
   !$OMP END PARALLEL
 
-  call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1))
+  call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1))
   
-  deallocate(U,Vt,S_half,D)
+  deallocate(U,Vt,S,D)
 end
 
 
diff --git a/src/Utils/constants.include.F b/src/Utils/constants.include.F
index 4974fd8e..7f1e71f3 100644
--- a/src/Utils/constants.include.F
+++ b/src/Utils/constants.include.F
@@ -1,6 +1,6 @@
 integer, parameter :: max_dim = 511
 integer, parameter :: SIMD_vector = 32
-integer, parameter :: N_int_max = 16
+integer, parameter :: N_int_max = 32
 
 double precision, parameter :: pi =  dacos(-1.d0)
 double precision, parameter :: sqpi =  dsqrt(dacos(-1.d0))