Merge branch 'dev' into cleaning_dft

2024-11-07 14:03:37 +01:00 · 2022-07-06 13:49:31 +02:00 · 2022-07-06 13:49:31 +02:00 · 4f610a49d6
commit 4f610a49d6
parent 17d20c6880 069d673e0a
68 changed files with 1851 additions and 1043 deletions
--- a/.drone.yml
+++ b/.drone.yml
@ -0,0 +1,52 @@
+---
+kind: pipeline
+type: docker
+name: gfortran-debug
+
+clone:
+  depth: 10
+
+steps:
+  - name: configure debug
+    image: scemama666/qp2_env
+    commands:
+    - ./configure -i all -c ./config/gfortran_debug.cfg
+    - bash -c "source quantum_package.rc ; exec qp_plugins download https://gitlab.com/scemama/qp_plugins_scemama"
+    - bash -c "source quantum_package.rc ; exec qp_plugins install champ"
+
+  - name: compile debug
+    image: scemama666/qp2_env
+    commands:
+    - bash -c "source quantum_package.rc ; exec ninja"
+
+  - name: testing debug
+    image: scemama666/qp2_env
+    commands:
+    - bash -c "source quantum_package.rc ; TRAVIS=1 exec qp_test -a"
+
+  - name: configure fast
+    image: scemama666/qp2_env
+    commands:
+    - ./configure -c ./config/gfortran_avx.cfg
+
+  - name: compile fast
+    image: scemama666/qp2_env
+    commands:
+    - bash -c "source quantum_package.rc ; exec ninja"
+
+  - name: testing fast
+    image: scemama666/qp2_env
+    commands:
+    - bash -c "source quantum_package.rc ; exec qp_test -a"
+
+  - name: notify
+    image: drillster/drone-email
+    settings:
+      host:
+        from_secret: hostname # irsamc.ups-tlse.fr
+      from:
+        from_secret: from # drone@irssv7.ups-tlse.fr
+      recipients:
+        from_secret: recipients # scemama@irsamc.ups-tlse.fr
+      when:
+        status: [changed, failure]
--- a/README.md
+++ b/README.md
@ -1,6 +1,7 @@
 # Quantum Package 2.2

-<img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250">
+<!--- img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250" --->
+<img src="https://trex-coe.eu/sites/default/files/styles/responsive_no_crop/public/2021-12/Risorsa%2014_0.png"  width="250">


 [![DOI](https://zenodo.org/badge/167513335.svg)](https://zenodo.org/badge/latestdoi/167513335)
@ -44,9 +45,19 @@ https://arxiv.org/abs/1902.08154

 # Credits

+* [TREX Center of Excellence](https://trex-coe.eu)
+* [ERC PTEROSOR](https://lcpq.github.io/PTEROSOR)
 * [CNRS](http://www.cnrs.fr)
 * [Laboratoire de Chimie et Physique Quantiques](http://lcpq.ups-tlse.fr)
 * [Laboratoire de Chimie Théorique](http://www.lct.jussieu.fr)
 * [Argonne  Leadership Computing  Facility](http://alcf.anl.gov)
 * [CALMIP](https://www.calmip.univ-toulouse.fr)

+
+------------------------------
+
+<img src="https://lcpq.github.io/PTEROSOR/img/ERC.png" width="200" />
+
+[TREX: Targeting Real Chemical Accuracy at the Exascale](https://trex-coe.eu) project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content.
+
+[PTEROSOR](https://lcpq.github.io/PTEROSOR) project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 863481).
--- a/bin/qp_convert_output_to_ezfio
+++ b/bin/qp_convert_output_to_ezfio
@ -146,6 +146,17 @@ def write_ezfio(res, filename):
    ezfio.set_ao_basis_ao_nucl(at)
    ezfio.set_ao_basis_ao_prim_num(num_prim)
    ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
+    try:
+      normf = res.normf
+      if normf == 0:
+        ezfio.set_ao_basis_ao_normalized(True)
+      elif normf == 1:
+        ezfio.set_ao_basis_ao_normalized(False)
+      else:
+        print("BUG in NORMF")
+        sys.exit(0)
+    except AttributeError:
+      ezfio.set_ao_basis_ao_normalized(True)

    # ~#~#~#~#~#~#~ #
    # P a r s i n g #
@ -224,7 +235,7 @@ def write_ezfio(res, filename):
          exponent += [p.expo for p in b.prim]
          ang_mom.append(str.count(s, "z"))
          shell_prim_num.append(len(b.prim))
-          shell_index += [nshell_tot+1] * len(b.prim)
+          shell_index += [nshell_tot] * len(b.prim)

    # ~#~#~#~#~ #
    # W r i t e #
--- a/bin/qp_test
+++ b/bin/qp_test
@ -60,19 +60,14 @@ def main(arguments):
        print("Running tests for %s"%(bats_file))
        print("")
        if arguments["-v"]:
-            p = None
            if arguments["TEST"]:
                test = "export TEST=%s ; "%arguments["TEST"]
            else:
                test = ""
-            try:
-                os.system(test+" python3 bats_to_sh.py "+bats_file+
+            os.system(test+" python3 bats_to_sh.py "+bats_file+
                                    "| bash")
-            except:
-                if p:
-                    p.terminate()
        else:
-            subprocess.check_call(["bats", bats_file], env=os.environ)
+            subprocess.check_call(["bats", "--verbose-run", "--trace", bats_file], env=os.environ)



--- a/config/ifort_2019_avx.cfg
+++ b/config/ifort_2019_avx.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel 
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED

--- a/config/ifort_2019_avx_mpi.cfg
+++ b/config/ifort_2019_avx_mpi.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL -DSET_NESTED

--- a/config/ifort_2019_mpi_rome.cfg
+++ b/config/ifort_2019_mpi_rome.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED

--- a/config/ifort_2019_rome.cfg
+++ b/config/ifort_2019_rome.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED

--- a/config/ifort_2019_sse4.cfg
+++ b/config/ifort_2019_sse4.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED

--- a/config/ifort_2019_sse4_mpi.cfg
+++ b/config/ifort_2019_sse4_mpi.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL -DSET_NESTED

--- a/config/ifort_2019_xHost.cfg
+++ b/config/ifort_2019_xHost.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=64 -DINTEL -DSET_NESTED

--- a/config/ifort_2021_avx.cfg
+++ b/config/ifort_2021_avx.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL 

--- a/config/ifort_2021_avx_mpi.cfg
+++ b/config/ifort_2021_avx_mpi.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL 

--- a/config/ifort_2021_mpi_rome.cfg
+++ b/config/ifort_2021_mpi_rome.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL 

--- a/config/ifort_2021_rome.cfg
+++ b/config/ifort_2021_rome.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL 

--- a/config/ifort_2021_sse4.cfg
+++ b/config/ifort_2021_sse4.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DINTEL 

--- a/config/ifort_2021_sse4_mpi.cfg
+++ b/config/ifort_2021_sse4_mpi.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL 

--- a/config/ifort_2021_xHost.cfg
+++ b/config/ifort_2021_xHost.cfg
@ -7,7 +7,7 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
+LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=64 -DINTEL 

--- a/4
+++ b/4
@ -281,8 +281,8 @@ EOF

            execute << EOF
              cd "\${QP_ROOT}"/external
-              tar -zxf qp2-dependencies/bats-v1.1.0.tar.gz
-              ( cd bats-core-1.1.0/ ; ./install.sh \${QP_ROOT})
+              tar -zxf qp2-dependencies/bats-v1.7.0.tar.gz
+              ( cd bats-core-1.7.0/ ; ./install.sh \${QP_ROOT})
 EOF

    else
--- a/etc/qp.rc
+++ b/etc/qp.rc
@ -80,8 +80,6 @@ function qp()
      if [[ -d $NAME ]] ; then
        [[ -d $EZFIO_FILE ]] && ezfio unset_file
        ezfio set_file $NAME
-      else
-        qp_create_ezfio -h | more
      fi
      unset _ARGS
      ;;
--- a/external/qp2-dependencies
+++ b/external/qp2-dependencies
@ -1 +1 @@
-Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0
+Subproject commit 242151e03d1d6bf042387226431d82d35845686a
--- a/ocaml/Command_line.ml
+++ b/ocaml/Command_line.ml
@ -1,3 +1,5 @@
+exception Error of string
+
 type short_opt     = char
 type long_opt      = string 
 type optional      = Mandatory | Optional 
@ -181,15 +183,16 @@ let set_specs specs_in =
    Getopt.parse_cmdline cmd_specs (fun x -> anon_args := !anon_args @ [x]);

    if show_help () then
-        (help () ; exit 0);
+        help ()
+    else
+        (* Check that all mandatory arguments are set *)
+        List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
+        |> List.iter (fun x -> 
+            match get x.long with
+            | Some _ -> ()
+            | None -> raise (Error ("--"^x.long^" option is missing."))
+            )

-    (* Check that all mandatory arguments are set *)
-    List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
-    |> List.iter (fun x -> 
-        match get x.long with
-        | Some _ -> ()
-        | None -> failwith ("Error: --"^x.long^" option is missing.")
-        )
 ;;


--- a/ocaml/Command_line.mli
+++ b/ocaml/Command_line.mli
@ -59,6 +59,8 @@ let () =
 *)


+exception Error of string
+
 type short_opt = char

 type long_opt = string
--- a/ocaml/Molecule.ml
+++ b/ocaml/Molecule.ml
@ -101,7 +101,7 @@ let to_string_general ~f m =
  |> String.concat "\n"

 let to_string =
-  to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x)
+  to_string_general ~f:(fun x -> Atom.to_string ~units:Units.Angstrom x)

 let to_xyz =
  to_string_general ~f:Atom.to_xyz
@ -113,7 +113,7 @@ let of_xyz_string
    s =
  let l = String_ext.split s ~on:'\n'
       |> List.filter (fun x -> x <> "")
-       |> list_map (fun x -> Atom.of_string units x)
+       |> list_map (fun x -> Atom.of_string ~units x)
  in
  let ne = ( get_charge {
        nuclei=l ;
--- a/ocaml/Qputils.ml
+++ b/ocaml/Qputils.ml
@ -56,3 +56,7 @@ let string_of_string s = s
 let list_map f l =
  List.rev_map f l
  |> List.rev
+
+let socket_convert socket =
+    ((Obj.magic (Obj.repr socket)) : [ `Xsub ] Zmq.Socket.t )
+
--- a/ocaml/qp_create_ezfio.ml
+++ b/ocaml/qp_create_ezfio.ml
@ -677,6 +677,7 @@ let run ?o b au c d m p cart xyz_file =

 let () =

+  try (

  let open Command_line in
  begin
@ -734,7 +735,7 @@ If a file with the same name as the basis set exists, this file will be read.  O

  let basis =
    match Command_line.get "basis" with
-    | None -> assert false
+    | None -> ""
    | Some x -> x
  in

@ -773,10 +774,14 @@ If a file with the same name as the basis set exists, this file will be read.  O

  let xyz_filename =
    match Command_line.anon_args () with
-    | [x] -> x
-    | _ -> (Command_line.help () ; failwith "input file is missing")
+    | []  -> failwith "input file is missing"
+    | x::_ -> x
  in

  run ?o:output basis au charge dummy multiplicity pseudo cart xyz_filename
+  )
+  with
+  | Failure txt  -> Printf.eprintf "Fatal error: %s\n%!" txt
+  | Command_line.Error txt  -> Printf.eprintf "Command line error: %s\n%!" txt


--- a/ocaml/qp_run.ml
+++ b/ocaml/qp_run.ml
@ -110,7 +110,7 @@ let run slave ?prefix exe ezfio_file =
  let task_thread =
     let thread =
      Thread.create ( fun () ->
-         TaskServer.run port_number )
+         TaskServer.run ~port:port_number )
     in
     thread ();
  in
--- a/ocaml/qp_tunnel.ml
+++ b/ocaml/qp_tunnel.ml
@ -2,7 +2,7 @@ open Qputils
 open Qptypes

 type ezfio_or_address = EZFIO of string | ADDRESS of string
-type req_or_sub = REQ | SUB 
+type req_or_sub = REQ | SUB

 let localport = 42379

@ -29,7 +29,7 @@ let () =
  end;

  let arg =
-    let x = 
+    let x =
      match Command_line.anon_args () with
      | [x] -> x
      | _ -> begin
@ -44,7 +44,7 @@ let () =
  in


-  let localhost = 
+  let localhost =
     Lazy.force TaskServer.ip_address
  in

@ -52,28 +52,28 @@ let () =
  let long_address  =
    match arg with
    | ADDRESS x -> x
-    | EZFIO   x -> 
-	let ic = 
+    | EZFIO   x ->
+	let ic =
 	  Filename.concat  (Qpackage.ezfio_work x)  "qp_run_address"
          |> open_in
 	in
-        let result = 
+        let result =
          input_line ic
          |> String.trim
        in
        close_in ic;
        result
  in
-  
+
  let protocol, address, port =
    match String.split_on_char ':' long_address with
    | t :: a :: p :: [] -> t, a, int_of_string p
-    | _ -> failwith @@ 
+    | _ -> failwith @@
           Printf.sprintf "%s : Malformed address" long_address
  in


-  let zmq_context = 
+  let zmq_context =
    Zmq.Context.create ()
  in

@ -105,10 +105,10 @@ let () =


  let create_socket  sock_type  bind_or_connect  addr =
-    let socket = 
+    let socket =
      Zmq.Socket.create zmq_context sock_type
    in
-    let () = 
+    let () =
      try
        bind_or_connect socket addr
      with
@ -131,37 +131,64 @@ let () =
  Sys.set_signal  Sys.sigint   handler;


-  let new_thread req_or_sub addr_in addr_out =
+  let new_thread_req addr_in addr_out =
    let socket_in, socket_out =
-      match req_or_sub with
-      | REQ ->
          create_socket  Zmq.Socket.router  Zmq.Socket.bind     addr_in,
          create_socket  Zmq.Socket.dealer  Zmq.Socket.connect  addr_out
-      | SUB ->
-          create_socket  Zmq.Socket.sub  Zmq.Socket.connect  addr_in,
-          create_socket  Zmq.Socket.pub  Zmq.Socket.bind     addr_out
    in

-    if req_or_sub = SUB then
-        Zmq.Socket.subscribe  socket_in  "";

-
-
-   let action_in = 
-      match req_or_sub with
-        | REQ -> (fun () -> Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out)
-        | SUB -> (fun () -> Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out)
+    let action_in =
+        fun () -> Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out
    in

-    let action_out = 
-      match req_or_sub with
-        | REQ -> (fun () -> Zmq.Socket.recv_all  socket_out |> Zmq.Socket.send_all  socket_in )
-        | SUB -> (fun () -> () )
+    let action_out =
+        fun () -> Zmq.Socket.recv_all  socket_out |> Zmq.Socket.send_all  socket_in
    in

    let pollitem =
      Zmq.Poll.mask_of
-        [| (socket_in,  Zmq.Poll.In) ; (socket_out, Zmq.Poll.In) |]
+        [| (socket_convert socket_in,  Zmq.Poll.In) ; (socket_convert socket_out, Zmq.Poll.In) |]
+    in
+
+    while !run_status do
+
+        let polling =
+          Zmq.Poll.poll  ~timeout:1000  pollitem
+        in
+
+        match polling with
+          | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |]     -> ( action_out () ; action_in () )
+          | [| _                ; Some Zmq.Poll.In |]     ->   action_out ()
+          | [| Some Zmq.Poll.In ; _                |]     ->   action_in  ()
+          | _                    -> ()
+    done;
+
+    Zmq.Socket.close  socket_in;
+    Zmq.Socket.close  socket_out;
+  in
+
+  let new_thread_sub addr_in addr_out =
+    let socket_in, socket_out =
+          create_socket  Zmq.Socket.sub  Zmq.Socket.connect  addr_in,
+          create_socket  Zmq.Socket.pub  Zmq.Socket.bind     addr_out
+    in
+
+    Zmq.Socket.subscribe  socket_in  "";
+
+
+
+   let action_in =
+        fun () -> Zmq.Socket.recv_all  socket_in  |> Zmq.Socket.send_all  socket_out
+    in
+
+    let action_out =
+        fun () -> ()
+    in
+
+    let pollitem =
+      Zmq.Poll.mask_of
+        [| (socket_convert socket_in,  Zmq.Poll.In) ; (socket_convert socket_out, Zmq.Poll.In) |]
    in


@ -173,8 +200,8 @@ let () =

        match polling with
          | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |]     -> ( action_out () ; action_in () )
-          | [| _                ; Some Zmq.Poll.In |]     ->   action_out () 
-          | [| Some Zmq.Poll.In ; _                |]     ->   action_in  () 
+          | [| _                ; Some Zmq.Poll.In |]     ->   action_out ()
+          | [| Some Zmq.Poll.In ; _                |]     ->   action_in  ()
          | _                    -> ()
    done;

@ -193,8 +220,8 @@ let () =
      Printf.sprintf "tcp://*:%d"  localport
    in

-    let f () = 
-      new_thread  REQ  addr_in  addr_out
+    let f () =
+      new_thread_req  addr_in  addr_out
    in

    (Thread.create f) ()
@ -211,8 +238,8 @@ let () =
      Printf.sprintf "tcp://*:%d" (localport+2)
    in

-    let f () = 
-      new_thread  REQ  addr_in  addr_out
+    let f () =
+      new_thread_req  addr_in  addr_out
    in
    (Thread.create f) ()
  in
@ -227,8 +254,8 @@ let () =
      Printf.sprintf "tcp://*:%d" (localport+1)
    in

-    let f () = 
-      new_thread  SUB  addr_in  addr_out
+    let f () =
+      new_thread_sub  addr_in  addr_out
    in
    (Thread.create f) ()
  in
@ -236,7 +263,7 @@ let () =


  let input_thread =
-    let f () = 
+    let f () =
      let addr_out =
        match arg with
        | EZFIO _ -> None
@ -248,22 +275,22 @@ let () =
        Printf.sprintf "tcp://*:%d"  (localport+9)
      in

-      let socket_in = 
+      let socket_in =
        create_socket  Zmq.Socket.rep  Zmq.Socket.bind     addr_in
      in

      let socket_out =
-        match addr_out with 
+        match addr_out with
        | Some addr_out -> Some (
            create_socket  Zmq.Socket.req  Zmq.Socket.connect  addr_out)
        | None -> None
      in

-      let temp_file = 
+      let temp_file =
        Filename.temp_file "qp_tunnel" ".tar.gz"
      in

-      let get_ezfio_filename () = 
+      let get_ezfio_filename () =
        match arg with
        | EZFIO x -> x
        | ADDRESS _ ->
@ -277,9 +304,9 @@ let () =
          end
      in

-      let get_input () = 
+      let get_input () =
        match arg with
-        | EZFIO x -> 
+        | EZFIO x ->
          begin
            Printf.sprintf "tar --exclude=\"*.gz.*\" -zcf %s %s" temp_file x
            |> Sys.command |> ignore;
@ -291,11 +318,11 @@ let () =
            in
            ignore @@ Unix.lseek  fd  0  Unix.SEEK_SET ;
            let bstr =
-                Unix.map_file  fd  Bigarray.char 
+                Unix.map_file  fd  Bigarray.char
                  Bigarray.c_layout false [| len |]
                |> Bigarray.array1_of_genarray
            in
-            let result = 
+            let result =
              String.init len (fun i -> bstr.{i}) ;
            in
            Unix.close fd;
@ -313,7 +340,7 @@ let () =
          end
      in

-      let () = 
+      let () =
        match socket_out with
        | None -> ()
        | Some socket_out ->
@ -329,7 +356,7 @@ let () =
          | ADDRESS _ ->
            begin
              Printf.printf "Getting input... %!";
-              let ezfio_filename = 
+              let ezfio_filename =
                get_ezfio_filename ()
              in
              Printf.printf "%s%!" ezfio_filename;
@ -343,7 +370,7 @@ let () =
              |> Sys.command |> ignore ;
              let oc =
                Filename.concat  (Qpackage.ezfio_work ezfio_filename)  "qp_run_address"
-                |> open_out 
+                |> open_out
              in
              Printf.fprintf oc "tcp://%s:%d\n"  localhost  localport;
              close_out oc;
@ -359,9 +386,9 @@ let () =
      let action () =
        match Zmq.Socket.recv socket_in with
        | "get_input" -> get_input ()
-                        |> Zmq.Socket.send socket_in 
+                        |> Zmq.Socket.send socket_in
        | "get_ezfio_filename" -> get_ezfio_filename ()
-                        |> Zmq.Socket.send socket_in 
+                        |> Zmq.Socket.send socket_in
        | "test" -> Zmq.Socket.send socket_in "OK"
        | x -> Printf.sprintf "Message '%s' not understood" x
               |> Zmq.Socket.send socket_in
@ -372,7 +399,7 @@ On remote hosts, create ssh tunnel using:
 ssh -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d %s &
 Or from this host connect to clients using:
 ssh -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d <host> &
-%!" 
+%!"
                (port  ) localhost (localport  )
                (port+1) localhost (localport+1)
                (port+2) localhost (localport+2)
@ -392,12 +419,12 @@ Or from this host connect to clients using:
          match polling.(0) with
            | Some Zmq.Poll.In     -> action ()
            | None                 -> ()
-            | Some Zmq.Poll.In_out 
+            | Some Zmq.Poll.In_out
            | Some Zmq.Poll.Out    -> ()

      done;

-      let () = 
+      let () =
        match socket_out with
        | Some socket_out -> Zmq.Socket.close  socket_out
        | None -> ()
@ -415,7 +442,7 @@ Or from this host connect to clients using:
  Thread.join ocaml_thread;
  Zmq.Context.terminate zmq_context;
  Printf.printf "qp_tunnel exited properly.\n"
-  
+



--- a/scripts/cipsi_save.sh
+++ b/scripts/cipsi_save.sh
@ -0,0 +1,27 @@
+#!/bin/bash
+#
+# This script runs a CIPSI calculation as a sequence of single CIPSI iterations.
+# After each iteration, the EZFIO directory is saved.
+#
+# Usage: cipsi_save [EZFIO_FILE] [NDET]
+#
+# Example: cipsi_save file.ezfio 10000
+
+EZ=$1
+NDETMAX=$2
+
+qp set_file ${EZ}
+qp reset -d
+qp set determinants read_wf true
+declare -i NDET
+NDET=1
+while [[ ${NDET} -lt ${NDETMAX} ]]
+do
+        NDET=$(($NDET + $NDET))
+        qp set determinants n_det_max $NDET
+        qp run fci > ${EZ}.out
+        NDET=$(qp get determinants n_det)
+	mv ${EZ}.out ${EZ}.${NDET}.out
+        cp -r ${EZ} ${EZ}.${NDET}
+done
+
--- a/src/cipsi/EZFIO.cfg
+++ b/src/cipsi/EZFIO.cfg
@ -34,3 +34,9 @@ doc: Maximum number of excitation for beta determinants with respect to the Hart
 interface: ezfio,ocaml,provider
 default: -1

+[twice_hierarchy_max]
+type: integer
+doc: Twice the maximum hierarchy parameter (excitation degree plus half the seniority number). Using -1 selects all determinants
+interface: ezfio,ocaml,provider
+default: -1
+
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@ -290,9 +290,13 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
      call set_multiple_levels_omp(.False.)


-      print '(A)', '========== ======================= ===================== ===================== ==========='
-      print '(A)', ' Samples          Energy                Variance               Norm^2            Seconds'
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      ! old
+      !print '(A)', '========== ======================= ===================== ===================== ==========='
+      !print '(A)', ' Samples          Energy                Variance               Norm^2            Seconds'
+      !print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
+      print '(A)', ' Samples           Energy             PT2            Variance          Norm^2       Convergence    Seconds'
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='

      PROVIDE global_selection_buffer

@ -316,7 +320,10 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
      call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2')
      call set_multiple_levels_omp(.True.)

-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      ! old
+      !print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
+

      do k=1,N_states
        pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate)
@ -414,6 +421,17 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
  double precision :: rss
  double precision, external :: memory_of_double, memory_of_int

+  character(len=20) :: format_str1, str_error1, format_str2, str_error2
+  character(len=20) :: format_str3, str_error3, format_str4, str_error4
+  character(len=20) :: format_value1, format_value2, format_value3, format_value4
+  character(len=20) :: str_value1, str_value2, str_value3, str_value4
+  character(len=20) :: str_conv
+  double precision  :: value1, value2, value3, value4
+  double precision  :: error1, error2, error3, error4
+  integer           :: size1,size2,size3,size4
+
+  double precision :: conv_crit
+
  sending =.False.

  rss  = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2)
@ -537,14 +555,60 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_

          if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
            time1 = time
-            print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1)', c, &
-              pt2_data     % pt2(pt2_stoch_istate) +E, &
-              pt2_data_err % pt2(pt2_stoch_istate), &
-              pt2_data     % variance(pt2_stoch_istate), &
-              pt2_data_err % variance(pt2_stoch_istate), &
-              pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
-              pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
-              time-time0
+
+            value1 =  pt2_data     % pt2(pt2_stoch_istate) + E
+            error1 =  pt2_data_err % pt2(pt2_stoch_istate)           
+            value2 =  pt2_data     % pt2(pt2_stoch_istate)
+            error2 =  pt2_data_err % pt2(pt2_stoch_istate)           
+            value3 =  pt2_data     % variance(pt2_stoch_istate)           
+            error3 =  pt2_data_err % variance(pt2_stoch_istate)           
+            value4 =  pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate)           
+            error4 =  pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate)           
+            
+            ! Max size of the values (FX.Y) with X=size
+            size1 = 15
+            size2 = 12
+            size3 = 12
+            size4 = 12
+
+            ! To generate the format: number(error)
+            call format_w_error(value1,error1,size1,8,format_value1,str_error1)
+            call format_w_error(value2,error2,size2,8,format_value2,str_error2)
+            call format_w_error(value3,error3,size3,8,format_value3,str_error3)
+            call format_w_error(value4,error4,size4,8,format_value4,str_error4)
+
+            ! value > string with the right format
+            write(str_value1,'('//format_value1//')') value1
+            write(str_value2,'('//format_value2//')') value2
+            write(str_value3,'('//format_value3//')') value3
+            write(str_value4,'('//format_value4//')') value4
+
+            ! Convergence criterion
+            conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
+                  (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
+            write(str_conv,'(G10.3)') conv_crit
+
+            write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,&
+            adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),&
+            adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),&
+            adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),&
+            adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),&
+            adjustl(str_conv),& 
+            time-time0
+
+            ! Old print
+            !print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1,1pE16.6,1pE16.6)', c, &
+            !  pt2_data     % pt2(pt2_stoch_istate) +E, &
+            !  pt2_data_err % pt2(pt2_stoch_istate), &
+            !  pt2_data     % variance(pt2_stoch_istate), &
+            !  pt2_data_err % variance(pt2_stoch_istate), &
+            !  pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            !  pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            !  time-time0, &
+            !  pt2_data     % pt2(pt2_stoch_istate), &
+            !  dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
+            !      (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
+
            if (stop_now .or. (                                      &
                  (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
                  (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then
--- a/src/cipsi/run_selection_slave.irp.f
+++ b/src/cipsi/run_selection_slave.irp.f
@ -61,10 +61,14 @@ subroutine run_selection_slave(thread,iproc,energy)
        if (N /= buf%N) then
          print *, 'N=', N
          print *, 'buf%N=', buf%N
-          print *, 'bug in ', irp_here
-          stop '-1'
+          print *, 'In ', irp_here, ': N /= buf%N'
+          stop -1
        end if
      end if
+      if (i_generator > N_det_generators) then
+          print *, 'In ', irp_here, ': i_generator > N_det_generators'
+          stop -1
+      endif
      call select_connected(i_generator,energy,pt2_data,buf,subset,pt2_F(i_generator))
    endif

--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@ -195,7 +195,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d

  integer                        :: l_a, nmax, idx
  integer, allocatable           :: indices(:), exc_degree(:), iorder(:)
-  double precision, parameter :: norm_thr = 1.d-16
+
+  ! Removed to avoid introducing determinants already presents in the wf
+  !double precision, parameter :: norm_thr = 1.d-16
+
  allocate (indices(N_det),                                          &
      exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))

@ -215,10 +218,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
      i = psi_bilinear_matrix_rows(l_a)
      if (nt + exc_degree(i) <= 4) then
        idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a))
-        if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
+        ! Removed to avoid introducing determinants already presents in the wf
+        !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
          indices(k) = idx
          k=k+1
-        endif
+        !endif
      endif
    enddo
  enddo
@ -242,10 +246,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
        idx = psi_det_sorted_order(                                  &
            psi_bilinear_matrix_order(                               &
            psi_bilinear_matrix_transp_order(l_a)))
-        if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
+        ! Removed to avoid introducing determinants already presents in the wf
+        !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
          indices(k) = idx
          k=k+1
-        endif
+        !endif
      endif
    enddo
  enddo
@ -253,8 +258,6 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
  deallocate(exc_degree)
  nmax=k-1

-  call isort_noidx(indices,nmax)
-
  ! Start with 32 elements. Size will double along with the filtering.
  allocate(preinteresting(0:32), prefullinteresting(0:32),     &
      interesting(0:32), fullinteresting(0:32))
@ -566,6 +569,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
  double precision, external :: diag_H_mat_elem_fock
  double precision :: E_shift
  double precision :: s_weight(N_states,N_states)
+  logical, external :: is_in_wavefunction
  PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs
  do jstate=1,N_states
    do istate=1,N_states
@ -707,6 +711,25 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
        if (do_cycle) cycle
      endif

+      if (twice_hierarchy_max >= 0) then
+        s = 0
+        do k=1,N_int
+          s = s + popcnt(ieor(det(k,1),det(k,2)))
+        enddo
+        if ( mod(s,2)>0 ) stop 'For now, hierarchy CI is defined only for an even number of electrons'
+        if (excitation_ref == 1) then
+          call get_excitation_degree(HF_bitmask,det(1,1),degree,N_int)
+        else if (excitation_ref == 2) then
+          stop 'For now, hierarchy CI is defined only for a single reference determinant'
+!         do k=1,N_dominant_dets_of_cfgs
+!           call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int)
+!         enddo
+        endif
+        integer :: twice_hierarchy
+        twice_hierarchy = degree + s/2
+        if (twice_hierarchy > twice_hierarchy_max) cycle
+      endif
+
      Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)

      w = 0d0
@ -777,7 +800,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d

        alpha_h_psi = mat(istate, p1, p2)

-        pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate)
+        do k=1,N_states
+          pt2_data % overlap(k,istate) = pt2_data % overlap(k,istate) + coef(k) * coef(istate)
+        end do
        pt2_data % variance(istate)  = pt2_data % variance(istate) + alpha_h_psi * alpha_h_psi
        pt2_data % pt2(istate)       = pt2_data % pt2(istate)      + e_pert(istate)

@ -828,8 +853,27 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
            endif

        end select
+
+        ! To force the inclusion of determinants with a positive pt2 contribution
+        if (e_pert(istate) > 1d-8) then
+          w = -huge(1.0)
+        endif
+
      end do

+!!!BEGIN_DEBUG
+!      ! To check if the pt2 is taking determinants already in the wf
+!      if (is_in_wavefunction(det(N_int,1),N_int)) then
+!        print*, 'A determinant contributing to the pt2 is already in'
+!        print*, 'the wave function:'
+!        call  print_det(det(N_int,1),N_int)
+!        print*,'contribution to the pt2 for the states:', e_pert(:)
+!        print*,'error in the filtering in'
+!        print*, 'cipsi/selection.irp.f sub:  selecte_singles_and_doubles'
+!        print*, 'abort'
+!        call abort
+!      endif
+!!!END_DEBUG

      integer(bit_kind) :: occ(N_int,2), n
      if (h0_type == 'CFG') then
--- a/src/cis/cis.irp.f
+++ b/src/cis/cis.irp.f
@ -74,7 +74,7 @@ subroutine run
    print*,'******************************************************'
    print*,'Excitation energies (au)                     (eV)'
    do i = 2, N_states
-      print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0
+      print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1)) * ha_to_ev
    enddo
    print*,''
  endif
--- a/src/cis_read/EZFIO.cfg
+++ b/src/cis_read/EZFIO.cfg
@ -0,0 +1,7 @@
+[energy]
+type: double precision
+doc: Variational |CIS| energy
+interface: ezfio
+size: (determinants.n_states)
+
+
--- a/src/cis_read/NEED
+++ b/src/cis_read/NEED
@ -0,0 +1,3 @@
+selectors_full
+generators_full
+davidson_undressed
--- a/src/cis_read/README.rst
+++ b/src/cis_read/README.rst
@ -0,0 +1,5 @@
+===
+cis_read
+===
+
+Reads the input WF and performs all singles on top of it.
--- a/src/cis_read/cis_read.irp.f
+++ b/src/cis_read/cis_read.irp.f
@ -0,0 +1,88 @@
+program cis
+  implicit none
+  BEGIN_DOC
+!
+! Configuration Interaction with Single excitations.
+!
+! This program takes a reference Slater determinant of ROHF-like
+! occupancy, and performs all single excitations on top of it.
+! Disregarding spatial symmetry, it computes the `n_states` lowest
+! eigenstates of that CI matrix. (see :option:`determinants n_states`)
+!
+! This program can be useful in many cases:
+!
+!
+! 1. Ground state calculation
+!
+!    To be sure to have the lowest |SCF| solution, perform an :ref:`scf`
+!    (see the :ref:`module_hartree_fock` module), then a :ref:`cis`, save the
+!    natural orbitals (see :ref:`save_natorb`) and re-run an :ref:`scf`
+!    optimization from this |MO| guess.
+!
+!
+! 2. Excited states calculations
+!
+!    The lowest excited states are much likely to be dominated by
+!    single-excitations. Therefore, running a :ref:`cis` will save the
+!    `n_states` lowest states within the |CIS| space in the |EZFIO|
+!    directory, which can afterwards be used as guess wave functions for
+!    a further multi-state |FCI| calculation if :option:`determinants
+!    read_wf` is set to |true| before running the :ref:`fci` executable.
+!
+!
+! If :option:`determinants s2_eig` is set to |true|, the |CIS|
+! will only retain states having the expected |S^2| value (see
+! :option:`determinants expected_s2`). Otherwise, the |CIS| will take
+! the lowest :option:`determinants n_states`, whatever multiplicity
+! they are.
+!
+! .. note::
+!
+!   To discard some orbitals, use the :ref:`qp_set_mo_class` 
+!   command to specify:
+!
+!   * *core* orbitals which will be always doubly occupied
+!
+!   * *act* orbitals where an electron can be either excited from or to
+!
+!   * *del* orbitals which will be never occupied
+!
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run
+end
+
+subroutine run
+  implicit none
+  integer                        :: i
+
+
+  if(pseudo_sym)then
+   call H_apply_cis_sym
+  else
+   call H_apply_cis
+  endif
+  print*,''
+  print *,  'N_det = ', N_det
+  print*,'******************************'
+  print *,  'Energies  of the states:'
+  do i = 1,N_states
+    print *,  i, CI_energy(i)
+  enddo
+  if (N_states > 1) then
+    print*,''
+    print*,'******************************************************'
+    print*,'Excitation energies (au)                     (eV)'
+    do i = 2, N_states
+      print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0
+    enddo
+    print*,''
+  endif
+
+  call ezfio_set_cis_energy(CI_energy)
+  psi_coef = ci_eigenvectors
+  SOFT_TOUCH psi_coef
+  call save_wavefunction_truncated(save_threshold)
+
+end
--- a/src/cis_read/h_apply.irp.f
+++ b/src/cis_read/h_apply.irp.f
@ -0,0 +1,14 @@
+! Generates subroutine H_apply_cis
+! --------------------------------
+
+BEGIN_SHELL [ /usr/bin/env python3 ]
+from generate_h_apply import H_apply
+H = H_apply("cis",do_double_exc=False)
+print(H)
+
+H = H_apply("cis_sym",do_double_exc=False)
+H.filter_only_connected_to_hf()
+print(H)
+
+END_SHELL
+
--- a/src/cisd/30.cisd.bats
+++ b/src/cisd/30.cisd.bats
@ -77,7 +77,7 @@ function run() {
  [[ -n $TRAVIS ]] && skip
  qp set_file ch4.ezfio
  qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]"
-  run -40.2403962667047 -39.8433221754964
+  run -40.2403962667047 -39.843315
 }

@test "SiH3" { # 20.2202s 1.38648m
--- a/src/cisd/cisd.irp.f
+++ b/src/cisd/cisd.irp.f
@ -47,6 +47,37 @@ program cisd
  PROVIDE N_states
  read_wf = .False.
  SOFT_TOUCH read_wf
+
+  integer :: i,k
+
+  if(pseudo_sym)then
+   call H_apply_cisd_sym
+  else
+   call H_apply_cisd
+  endif
+  double precision :: r1, r2
+  double precision, allocatable :: U_csf(:,:)
+
+  allocate(U_csf(N_csf,N_states))
+  U_csf = 0.d0
+  U_csf(1,1) = 1.d0
+  do k=2,N_states
+    do i=1,N_csf
+        call random_number(r1)
+        call random_number(r2)
+        r1 = dsqrt(-2.d0*dlog(r1))
+        r2 = dacos(-1.d0)*2.d0*r2
+        U_csf(i,k) = r1*dcos(r2)
+    enddo
+    U_csf(k,k) = U_csf(k,k) +100.d0
+  enddo
+  do k=1,N_states
+    call normalize(U_csf(1,k),N_csf)
+  enddo
+  call convertWFfromCSFtoDET(N_states,U_csf(1,1),psi_coef(1,1))
+  deallocate(U_csf)
+  SOFT_TOUCH psi_coef
+
  call run
 end

@ -56,20 +87,16 @@ subroutine run
  double precision               :: cisdq(N_states), delta_e
  double precision,external      :: diag_h_mat_elem

-  if(pseudo_sym)then
-   call H_apply_cisd_sym
-  else
-   call H_apply_cisd
-  endif
  psi_coef = ci_eigenvectors
-  SOFT_TOUCH psi_coef
  call save_wavefunction_truncated(save_threshold)
  call ezfio_set_cisd_energy(CI_energy)

  do i = 1,N_states
    k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1)
    delta_E  = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
-    cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
+    if (elec_alpha_num + elec_beta_num >= 4) then
+      cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
+    endif
  enddo
  print *,  'N_det = ', N_det
  print*,''
@ -78,26 +105,43 @@ subroutine run
  do i = 1,N_states
    print *,  i, CI_energy(i)
  enddo
-  print*,''
-  print*,'******************************'
-  print *,  'CISD+Q Energies'
-  do i = 1,N_states
-    print *,  i, cisdq(i)
-  enddo
+  if (elec_alpha_num + elec_beta_num >= 4) then
+    print*,''
+    print*,'******************************'
+    print *,  'CISD+Q Energies'
+    do i = 1,N_states
+      print *,  i, cisdq(i)
+    enddo
+  endif
  if (N_states > 1) then
-    print*,''
-    print*,'******************************'
-    print*,'Excitation energies (au)    (CISD+Q)'
-    do i = 2, N_states
-      print*, i ,CI_energy(i) - CI_energy(1), cisdq(i) - cisdq(1)
-    enddo
-    print*,''
-    print*,'******************************'
-    print*,'Excitation energies (eV)    (CISD+Q)'
-    do i = 2, N_states
-      print*, i ,(CI_energy(i) - CI_energy(1))/0.0367502d0, &
-        (cisdq(i) - cisdq(1)) / 0.0367502d0
-    enddo
+    if (elec_alpha_num + elec_beta_num >= 4) then
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (au)    (CISD+Q)'
+      do i = 2, N_states
+        print*, i ,CI_energy(i) - CI_energy(1), cisdq(i) - cisdq(1)
+      enddo
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (eV)    (CISD+Q)'
+      do i = 2, N_states
+        print*, i ,(CI_energy(i) - CI_energy(1)) * ha_to_ev, &
+          (cisdq(i) - cisdq(1)) * ha_to_ev
+      enddo
+    else
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (au)    (CISD)'
+      do i = 2, N_states
+        print*, i ,CI_energy(i) - CI_energy(1)
+      enddo
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (eV)    (CISD)'
+      do i = 2, N_states
+        print*, i ,(CI_energy(i) - CI_energy(1)) * ha_to_ev
+      enddo
+    endif
  endif

 end
--- a/src/csf/cfgCI_utils.c
+++ b/src/csf/cfgCI_utils.c
@ -68,10 +68,16 @@ void getBFIndexList(int NSOMO, int *BF1, int *IdxListBF1){
                    break;
                }
            }
-            BFcopy[Iidx] = -1;
-            BFcopy[Jidx] = -1;
-            IdxListBF1[Jidx] = Iidx;
-            IdxListBF1[Iidx] = Jidx;
+            if(countN1 <= 0){
+                BFcopy[Iidx] = -1;
+                IdxListBF1[Iidx] = Iidx;
+            }
+            else{
+                BFcopy[Iidx] = -1;
+                BFcopy[Jidx] = -1;
+                IdxListBF1[Jidx] = Iidx;
+                IdxListBF1[Iidx] = Jidx;
+            }
        }
    }

@ -328,10 +334,21 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
                 Get Overlap
    ************************************/
    // Fill matrix
+
+    int rowsbftodetI, colsbftodetI;
+
+    /***********************************
+         Get BFtoDeterminant Matrix
+    ************************************/
+
+    printf("In convertcsftodet\n");
+    convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
+
    int rowsI = 0;
    int colsI = 0;

-    getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+    //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+    getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);


    /***********************************
@ -342,14 +359,6 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl

    gramSchmidt(overlapMatrixI, rowsI, colsI, orthoMatrixI);

-    /***********************************
-         Get BFtoDeterminant Matrix
-    ************************************/
-
-    int rowsbftodetI, colsbftodetI;
-
-    convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
-
    /***********************************
         Get Final CSF to Det Matrix
    ************************************/
@ -1297,12 +1306,16 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
    double fac = 1.0;
    for(int i = 0; i < NSOMO; i++)
        donepq[i] = 0.0;
+    for(int i=0;i<npairs;++i){
+      for(int j=0;j<NSOMO;++j)
+        detslist[i*NSOMO + j]=0;
+    }

    for(int i = 0; i < NSOMO; i++){
        idxp = BF1[i];
        idxq = BF1[idxp];
        // Do one pair only once
-        if(donepq[idxp] > 0.0 || donepq[idxq] > 0.0) continue;
+        if(donepq[idxp] > 0.0 || donepq[idxq] > 0.0 || idxp == idxq) continue;
        fac *= 2.0;
        donepq[idxp] = 1.0;
        donepq[idxq] = 1.0;
@ -1319,15 +1332,19 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
        }
        shft /= 2;
    }
-
+    
    // Now get the addresses
    int inpdet[NSOMO];
    int phase_cfg_to_qp=1;
    int addr = -1;
    for(int i = 0; i < npairs; i++){
-        for(int j = 0; j < NSOMO; j++)
+        for(int j = 0; j < NSOMO; j++) {
            inpdet[j] = detslist[i*NSOMO + j];
+            printf(" %d ",inpdet[j]);
+        }
+        printf("\n");
        findAddofDetDriver(dettree, NSOMO, inpdet, &addr);
+        printf("(%d) - addr  = %d\n",i,addr);
        // Calculate the phase for cfg to QP2 conversion
        //get_phase_cfg_to_qp_inpList(inpdet, NSOMO, &phase_cfg_to_qp);
        //rowvec[addr] = 1.0 * phaselist[i]*phase_cfg_to_qp/sqrt(fac);
@ -1416,7 +1433,6 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
        getIthBFDriver(&bftree, NSOMO, addI, BF1);
        getBFIndexList(NSOMO, BF1, IdxListBF1);

-
        // Get ith row
        getbftodetfunction(&dettree, NSOMO, MS, IdxListBF1, rowvec);

@ -1425,6 +1441,11 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *

        for(int k=0;k<ndets;k++)
            rowvec[k]=0.0;
+
+        for(int j=0;j<NSOMO;++j){
+          BF1[j]=0;
+          IdxListBF1[j]=0;
+        }
    }

    // Garbage collection
@ -1680,7 +1701,6 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in

    gramSchmidt(overlapMatrixJ, rowsJ, colsJ, orthoMatrixJ);

-
    int rowsA = 0;
    int colsA = 0;

--- a/src/csf/tree_utils.c
+++ b/src/csf/tree_utils.c
@ -17,27 +17,31 @@ void buildTree(Tree *bftree,
    if(isomo > NSOMOMax || icpl < 0 || izeros > zeromax ) return;

    // If we find a valid BF assign its address
-    if(isomo == NSOMOMax){
+    if(isomo == NSOMOMax && icpl == MSmax){
        (*inode)->addr = bftree->rootNode->addr;
        bftree->rootNode->addr += 1;
        return;
    }

    // Call 0 branch
-    if(((*inode)->C0) == NULL && izeros+1 <= zeromax){
-        ((*inode)->C0) = malloc(sizeof(Node));
-        (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
-        buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
+    if(izeros+1 <= zeromax){
+      if(((*inode)->C0) == NULL){
+          ((*inode)->C0) = malloc(sizeof(Node));
+          (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
+          buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
+      }
+      else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
    }
-    else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);

    // Call 1 branch
-    if(((*inode)->C1) == NULL && icpl-1 >= 0){
-        ((*inode)->C1) = malloc(sizeof(Node));
-        (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
-        buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
+    if(icpl-1 >=0){
+      if(((*inode)->C1) == NULL){
+          ((*inode)->C1) = malloc(sizeof(Node));
+          (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
+          buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
+      }
+      else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
    }
-    else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);

    return;
 }
--- a/src/davidson/diagonalization_hcsf_dressed.irp.f
+++ b/src/davidson/diagonalization_hcsf_dressed.irp.f
@ -124,7 +124,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
    stop -1
  endif

-  itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1
+  itermax = max(2,min(davidson_sze_max, sze_csf/N_st_diag))+1
  itertot = 0

  if (state_following) then
@ -263,29 +263,20 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
  ! ===================

  converged = .False.
-
+  call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),U_csf(1,1))
  do k=N_st+1,N_st_diag
-    do i=1,sze
+    do i=1,sze_csf
        call random_number(r1)
        call random_number(r2)
        r1 = dsqrt(-2.d0*dlog(r1))
        r2 = dtwo_pi*r2
-        u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
+        U_csf(i,k) = r1*dcos(r2) * u_csf(i,k-N_st)
    enddo
-    u_in(k,k) = u_in(k,k) + 10.d0
+    U_csf(k,k) = u_csf(k,k) + 10.d0
  enddo
  do k=1,N_st_diag
-    call normalize(u_in(1,k),sze)
+    call normalize(U_csf(1,k),sze_csf)
  enddo
-
-  do k=1,N_st_diag
-    do i=1,sze
-      U(i,k) = u_in(i,k)
-    enddo
-  enddo
-
-  ! Make random verctors eigenstates of S2
-  call convertWFfromDETtoCSF(N_st_diag,U(1,1),U_csf(1,1))
  call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1))

  do while (.not.converged)
--- a/src/davidson/diagonalize_ci.irp.f
+++ b/src/davidson/diagonalize_ci.irp.f
@ -1,3 +1,13 @@
+BEGIN_PROVIDER [ character*(3), sigma_vector_algorithm ]
+ implicit none
+ BEGIN_DOC
+ ! If 'det', use <Psi_det|H|Psi_det> in Davidson
+ !
+ ! If 'cfg', use <Psi_csf|H|Psi_csf> in Davidson
+ END_DOC
+ sigma_vector_algorithm = 'det'
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
  implicit none
  BEGIN_DOC
@ -61,9 +71,18 @@ END_PROVIDER
   if (diag_algorithm == "Davidson") then

     if (do_csf) then
-       call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
-         size(CI_eigenvectors,1),CI_electronic_energy,               &
-         N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+       if (sigma_vector_algorithm == 'det') then
+         call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
+           size(CI_eigenvectors,1),CI_electronic_energy,               &
+           N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+!       else if (sigma_vector_algorithm == 'cfg') then
+!       call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
+!         size(CI_eigenvectors,1),CI_electronic_energy,               &
+!         N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+!       else
+!         print *, irp_here
+!         stop 'bug'
+       endif
     else
       call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_s2, &
         size(CI_eigenvectors,1),CI_electronic_energy,               &
--- a/src/determinants/determinants.irp.f
+++ b/src/determinants/determinants.irp.f
@ -77,28 +77,31 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
  END_DOC
  PROVIDE ezfio_filename
  logical                        :: exists
-  if (mpi_master) then
-    call ezfio_has_determinants_n_det(exists)
-    if (exists) then
-      call ezfio_get_determinants_n_det(psi_det_size)
-    else
-      psi_det_size = 1
+  psi_det_size = N_states
+  PROVIDE mpi_master
+  if (read_wf) then
+    if (mpi_master) then
+      call ezfio_has_determinants_n_det(exists)
+      if (exists) then
+        call ezfio_get_determinants_n_det(psi_det_size)
+      else
+        psi_det_size = N_states
+      endif
+      call write_int(6,psi_det_size,'Dimension of the psi arrays')
    endif
-    call write_int(6,psi_det_size,'Dimension of the psi arrays')
+    IRP_IF MPI_DEBUG
+      print *,  irp_here, mpi_rank
+      call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+    IRP_ENDIF
+    IRP_IF MPI
+      include 'mpif.h'
+      integer                        :: ierr
+      call MPI_BCAST( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read psi_det_size with MPI'
+      endif
+    IRP_ENDIF
  endif
-  IRP_IF MPI_DEBUG
-    print *,  irp_here, mpi_rank
-    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
-  IRP_ENDIF
-  IRP_IF MPI
-    include 'mpif.h'
-    integer                        :: ierr
-    call MPI_BCAST( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    if (ierr /= MPI_SUCCESS) then
-      stop 'Unable to read psi_det_size with MPI'
-    endif
-  IRP_ENDIF
-

 END_PROVIDER

@ -539,7 +542,7 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
  integer                        :: i,j,k, ndet_qp_edit

  if (mpi_master) then
-    ndet_qp_edit = min(ndet,N_det_qp_edit)
+    ndet_qp_edit = min(ndet,10000)

    call ezfio_set_determinants_N_int(N_int)
    call ezfio_set_determinants_bit_kind(bit_kind)
--- a/src/determinants/dipole_moments.irp.f
+++ b/src/determinants/dipole_moments.irp.f
@ -66,9 +66,13 @@ END_PROVIDER
    write(*,'(i16)',advance='no') i
  end do
  write(*,*) ''
-  write(*,'(A17,100(1pE16.8))') 'x_dipole_moment = ',x_dipole_moment
-  write(*,'(A17,100(1pE16.8))') 'y_dipole_moment = ',y_dipole_moment
-  write(*,'(A17,100(1pE16.8))') 'z_dipole_moment = ',z_dipole_moment
+  write(*,'(A23,100(1pE16.8))') 'x_dipole_moment (au) = ',x_dipole_moment
+  write(*,'(A23,100(1pE16.8))') 'y_dipole_moment (au) = ',y_dipole_moment
+  write(*,'(A23,100(1pE16.8))') 'z_dipole_moment (au) = ',z_dipole_moment
+  write(*,*) ''
+  write(*,'(A23,100(1pE16.8))') 'x_dipole_moment (D) = ',x_dipole_moment * au_to_D
+  write(*,'(A23,100(1pE16.8))') 'y_dipole_moment (D) = ',y_dipole_moment * au_to_D
+  write(*,'(A23,100(1pE16.8))') 'z_dipole_moment (D) = ',z_dipole_moment * au_to_D
  !print*, 'x_dipole_moment                  = ',x_dipole_moment
  !print*, 'y_dipole_moment                  = ',y_dipole_moment
  !print*, 'z_dipole_moment                  = ',z_dipole_moment
--- a/src/determinants/slater_rules.irp.f
+++ b/src/determinants/slater_rules.irp.f
@ -623,7 +623,8 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
  integer                        :: occ(Nint*bit_kind_size,2)
  double precision               :: diag_H_mat_elem, phase
  integer                        :: n_occ_ab(2)
-  PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
+  PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
+  PROVIDE ao_one_e_integrals mo_one_e_integrals

  ASSERT (Nint > 0)
  ASSERT (Nint == N_int)
@ -681,7 +682,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
    case (1)
      call get_single_excitation(key_i,key_j,exc,phase,Nint)
      !DIR$ FORCEINLINE
-      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
      if (exc(0,1,1) == 1) then
        ! Single alpha
        m = exc(1,1,1)
@ -700,10 +700,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
  end select
 end

-
-
-
-
 subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble,phase)
  use bitmasks
  implicit none
@ -1038,7 +1034,6 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)

 end

-
 subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
  use bitmasks
  implicit none
--- a/src/determinants/slater_rules_wee_mono.irp.f
+++ b/src/determinants/slater_rules_wee_mono.irp.f
@ -282,9 +282,7 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij)
  double precision               :: get_two_e_integral
  integer                        :: m,n,p,q
  integer                        :: i,j,k
-  integer                        :: occ(Nint*bit_kind_size,2)
  double precision               :: diag_H_mat_elem, phase,phase_2
-  integer                        :: n_occ_ab(2)
  PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals ref_bitmask_two_e_energy

  ASSERT (Nint > 0)
@ -342,7 +340,6 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij)
    case (1)
      call get_single_excitation(key_i,key_j,exc,phase,Nint)
      !DIR$ FORCEINLINE
-      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
      if (exc(0,1,1) == 1) then
        ! Mono alpha
        m = exc(1,1,1)
--- a/src/fci/40.fci.bats
+++ b/src/fci/40.fci.bats
@ -8,12 +8,12 @@ function run() {
  test_exe fci || skip
  qp edit --check
  qp set perturbation do_pt2 False
-  qp set determinants n_det_max 8000
+  qp set determinants n_det_max $3
  qp set determinants n_states  1
  qp set davidson threshold_davidson 1.e-10
  qp set davidson n_states_diag 8
  qp run fci
-  energy1="$(ezfio get fci energy | tr '[]' ' ' | cut -d ',' -f 1)"
+  energy1="$(qp get fci energy | tr '[]' ' ' | cut -d ',' -f 1)"
  eq $energy1 $1 $thresh
 }

@ -22,166 +22,167 @@ function run_stoch() {
  thresh=$2
  test_exe fci || skip
  qp set perturbation do_pt2 True
+  qp set perturbation pt2_relative_error 0.005
  qp set determinants n_det_max $3
  qp set determinants n_states  1
  qp set davidson threshold_davidson 1.e-10
  qp set davidson n_states_diag 1
  qp run fci
-  energy1="$(ezfio get fci energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)"
+  energy1="$(qp get fci energy_extrapolated | tr '[]' ' ' | cut -d ',' -f 1)"
  eq $energy1 $1 $thresh
 }

-@test "B-B" {
+
+@test "B-B" {  # 0:00:10
  qp set_file b2_stretched.ezfio
  qp set determinants n_det_max 10000
  qp set_frozen_core
-  run_stoch -49.14103054419 3.e-4 10000
+  run_stoch -49.14097596  0.0001  10000
 }

-@test "F2" { # 4.07m
-  [[ -n $TRAVIS ]] && skip
-  qp set_file f2.ezfio
-  qp set_frozen_core
-  run_stoch -199.304922384814 3.e-4  100000
-}
-
-@test "NH3" { # 10.6657s
+@test "NH3" { # 0:00:11
  qp set_file nh3.ezfio
  qp set_mo_class --core="[1-4]" --act="[5-72]"
-  run -56.244753429144986  3.e-4  100000
+  run -56.24474908 1.e-5  10000
 }

-@test "DHNO" { # 11.4721s
+@test "DHNO" { # 0:00:10 
  qp set_file dhno.ezfio
  qp set_mo_class --core="[1-7]" --act="[8-64]"
-  run -130.459020029816 3.e-4  100000
+  run -130.45904647  1.e-4  10000
 }

-@test "HCO" { # 12.2868s
+@test "HCO" { # 0:01:16 
  qp set_file hco.ezfio
-  run -113.389297812482 6.e-4  100000
+  run_stoch -113.41448940  2.e-3  50000
 }

-@test "H2O2" { # 12.9214s
+@test "H2O2" { # 0:01:48 
  qp set_file h2o2.ezfio
  qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]"
-  run -151.00467 1.e-4  100000
+  run_stoch -151.02437936 2.e-3  100000
 }

-@test "HBO" { # 13.3144s
+@test "HBO" { # 0:00:46 
  [[ -n $TRAVIS ]] && skip
  qp set_file hbo.ezfio
-  run -100.212560384678 1.e-3  100000
+  run_stoch -100.221198108988 2.e-3  50000
 }

-@test "H2O" { # 11.3727s
+@test "H2O" { # 0:01:05 
  [[ -n $TRAVIS ]] && skip
  qp set_file h2o.ezfio
-  run -76.2361605151999 3.e-4  100000
+  run_stoch -76.241332121813  1.e-3  100000
 }

-@test "ClO" { # 13.3755s
+@test "ClO" { # 0:03:07 
  [[ -n $TRAVIS ]] && skip
  qp set_file clo.ezfio
-  run -534.545616787223 3.e-4  100000
+  run_stoch -534.573564655419 1.e-3  100000
 }

-@test "SO" { # 13.4952s
+@test "SO" { # 0:01:49 
  [[ -n $TRAVIS ]] && skip
  qp set_file so.ezfio
-  run -26.0096209515081 1.e-3  100000
+  run_stoch -26.04335528  5.e-3  100000
 }

-@test "H2S" { # 13.6745s
+@test "H2S" { # 0:01:12
  [[ -n $TRAVIS ]] && skip
  qp set_file h2s.ezfio
-  run -398.859168655255 3.e-4  100000
+  run_stoch -398.865173546866 1.e-3  50000
 }

-@test "OH" { # 13.865s
+@test "OH" { # 0:00:41
  [[ -n $TRAVIS ]] && skip
  qp set_file oh.ezfio
-  run -75.6121856748294 3.e-4   100000
+  run_stoch -75.6193013819546 1.e-3  50000
 }

-@test "SiH2_3B1" { # 13.938ss
+@test "SiH2_3B1" { # 0:00:50
  [[ -n $TRAVIS ]] && skip
  qp set_file sih2_3b1.ezfio
-  run -290.0175411299477 3.e-4  100000
+  run_stoch -290.01754869  3.e-5  50000
 }

-@test "H3COH" { # 14.7299s
+@test "H3COH" { # 0:01:05
  [[ -n $TRAVIS ]] && skip
  qp set_file h3coh.ezfio
-  run -115.205191406072 3.e-4  100000
+  run_stoch -115.224147057725 2.e-3  50000
 }

-@test "SiH3" { # 15.99s
+@test "SiH3" { # 0:01:09
  [[ -n $TRAVIS ]] && skip
  qp set_file sih3.ezfio
-  run -5.57241217753818 3.e-4   100000
+  run_stoch -5.57812512359276 1.e-3  50000
 }

-@test "CH4" { # 16.1612s
+@test "CH4" { # 0:02:06
  [[ -n $TRAVIS ]] && skip
  qp set_file ch4.ezfio
  qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]"
-  run -40.2409678239136 3.e-4  100000
+  run_stoch -40.2419474611994 1.e-4  100000
 }

-@test "ClF" { # 16.8864s
+@test "ClF" { # 0:01:55 
  [[ -n $TRAVIS ]] && skip
  qp set_file clf.ezfio
-  run -559.169313755572 3.e-4  100000
+  run_stoch -559.20666465 1.e-2  50000
 }

-@test "SO2" { # 17.5645s
+@test "SO2" { # 0:00:24
  [[ -n $TRAVIS ]] && skip
  qp set_file so2.ezfio
  qp set_mo_class --core="[1-8]" --act="[9-87]"
-  run -41.5746738713298 3.e-4  100000
+  run_stoch -41.57468756 1.e-4  50000
 }

-@test "C2H2" { # 17.6827s
+@test "C2H2" { # 0:00:57 
  [[ -n $TRAVIS ]] && skip
  qp set_file c2h2.ezfio
  qp set_mo_class --act="[1-30]" --del="[31-36]"
-  run -12.3685464085969 3.e-4  100000
+  run_stoch -12.3862664765532 1.e-3  50000
 }

-@test "N2" { # 18.0198s
+@test "N2" { # 0:01:15
  [[ -n $TRAVIS ]] && skip
  qp set_file n2.ezfio
  qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]"
-  run -109.28681540699360 3.e-4  100000
+  run_stoch -109.311954243348 2.e-3  50000
 }

-@test "N2H4" { # 18.5006s
+@test "N2H4" { # 0:00:51
  [[ -n $TRAVIS ]] && skip
  qp set_file n2h4.ezfio
  qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]"
-  run -111.367332681559 3.e-4  100000
+  run_stoch -111.38119165053 1.e-3 50000
 }

-@test "CO2" { # 21.1748s
+@test "CO2" { # 0:01:00
  [[ -n $TRAVIS ]] && skip
  qp set_file co2.ezfio
  qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]"
-  run -187.968547952413 3.e-4  100000
+  run_stoch -188.002190327443 2.e-3  50000
 }

-
-@test "[Cu(NH3)4]2+" { # 25.0417s
+@test "[Cu(NH3)4]2+" { # 0:01:53
  [[ -n $TRAVIS ]] && skip
  qp set_file cu_nh3_4_2plus.ezfio
  qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
-  run -1862.9869374387192  3.e-04  100000
+  run_stoch -1862.98705340328 1.e-05  50000
 }

-@test "HCN" { # 20.3273s
+@test "HCN" { # 0:01:26
  [[ -n $TRAVIS ]] && skip
  qp set_file hcn.ezfio
  qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]"
-  run -93.0771143355433 3.e-4  100000
+  run_stoch -93.0980746734051 5.e-4  50000
+}
+
+@test "F2" { # 0:03:34
+  [[ -n $TRAVIS ]] && skip
+  qp set_file f2.ezfio
+  qp set_frozen_core
+  run_stoch -199.307512211742 0.002  100000
 }

--- a/src/fci/EZFIO.cfg
+++ b/src/fci/EZFIO.cfg
@ -10,3 +10,9 @@ doc: Calculated |FCI| energy + |PT2|
 interface: ezfio
 size: (determinants.n_states)

+[energy_extrapolated]
+type: double precision
+doc: Calculated |FCI| energy + |PT2|
+interface: ezfio
+size: (determinants.n_states)
+
--- a/src/iterations/print_extrapolation.irp.f
+++ b/src/iterations/print_extrapolation.irp.f
@ -35,12 +35,13 @@ subroutine print_extrapolated_energy
    do k=2,min(N_iter,8)
      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
          extrapolated_energy(k,i) - extrapolated_energy(k,1), &
-          (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
+          (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * ha_to_ev
    enddo
    write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
  enddo

  print *,  ''
+  call ezfio_set_fci_energy_extrapolated(extrapolated_energy(min(N_iter,3),1:N_states))

 end subroutine

--- a/src/iterations/print_summary.irp.f
+++ b/src/iterations/print_summary.irp.f
@ -36,7 +36,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
  write(*,fmt) '# E          ', e_(1:N_states_p)
  if (N_states_p > 1) then
    write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
-    write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
+    write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*ha_to_ev
  endif
  write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
  write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p)
@ -47,8 +47,8 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
  if (N_states_p > 1) then
    write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), &
      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p)
-    write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, &
-      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p)
+    write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*ha_to_ev, &
+      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*ha_to_ev, k=1,N_states_p)
  endif
  write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
  write(*,fmt)
@ -82,19 +82,19 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
    print *, 'Variational Energy difference (au | eV)'
    do i=2, N_states_p
      print*,'Delta E = ', (e_(i) - e_(1)), &
-        (e_(i) - e_(1)) * 27.211396641308d0
+        (e_(i) - e_(1)) * ha_to_ev
    enddo
    print *,  '-----'
    print*, 'Variational + perturbative Energy difference (au | eV)'
    do i=2, N_states_p
      print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), &
-        (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0
+        (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * ha_to_ev
    enddo
    print *,  '-----'
    print*, 'Variational + renormalized perturbative Energy difference (au | eV)'
    do i=2, N_states_p
      print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), &
-        (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0
+        (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * ha_to_ev
    enddo
  endif

--- a/src/mo_two_e_erf_ints/map_integrals_erf.irp.f
+++ b/src/mo_two_e_erf_ints/map_integrals_erf.irp.f
@ -235,11 +235,11 @@ subroutine get_mo_two_e_integrals_erf_ij(k,l,sze,out_array,map)

  logical :: integral_is_in_map
  if (key_kind == 8) then
-    call i8radix_sort(hash,iorder,kk,-1)
+    call i8sort(hash,iorder,kk)
  else if (key_kind == 4) then
-    call iradix_sort(hash,iorder,kk,-1)
+    call isort(hash,iorder,kk)
  else if (key_kind == 2) then
-    call i2radix_sort(hash,iorder,kk,-1)
+    call i2sort(hash,iorder,kk)
  endif

  call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk)
@ -290,11 +290,11 @@ subroutine get_mo_two_e_integrals_erf_i1j1(k,l,sze,out_array,map)

  logical :: integral_is_in_map
  if (key_kind == 8) then
-    call i8radix_sort(hash,iorder,kk,-1)
+    call i8sort(hash,iorder,kk)
  else if (key_kind == 4) then
-    call iradix_sort(hash,iorder,kk,-1)
+    call isort(hash,iorder,kk)
  else if (key_kind == 2) then
-    call i2radix_sort(hash,iorder,kk,-1)
+    call i2sort(hash,iorder,kk)
  endif

  call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk)
--- a/src/mo_two_e_ints/mo_bi_integrals.irp.f
+++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f
@ -53,7 +53,11 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
 !    call four_idx_novvvv
    call four_idx_novvvv_old
  else
-    call add_integrals_to_map(full_ijkl_bitmask_4)
+    if (32.d-9*dble(ao_num)**4 < dble(qp_max_mem)) then
+      call four_idx_dgemm
+    else
+      call add_integrals_to_map(full_ijkl_bitmask_4)
+    endif
  endif

  call wall_time(wall_2)
@ -77,6 +81,94 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]

 END_PROVIDER

+subroutine four_idx_dgemm
+  implicit none
+  integer :: p,q,r,s,i,j,k,l
+  double precision, allocatable :: a1(:,:,:,:)
+  double precision, allocatable :: a2(:,:,:,:)
+
+  allocate (a1(ao_num,ao_num,ao_num,ao_num))
+
+  print *, 'Getting AOs'
+  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,r,s)
+  do s=1,ao_num
+    do r=1,ao_num
+      do q=1,ao_num
+        call get_ao_two_e_integrals(q,r,s,ao_num,a1(1,q,r,s))
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+  print *, '1st transformation'
+  ! 1st transformation
+  allocate (a2(ao_num,ao_num,ao_num,mo_num))
+  call dgemm('T','N', (ao_num*ao_num*ao_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*ao_num*ao_num))
+
+  ! 2nd transformation
+  print *, '2nd transformation'
+  deallocate (a1)
+  allocate (a1(ao_num,ao_num,mo_num,mo_num))
+  call dgemm('T','N', (ao_num*ao_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (ao_num*ao_num*mo_num))
+
+  ! 3rd transformation
+  print *, '3rd transformation'
+  deallocate (a2)
+  allocate (a2(ao_num,mo_num,mo_num,mo_num))
+  call dgemm('T','N', (ao_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*mo_num*mo_num))
+
+  ! 4th transformation
+  print *, '4th transformation'
+  deallocate (a1)
+  allocate (a1(mo_num,mo_num,mo_num,mo_num))
+  call dgemm('T','N', (mo_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (mo_num*mo_num*mo_num))
+
+  deallocate (a2)
+
+  integer :: n_integrals, size_buffer
+  integer(key_kind)  , allocatable :: buffer_i(:)
+  real(integral_kind), allocatable :: buffer_value(:)
+  size_buffer = min(ao_num*ao_num*ao_num,16000000)
+
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals)
+  allocate ( buffer_i(size_buffer), buffer_value(size_buffer) )
+
+  n_integrals = 0
+  !$OMP DO
+  do l=1,mo_num
+    do k=1,mo_num
+      do j=1,l
+        do i=1,k
+            if (abs(a1(i,j,k,l)) < mo_integrals_threshold) then
+              cycle
+            endif
+            n_integrals += 1
+            buffer_value(n_integrals) = a1(i,j,k,l)
+            !DIR$ FORCEINLINE
+            call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals))
+            if (n_integrals == size_buffer) then
+              call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
+              n_integrals = 0
+            endif
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+
+  call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
+
+  deallocate(buffer_i, buffer_value)
+  !$OMP END PARALLEL
+
+  deallocate (a1)
+
+  call map_unique(mo_integrals_map)
+
+  integer*8                      :: get_mo_map_size, mo_map_size
+  mo_map_size = get_mo_map_size()
+
+end subroutine

 subroutine add_integrals_to_map(mask_ijkl)
  use bitmasks
--- a/src/tools/truncate_wf.irp.f
+++ b/src/tools/truncate_wf.irp.f
@ -54,11 +54,23 @@ subroutine routine_s2
 double precision, allocatable :: psi_coef_tmp(:,:)
 integer :: i,j,k
 double precision :: accu(N_states)
+ integer :: weights(0:16), ix
+ double precision :: x

- print *,  'Weights of the CFG'
+ weights(:) = 0
 do i=1,N_det
-   print *,  i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
+   x = -dlog(1.d-32+sum(weight_configuration(det_to_configuration(i),:)))/dlog(10.d0)
+   ix = min(int(x), 16)
+   weights(ix) += 1
 enddo
+
+ print *,  'Histogram of the weights of the CFG'
+ do i=0,15
+   print *, '  10^{-', i, '}   ', weights(i) 
+ end do
+ print *, '< 10^{-', 15, '}   ', weights(16) 
+
+
 print*, 'Min weight of the configuration?'
 read(5,*) wmin

--- a/src/utils/format_w_error.irp.f
+++ b/src/utils/format_w_error.irp.f
@ -0,0 +1,71 @@
+subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error)
+
+  implicit none
+ 
+  BEGIN_DOC
+  ! Format for double precision, value(error)
+  END_DOC
+
+  ! in
+  ! | value         | double precision | value...      |
+  ! | error         | double precision | error...      |
+  ! | size_nb       | integer          | X in FX.Y     |
+  ! | max_nb_digits | integer          | Max Y in FX.Y |
+
+  ! out
+  ! | format_value  | character | string FX.Y for the format |
+  ! | str_error     | character | string of the error |
+
+  ! internal
+  ! | str_size      | character | size in string format                             |
+  ! | nb_digits     | integer   | number of digits Y in FX.Y depending of the error |
+  ! | str_nb_digits | character | nb_digits in string format                        |
+  ! | str_exp       | character | string of the value in exponential format         |
+
+  ! in
+  double precision, intent(in)   :: error, value
+  integer, intent(in)            :: size_nb, max_nb_digits
+
+  ! out
+  character(len=20), intent(out) :: str_error, format_value
+
+  ! internal
+  character(len=20)              :: str_size, str_nb_digits, str_exp
+  integer                        :: nb_digits
+
+  ! max_nb_digit: Y max
+  ! size_nb = Size of the double: X (FX.Y)
+  write(str_size,'(I3)') size_nb
+
+  ! Error
+  write(str_exp,'(1pE20.0)') error
+  str_error = trim(adjustl(str_exp))
+  
+  ! Number of digit: Y (FX.Y) from the exponent
+  str_nb_digits = str_exp(19:20)
+  read(str_nb_digits,*) nb_digits
+ 
+  ! If the error is 0d0
+  if (error <= 1d-16) then 
+    write(str_nb_digits,*) max_nb_digits
+  endif
+
+  ! If the error is too small 
+  if (nb_digits > max_nb_digits) then
+      write(str_nb_digits,*) max_nb_digits
+      str_error(1:1) = '0'
+  endif
+
+  ! If the error is too big (>= 0.5)
+  if (error >= 0.5d0) then
+    str_nb_digits = '1'
+    str_error(1:1) = '*'
+  endif
+
+  ! FX.Y,A1,A1,A1 for value(str_error)
+  !string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1'
+
+  ! FX.Y just for the value 
+  format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))
+
+end
--- a/src/utils/map_module.f90
+++ b/src/utils/map_module.f90
@ -238,11 +238,11 @@ subroutine cache_map_sort(map)
      iorder(i) = i
    enddo
    if (cache_key_kind == 2) then
-      call i2radix_sort(map%key,iorder,map%n_elements,-1)
+      call i2sort(map%key,iorder,map%n_elements,-1)
    else if (cache_key_kind == 4) then
-      call iradix_sort(map%key,iorder,map%n_elements,-1)
+      call isort(map%key,iorder,map%n_elements,-1)
    else if (cache_key_kind == 8) then
-      call i8radix_sort(map%key,iorder,map%n_elements,-1)
+      call i8sort(map%key,iorder,map%n_elements,-1)
    endif
    if (integral_kind == 4) then
      call set_order(map%value,iorder,map%n_elements)
--- a/src/utils/qsort.c
+++ b/src/utils/qsort.c
@ -0,0 +1,373 @@
+/* [[file:~/qp2/src/utils/qsort.org::*Generated%20C%20file][Generated C file:1]] */
+#include <stdlib.h>
+#include <stdint.h>
+
+struct int16_t_comp {
+  int16_t    x;
+  int32_t i;
+};
+
+int compare_int16_t( const void * l, const void * r )
+{
+  const int16_t * restrict _l= l;
+  const int16_t * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_int16_t(int16_t* restrict A_in, int32_t* restrict iorder, int32_t isize) {
+  struct int16_t_comp* A = malloc(isize * sizeof(struct int16_t_comp));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct int16_t_comp), compare_int16_t);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_int16_t_noidx(int16_t* A, int32_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(int16_t), compare_int16_t);
+}
+
+
+struct int16_t_comp_big {
+  int16_t    x;
+  int64_t i;
+};
+
+int compare_int16_t_big( const void * l, const void * r )
+{
+  const int16_t * restrict _l= l;
+  const int16_t * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_int16_t_big(int16_t* restrict A_in, int64_t* restrict iorder, int64_t isize) {
+  struct int16_t_comp_big* A = malloc(isize * sizeof(struct int16_t_comp_big));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct int16_t_comp_big), compare_int16_t_big);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_int16_t_noidx_big(int16_t* A, int64_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(int16_t), compare_int16_t_big);
+}
+
+
+struct int32_t_comp {
+  int32_t    x;
+  int32_t i;
+};
+
+int compare_int32_t( const void * l, const void * r )
+{
+  const int32_t * restrict _l= l;
+  const int32_t * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_int32_t(int32_t* restrict A_in, int32_t* restrict iorder, int32_t isize) {
+  struct int32_t_comp* A = malloc(isize * sizeof(struct int32_t_comp));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct int32_t_comp), compare_int32_t);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_int32_t_noidx(int32_t* A, int32_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(int32_t), compare_int32_t);
+}
+
+
+struct int32_t_comp_big {
+  int32_t    x;
+  int64_t i;
+};
+
+int compare_int32_t_big( const void * l, const void * r )
+{
+  const int32_t * restrict _l= l;
+  const int32_t * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_int32_t_big(int32_t* restrict A_in, int64_t* restrict iorder, int64_t isize) {
+  struct int32_t_comp_big* A = malloc(isize * sizeof(struct int32_t_comp_big));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct int32_t_comp_big), compare_int32_t_big);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_int32_t_noidx_big(int32_t* A, int64_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(int32_t), compare_int32_t_big);
+}
+
+
+struct int64_t_comp {
+  int64_t    x;
+  int32_t i;
+};
+
+int compare_int64_t( const void * l, const void * r )
+{
+  const int64_t * restrict _l= l;
+  const int64_t * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_int64_t(int64_t* restrict A_in, int32_t* restrict iorder, int32_t isize) {
+  struct int64_t_comp* A = malloc(isize * sizeof(struct int64_t_comp));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct int64_t_comp), compare_int64_t);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_int64_t_noidx(int64_t* A, int32_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(int64_t), compare_int64_t);
+}
+
+
+struct int64_t_comp_big {
+  int64_t    x;
+  int64_t i;
+};
+
+int compare_int64_t_big( const void * l, const void * r )
+{
+  const int64_t * restrict _l= l;
+  const int64_t * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_int64_t_big(int64_t* restrict A_in, int64_t* restrict iorder, int64_t isize) {
+  struct int64_t_comp_big* A = malloc(isize * sizeof(struct int64_t_comp_big));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct int64_t_comp_big), compare_int64_t_big);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_int64_t_noidx_big(int64_t* A, int64_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(int64_t), compare_int64_t_big);
+}
+
+
+struct double_comp {
+  double    x;
+  int32_t i;
+};
+
+int compare_double( const void * l, const void * r )
+{
+  const double * restrict _l= l;
+  const double * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_double(double* restrict A_in, int32_t* restrict iorder, int32_t isize) {
+  struct double_comp* A = malloc(isize * sizeof(struct double_comp));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct double_comp), compare_double);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_double_noidx(double* A, int32_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(double), compare_double);
+}
+
+
+struct double_comp_big {
+  double    x;
+  int64_t i;
+};
+
+int compare_double_big( const void * l, const void * r )
+{
+  const double * restrict _l= l;
+  const double * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_double_big(double* restrict A_in, int64_t* restrict iorder, int64_t isize) {
+  struct double_comp_big* A = malloc(isize * sizeof(struct double_comp_big));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct double_comp_big), compare_double_big);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_double_noidx_big(double* A, int64_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(double), compare_double_big);
+}
+
+
+struct float_comp {
+  float    x;
+  int32_t i;
+};
+
+int compare_float( const void * l, const void * r )
+{
+  const float * restrict _l= l;
+  const float * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_float(float* restrict A_in, int32_t* restrict iorder, int32_t isize) {
+  struct float_comp* A = malloc(isize * sizeof(struct float_comp));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct float_comp), compare_float);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_float_noidx(float* A, int32_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(float), compare_float);
+}
+
+
+struct float_comp_big {
+  float    x;
+  int64_t i;
+};
+
+int compare_float_big( const void * l, const void * r )
+{
+  const float * restrict _l= l;
+  const float * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_float_big(float* restrict A_in, int64_t* restrict iorder, int64_t isize) {
+  struct float_comp_big* A = malloc(isize * sizeof(struct float_comp_big));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct float_comp_big), compare_float_big);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_float_noidx_big(float* A, int64_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(float), compare_float_big);
+}
+/* Generated C file:1 ends here */
--- a/src/utils/qsort.org
+++ b/src/utils/qsort.org
@ -0,0 +1,169 @@
+#+TITLE: Quick sort binding for Fortran
+
+* C template
+
+#+NAME: c_template
+#+BEGIN_SRC c
+struct TYPE_comp_big {
+  TYPE    x;
+  int32_t i;
+};
+
+int compare_TYPE_big( const void * l, const void * r )
+{
+  const TYPE * restrict _l= l;
+  const TYPE * restrict _r= r;
+  if( *_l > *_r ) return 1;
+  if( *_l < *_r ) return -1;
+  return 0;
+}
+
+void qsort_TYPE_big(TYPE* restrict A_in, int32_t* restrict iorder, int32_t isize) {
+  struct TYPE_comp_big* A = malloc(isize * sizeof(struct TYPE_comp_big));
+  if (A == NULL) return;
+
+  for (int i=0 ; i<isize ; ++i) {
+    A[i].x = A_in[i];
+    A[i].i = iorder[i];
+  }
+
+  qsort( (void*) A, (size_t) isize, sizeof(struct TYPE_comp_big), compare_TYPE_big);
+
+  for (int i=0 ; i<isize ; ++i) {
+    A_in[i] = A[i].x;
+    iorder[i] = A[i].i;
+  }
+  free(A);
+}
+
+void qsort_TYPE_noidx_big(TYPE* A, int32_t isize) {
+  qsort( (void*) A, (size_t) isize, sizeof(TYPE), compare_TYPE_big);
+}
+#+END_SRC
+
+* Fortran template
+
+#+NAME:f_template
+#+BEGIN_SRC f90
+subroutine Lsort_big_c(A, iorder, isize) bind(C, name="qsort_TYPE_big")
+  use iso_c_binding
+  integer(c_int32_t), value :: isize
+  integer(c_int32_t)        :: iorder(isize)
+  real   (c_TYPE)         :: A(isize)
+end subroutine Lsort_big_c
+
+subroutine Lsort_noidx_big_c(A, isize) bind(C, name="qsort_TYPE_noidx_big")
+  use iso_c_binding
+  integer(c_int32_t), value :: isize
+  real   (c_TYPE)         :: A(isize)
+end subroutine Lsort_noidx_big_c
+
+#+END_SRC
+
+#+NAME:f_template2
+#+BEGIN_SRC f90
+subroutine Lsort_big(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int32_t)        :: isize
+  integer(c_int32_t)        :: iorder(isize)
+  real   (c_TYPE)         :: A(isize)
+  call Lsort_big_c(A, iorder, isize)
+end subroutine Lsort_big
+
+subroutine Lsort_noidx_big(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int32_t) :: isize
+  real   (c_TYPE)    :: A(isize)
+  call Lsort_noidx_big_c(A, isize)
+end subroutine Lsort_noidx_big
+
+#+END_SRC
+
+* Python scripts for type replacements
+
+#+NAME: replaced
+#+begin_src python :results output :noweb yes
+data = """
+<<c_template>>
+"""
+for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]:
+    print( data.replace("TYPE", typ).replace("_big", "") )
+    print( data.replace("int32_t", "int64_t").replace("TYPE", typ) )
+#+end_src
+
+#+NAME: replaced_f
+#+begin_src python :results output :noweb yes
+data = """
+<<f_template>>
+"""
+c1 = {
+    "int16_t": "i2",
+    "int32_t": "i",
+    "int64_t": "i8",
+    "double": "d",
+    "float": ""
+}
+c2 = {
+    "int16_t": "integer",
+    "int32_t": "integer",
+    "int64_t": "integer",
+    "double": "real",
+    "float": "real"
+}
+
+for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]:
+    print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("TYPE", typ).replace("_big", "") )
+    print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("int32_t", "int64_t").replace("TYPE", typ) )
+#+end_src
+
+#+NAME: replaced_f2
+#+begin_src python :results output :noweb yes
+data = """
+<<f_template2>>
+"""
+c1 = {
+    "int16_t": "i2",
+    "int32_t": "i",
+    "int64_t": "i8",
+    "double": "d",
+    "float": ""
+}
+c2 = {
+    "int16_t": "integer",
+    "int32_t": "integer",
+    "int64_t": "integer",
+    "double": "real",
+    "float": "real"
+}
+
+for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]:
+    print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("TYPE", typ).replace("_big", "") )
+    print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("int32_t", "int64_t").replace("TYPE", typ) )
+#+end_src
+
+* Generated C file
+  
+#+BEGIN_SRC c :comments link :tangle qsort.c :noweb yes
+#include <stdlib.h>
+#include <stdint.h>
+<<replaced()>>
+#+END_SRC
+
+* Generated Fortran file
+
+#+BEGIN_SRC f90 :tangle qsort_module.f90 :noweb yes
+module qsort_module
+  use iso_c_binding
+  
+  interface
+     <<replaced_f()>>
+  end interface
+
+end module qsort_module
+
+<<replaced_f2()>>
+
+#+END_SRC
+
--- a/src/utils/qsort_module.f90
+++ b/src/utils/qsort_module.f90
@ -0,0 +1,347 @@
+module qsort_module
+  use iso_c_binding
+  
+  interface
+     
+     subroutine i2sort_c(A, iorder, isize) bind(C, name="qsort_int16_t")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer(c_int32_t)        :: iorder(isize)
+       integer   (c_int16_t)         :: A(isize)
+     end subroutine i2sort_c
+     
+     subroutine i2sort_noidx_c(A, isize) bind(C, name="qsort_int16_t_noidx")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer   (c_int16_t)         :: A(isize)
+     end subroutine i2sort_noidx_c
+     
+     
+     
+     subroutine i2sort_big_c(A, iorder, isize) bind(C, name="qsort_int16_t_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer(c_int64_t)        :: iorder(isize)
+       integer   (c_int16_t)         :: A(isize)
+     end subroutine i2sort_big_c
+     
+     subroutine i2sort_noidx_big_c(A, isize) bind(C, name="qsort_int16_t_noidx_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer   (c_int16_t)         :: A(isize)
+     end subroutine i2sort_noidx_big_c
+     
+     
+     
+     subroutine isort_c(A, iorder, isize) bind(C, name="qsort_int32_t")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer(c_int32_t)        :: iorder(isize)
+       integer   (c_int32_t)         :: A(isize)
+     end subroutine isort_c
+     
+     subroutine isort_noidx_c(A, isize) bind(C, name="qsort_int32_t_noidx")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer   (c_int32_t)         :: A(isize)
+     end subroutine isort_noidx_c
+     
+     
+     
+     subroutine isort_big_c(A, iorder, isize) bind(C, name="qsort_int32_t_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer(c_int64_t)        :: iorder(isize)
+       integer   (c_int32_t)         :: A(isize)
+     end subroutine isort_big_c
+     
+     subroutine isort_noidx_big_c(A, isize) bind(C, name="qsort_int32_t_noidx_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer   (c_int32_t)         :: A(isize)
+     end subroutine isort_noidx_big_c
+     
+     
+     
+     subroutine i8sort_c(A, iorder, isize) bind(C, name="qsort_int64_t")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer(c_int32_t)        :: iorder(isize)
+       integer   (c_int64_t)         :: A(isize)
+     end subroutine i8sort_c
+     
+     subroutine i8sort_noidx_c(A, isize) bind(C, name="qsort_int64_t_noidx")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer   (c_int64_t)         :: A(isize)
+     end subroutine i8sort_noidx_c
+     
+     
+     
+     subroutine i8sort_big_c(A, iorder, isize) bind(C, name="qsort_int64_t_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer(c_int64_t)        :: iorder(isize)
+       integer   (c_int64_t)         :: A(isize)
+     end subroutine i8sort_big_c
+     
+     subroutine i8sort_noidx_big_c(A, isize) bind(C, name="qsort_int64_t_noidx_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer   (c_int64_t)         :: A(isize)
+     end subroutine i8sort_noidx_big_c
+     
+     
+     
+     subroutine dsort_c(A, iorder, isize) bind(C, name="qsort_double")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer(c_int32_t)        :: iorder(isize)
+       real   (c_double)         :: A(isize)
+     end subroutine dsort_c
+     
+     subroutine dsort_noidx_c(A, isize) bind(C, name="qsort_double_noidx")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       real   (c_double)         :: A(isize)
+     end subroutine dsort_noidx_c
+     
+     
+     
+     subroutine dsort_big_c(A, iorder, isize) bind(C, name="qsort_double_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer(c_int64_t)        :: iorder(isize)
+       real   (c_double)         :: A(isize)
+     end subroutine dsort_big_c
+     
+     subroutine dsort_noidx_big_c(A, isize) bind(C, name="qsort_double_noidx_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       real   (c_double)         :: A(isize)
+     end subroutine dsort_noidx_big_c
+     
+     
+     
+     subroutine sort_c(A, iorder, isize) bind(C, name="qsort_float")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       integer(c_int32_t)        :: iorder(isize)
+       real   (c_float)         :: A(isize)
+     end subroutine sort_c
+     
+     subroutine sort_noidx_c(A, isize) bind(C, name="qsort_float_noidx")
+       use iso_c_binding
+       integer(c_int32_t), value :: isize
+       real   (c_float)         :: A(isize)
+     end subroutine sort_noidx_c
+     
+     
+     
+     subroutine sort_big_c(A, iorder, isize) bind(C, name="qsort_float_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       integer(c_int64_t)        :: iorder(isize)
+       real   (c_float)         :: A(isize)
+     end subroutine sort_big_c
+     
+     subroutine sort_noidx_big_c(A, isize) bind(C, name="qsort_float_noidx_big")
+       use iso_c_binding
+       integer(c_int64_t), value :: isize
+       real   (c_float)         :: A(isize)
+     end subroutine sort_noidx_big_c
+     
+     
+     
+  end interface
+
+end module qsort_module
+
+
+subroutine i2sort(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int32_t)        :: isize
+  integer(c_int32_t)        :: iorder(isize)
+  integer   (c_int16_t)         :: A(isize)
+  call i2sort_c(A, iorder, isize)
+end subroutine i2sort
+
+subroutine i2sort_noidx(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int32_t) :: isize
+  integer   (c_int16_t)    :: A(isize)
+  call i2sort_noidx_c(A, isize)
+end subroutine i2sort_noidx
+
+
+
+subroutine i2sort_big(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int64_t)        :: isize
+  integer(c_int64_t)        :: iorder(isize)
+  integer   (c_int16_t)         :: A(isize)
+  call i2sort_big_c(A, iorder, isize)
+end subroutine i2sort_big
+
+subroutine i2sort_noidx_big(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int64_t) :: isize
+  integer   (c_int16_t)    :: A(isize)
+  call i2sort_noidx_big_c(A, isize)
+end subroutine i2sort_noidx_big
+
+
+
+subroutine isort(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int32_t)        :: isize
+  integer(c_int32_t)        :: iorder(isize)
+  integer   (c_int32_t)         :: A(isize)
+  call isort_c(A, iorder, isize)
+end subroutine isort
+
+subroutine isort_noidx(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int32_t) :: isize
+  integer   (c_int32_t)    :: A(isize)
+  call isort_noidx_c(A, isize)
+end subroutine isort_noidx
+
+
+
+subroutine isort_big(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int64_t)        :: isize
+  integer(c_int64_t)        :: iorder(isize)
+  integer   (c_int32_t)         :: A(isize)
+  call isort_big_c(A, iorder, isize)
+end subroutine isort_big
+
+subroutine isort_noidx_big(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int64_t) :: isize
+  integer   (c_int32_t)    :: A(isize)
+  call isort_noidx_big_c(A, isize)
+end subroutine isort_noidx_big
+
+
+
+subroutine i8sort(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int32_t)        :: isize
+  integer(c_int32_t)        :: iorder(isize)
+  integer   (c_int64_t)         :: A(isize)
+  call i8sort_c(A, iorder, isize)
+end subroutine i8sort
+
+subroutine i8sort_noidx(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int32_t) :: isize
+  integer   (c_int64_t)    :: A(isize)
+  call i8sort_noidx_c(A, isize)
+end subroutine i8sort_noidx
+
+
+
+subroutine i8sort_big(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int64_t)        :: isize
+  integer(c_int64_t)        :: iorder(isize)
+  integer   (c_int64_t)         :: A(isize)
+  call i8sort_big_c(A, iorder, isize)
+end subroutine i8sort_big
+
+subroutine i8sort_noidx_big(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int64_t) :: isize
+  integer   (c_int64_t)    :: A(isize)
+  call i8sort_noidx_big_c(A, isize)
+end subroutine i8sort_noidx_big
+
+
+
+subroutine dsort(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int32_t)        :: isize
+  integer(c_int32_t)        :: iorder(isize)
+  real   (c_double)         :: A(isize)
+  call dsort_c(A, iorder, isize)
+end subroutine dsort
+
+subroutine dsort_noidx(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int32_t) :: isize
+  real   (c_double)    :: A(isize)
+  call dsort_noidx_c(A, isize)
+end subroutine dsort_noidx
+
+
+
+subroutine dsort_big(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int64_t)        :: isize
+  integer(c_int64_t)        :: iorder(isize)
+  real   (c_double)         :: A(isize)
+  call dsort_big_c(A, iorder, isize)
+end subroutine dsort_big
+
+subroutine dsort_noidx_big(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int64_t) :: isize
+  real   (c_double)    :: A(isize)
+  call dsort_noidx_big_c(A, isize)
+end subroutine dsort_noidx_big
+
+
+
+subroutine sort(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int32_t)        :: isize
+  integer(c_int32_t)        :: iorder(isize)
+  real   (c_float)         :: A(isize)
+  call sort_c(A, iorder, isize)
+end subroutine sort
+
+subroutine sort_noidx(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int32_t) :: isize
+  real   (c_float)    :: A(isize)
+  call sort_noidx_c(A, isize)
+end subroutine sort_noidx
+
+
+
+subroutine sort_big(A, iorder, isize) 
+  use qsort_module
+  use iso_c_binding
+  integer(c_int64_t)        :: isize
+  integer(c_int64_t)        :: iorder(isize)
+  real   (c_float)         :: A(isize)
+  call sort_big_c(A, iorder, isize)
+end subroutine sort_big
+
+subroutine sort_noidx_big(A, isize)
+  use iso_c_binding
+  use qsort_module
+  integer(c_int64_t) :: isize
+  real   (c_float)    :: A(isize)
+  call sort_noidx_big_c(A, isize)
+end subroutine sort_noidx_big
--- a/src/utils/sort.irp.f
+++ b/src/utils/sort.irp.f
@ -1,222 +1,4 @@
 BEGIN_TEMPLATE
- subroutine insertion_$Xsort (x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize) using the insertion sort algorithm.
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-  $type                          :: xtmp
-  integer                        :: i, i0, j, jmax
-
-  do i=2,isize
-    xtmp = x(i)
-    i0 = iorder(i)
-    j=i-1
-    do while (j>0)
-      if ((x(j) <= xtmp)) exit
-      x(j+1) = x(j)
-      iorder(j+1) = iorder(j)
-      j=j-1
-    enddo
-    x(j+1) = xtmp
-    iorder(j+1) = i0
-  enddo
- end subroutine insertion_$Xsort
-
- subroutine quick_$Xsort(x, iorder, isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize) using the quicksort algorithm.
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-  integer, external              :: omp_get_num_threads
-  call rec_$X_quicksort(x,iorder,isize,1,isize,nproc)
- end
-
- recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level)
-  implicit none
-  integer, intent(in)            :: isize, first, last, level
-  integer,intent(inout)          :: iorder(isize)
-  $type, intent(inout)           :: x(isize)
-  $type                          :: c, tmp
-  integer                        :: itmp
-  integer                        :: i, j
-
-  if(isize<2)return
-
-  c = x( shiftr(first+last,1) )
-  i = first
-  j = last
-  do
-    do while (x(i) < c)
-      i=i+1
-    end do
-    do while (c < x(j))
-      j=j-1
-    end do
-    if (i >= j) exit
-    tmp  = x(i)
-    x(i) = x(j)
-    x(j) = tmp
-    itmp      = iorder(i)
-    iorder(i) = iorder(j)
-    iorder(j) = itmp
-    i=i+1
-    j=j-1
-  enddo
-  if ( ((i-first <= 10000).and.(last-j <= 10000)).or.(level<=0) ) then
-    if (first < i-1) then
-      call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2)
-    endif
-    if (j+1 < last) then
-      call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2)
-    endif
-  else
-    if (first < i-1) then
-      call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2)
-    endif
-    if (j+1 < last) then
-      call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2)
-    endif
-  endif
- end
-
- subroutine heap_$Xsort(x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize) using the heap sort algorithm.
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-
-  integer                        :: i, k, j, l, i0
-  $type                          :: xtemp
-
-  l = isize/2+1
-  k = isize
-  do while (.True.)
-    if (l>1) then
-      l=l-1
-      xtemp = x(l)
-      i0 = iorder(l)
-    else
-      xtemp = x(k)
-      i0 = iorder(k)
-      x(k) = x(1)
-      iorder(k) = iorder(1)
-      k = k-1
-      if (k == 1) then
-        x(1) = xtemp
-        iorder(1) = i0
-        exit
-      endif
-    endif
-    i=l
-    j = shiftl(l,1)
-    do while (j<k)
-      if ( x(j) < x(j+1) ) then
-        j=j+1
-      endif
-      if (xtemp < x(j)) then
-        x(i) = x(j)
-        iorder(i) = iorder(j)
-        i = j
-        j = shiftl(j,1)
-      else
-        j = k+1
-      endif
-    enddo
-    if (j==k) then
-      if (xtemp < x(j)) then
-        x(i) = x(j)
-        iorder(i) = iorder(j)
-        i = j
-        j = shiftl(j,1)
-      else
-        j = k+1
-      endif
-    endif
-    x(i) = xtemp
-    iorder(i) = i0
-  enddo
- end subroutine heap_$Xsort
-
- subroutine heap_$Xsort_big(x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize) using the heap sort algorithm.
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  ! This is a version for very large arrays where the indices need
-  ! to be in integer*8 format
-  END_DOC
-  integer*8,intent(in)           :: isize
-  $type,intent(inout)            :: x(isize)
-  integer*8,intent(inout)        :: iorder(isize)
-
-  integer*8                      :: i, k, j, l, i0
-  $type                          :: xtemp
-
-  l = isize/2+1
-  k = isize
-  do while (.True.)
-    if (l>1) then
-      l=l-1
-      xtemp = x(l)
-      i0 = iorder(l)
-    else
-      xtemp = x(k)
-      i0 = iorder(k)
-      x(k) = x(1)
-      iorder(k) = iorder(1)
-      k = k-1
-      if (k == 1) then
-        x(1) = xtemp
-        iorder(1) = i0
-        exit
-      endif
-    endif
-    i=l
-    j = shiftl(l,1)
-    do while (j<k)
-      if ( x(j) < x(j+1) ) then
-        j=j+1
-      endif
-      if (xtemp < x(j)) then
-        x(i) = x(j)
-        iorder(i) = iorder(j)
-        i = j
-        j = shiftl(j,1)
-      else
-        j = k+1
-      endif
-    enddo
-    if (j==k) then
-      if (xtemp < x(j)) then
-        x(i) = x(j)
-        iorder(i) = iorder(j)
-        i = j
-        j = shiftl(j,1)
-      else
-        j = k+1
-      endif
-    endif
-    x(i) = xtemp
-    iorder(i) = i0
-  enddo
-
- end subroutine heap_$Xsort_big

 subroutine sorted_$Xnumber(x,isize,n)
  implicit none
@ -250,222 +32,6 @@ SUBST [ X, type ]
 END_TEMPLATE


-!---------------------- INTEL
-IRP_IF INTEL
-
-BEGIN_TEMPLATE
- subroutine $Xsort(x,iorder,isize)
-  use intel
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize).
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-  integer                        :: n
-  character, allocatable         :: tmp(:)
-  if (isize < 2) return
-  call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
-  allocate(tmp(n))
-  call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp)
-  deallocate(tmp)
-  iorder(1:isize) = iorder(1:isize)+1
-  call $Xset_order(x,iorder,isize)
- end
-
- subroutine $Xsort_noidx(x,isize)
-  use intel
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize).
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer                        :: n
-  character, allocatable         :: tmp(:)
-  if (isize < 2) return
-  call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
-  allocate(tmp(n))
-  call ippsSortRadixAscend_$ityp_I(x, isize, tmp)
-  deallocate(tmp)
- end
-
-SUBST [ X, type, ityp, n, ippsz ]
-   ; real ; 32f ; 4 ; 13 ;;
- i ; integer ; 32s ; 4 ; 11 ;;
- i2 ; integer*2 ; 16s ; 2 ; 7 ;;
-END_TEMPLATE
-
-BEGIN_TEMPLATE
-
- subroutine $Xsort(x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize).
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-  integer                        :: n
-  if (isize < 2) then
-    return
-  endif
-!  call sorted_$Xnumber(x,isize,n)
-!  if (isize == n) then
-!    return
-!  endif
-  if ( isize < 32) then
-    call insertion_$Xsort(x,iorder,isize)
-  else
-!    call heap_$Xsort(x,iorder,isize)
-    call quick_$Xsort(x,iorder,isize)
-  endif
- end subroutine $Xsort
-
-SUBST [ X, type ]
- d ; double precision ;;
-END_TEMPLATE
-
-BEGIN_TEMPLATE
-
- subroutine $Xsort(x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize).
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-  integer                        :: n
-  if (isize < 2) then
-    return
-  endif
-  call sorted_$Xnumber(x,isize,n)
-  if (isize == n) then
-    return
-  endif
-  if ( isize < 32) then
-    call insertion_$Xsort(x,iorder,isize)
-  else
-!    call $Xradix_sort(x,iorder,isize,-1)
-    call quick_$Xsort(x,iorder,isize)
-  endif
- end subroutine $Xsort
-
-SUBST [ X, type ]
- i8 ; integer*8 ;;
-END_TEMPLATE
-
-!---------------------- END INTEL
-IRP_ELSE
-!---------------------- NON-INTEL
-BEGIN_TEMPLATE
-
- subroutine $Xsort_noidx(x,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize).
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer, allocatable           :: iorder(:)
-  integer                        :: i
-  allocate(iorder(isize))
-  do i=1,isize
-   iorder(i)=i
-  enddo
-  call $Xsort(x,iorder,isize)
-  deallocate(iorder)
- end subroutine $Xsort_noidx
-
-SUBST [ X, type ]
-   ; real ;;
- d ; double precision ;;
- i ; integer ;;
- i8 ; integer*8 ;;
- i2 ; integer*2 ;;
-END_TEMPLATE
-
-BEGIN_TEMPLATE
-
- subroutine $Xsort(x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize).
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-  integer                        :: n
-  if (isize < 2) then
-    return
-  endif
-!  call sorted_$Xnumber(x,isize,n)
-!  if (isize == n) then
-!    return
-!  endif
-  if ( isize < 32) then
-    call insertion_$Xsort(x,iorder,isize)
-  else
-!    call heap_$Xsort(x,iorder,isize)
-    call quick_$Xsort(x,iorder,isize)
-  endif
- end subroutine $Xsort
-
-SUBST [ X, type ]
-   ; real ;;
- d ; double precision ;;
-END_TEMPLATE
-
-BEGIN_TEMPLATE
-
- subroutine $Xsort(x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize).
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  END_DOC
-  integer,intent(in)             :: isize
-  $type,intent(inout)            :: x(isize)
-  integer,intent(inout)          :: iorder(isize)
-  integer                        :: n
-  if (isize < 2) then
-    return
-  endif
-  call sorted_$Xnumber(x,isize,n)
-  if (isize == n) then
-    return
-  endif
-  if ( isize < 32) then
-    call insertion_$Xsort(x,iorder,isize)
-  else
-!    call $Xradix_sort(x,iorder,isize,-1)
-    call quick_$Xsort(x,iorder,isize)
-  endif
- end subroutine $Xsort
-
-SUBST [ X, type ]
- i ; integer ;;
- i8 ; integer*8 ;;
- i2 ; integer*2 ;;
-END_TEMPLATE
-
-IRP_ENDIF
-!---------------------- END NON-INTEL
-
-

 BEGIN_TEMPLATE
 subroutine $Xset_order(x,iorder,isize)
@ -491,47 +57,6 @@ BEGIN_TEMPLATE
  deallocate(xtmp)
 end

-SUBST [ X, type ]
-   ; real ;;
- d ; double precision ;;
- i ; integer ;;
- i8; integer*8 ;;
- i2; integer*2 ;;
-END_TEMPLATE
-
-
-BEGIN_TEMPLATE
- subroutine insertion_$Xsort_big (x,iorder,isize)
-  implicit none
-  BEGIN_DOC
-  ! Sort array x(isize) using the insertion sort algorithm.
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  ! This is a version for very large arrays where the indices need
-  ! to be in integer*8 format
-  END_DOC
-  integer*8,intent(in)           :: isize
-  $type,intent(inout)            :: x(isize)
-  integer*8,intent(inout)        :: iorder(isize)
-  $type                          :: xtmp
-  integer*8                      :: i, i0, j, jmax
-
-  do i=2_8,isize
-    xtmp = x(i)
-    i0 = iorder(i)
-    j = i-1_8
-    do while (j>0_8)
-      if (x(j)<=xtmp) exit
-      x(j+1_8) = x(j)
-      iorder(j+1_8) = iorder(j)
-      j = j-1_8
-    enddo
-    x(j+1_8) = xtmp
-    iorder(j+1_8) = i0
-  enddo
-
- end subroutine insertion_$Xsort_big
-
 subroutine $Xset_order_big(x,iorder,isize)
  implicit none
  BEGIN_DOC
@ -565,223 +90,3 @@ SUBST [ X, type ]
 END_TEMPLATE


-BEGIN_TEMPLATE
-
-recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix)
-  implicit none
-
-  BEGIN_DOC
-  ! Sort integer array x(isize) using the radix sort algorithm.
-  ! iorder in input should be (1,2,3,...,isize), and in output
-  ! contains the new order of the elements.
-  ! iradix should be -1 in input.
-  END_DOC
-  integer*$int_type, intent(in)  :: isize
-  integer*$int_type, intent(inout) :: iorder(isize)
-  integer*$type, intent(inout)   :: x(isize)
-  integer, intent(in)            :: iradix
-  integer                        :: iradix_new
-  integer*$type, allocatable     :: x2(:), x1(:)
-  integer*$type                  :: i4               ! data type
-  integer*$int_type, allocatable :: iorder1(:),iorder2(:)
-  integer*$int_type              :: i0, i1, i2, i3, i ! index type
-  integer*$type                  :: mask
-  integer                        :: err
-  !DIR$ ATTRIBUTES ALIGN : 128   :: iorder1,iorder2, x2, x1
-
-  if (isize < 2) then
-    return
-  endif
-
-  if (iradix == -1) then ! Sort Positive and negative
-
-    allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err)
-    if (err /= 0) then
-      print *,  irp_here, ': Unable to allocate arrays'
-      stop
-    endif
-
-    i1=1_$int_type
-    i2=1_$int_type
-    do i=1_$int_type,isize
-      if (x(i) < 0_$type) then
-        iorder1(i1) = iorder(i)
-        x1(i1) = -x(i)
-        i1 = i1+1_$int_type
-      else
-        iorder2(i2) = iorder(i)
-        x2(i2) = x(i)
-        i2 = i2+1_$int_type
-      endif
-    enddo
-    i1=i1-1_$int_type
-    i2=i2-1_$int_type
-
-    do i=1_$int_type,i2
-      iorder(i1+i) = iorder2(i)
-      x(i1+i) = x2(i)
-    enddo
-    deallocate(x2,iorder2,stat=err)
-    if (err /= 0) then
-      print *,  irp_here, ': Unable to deallocate arrays x2, iorder2'
-      stop
-    endif
-
-
-    if (i1 > 1_$int_type) then
-      call $Xradix_sort$big(x1,iorder1,i1,-2)
-      do i=1_$int_type,i1
-        x(i) = -x1(1_$int_type+i1-i)
-        iorder(i) = iorder1(1_$int_type+i1-i)
-      enddo
-    endif
-
-    if (i2>1_$int_type) then
-      call $Xradix_sort$big(x(i1+1_$int_type),iorder(i1+1_$int_type),i2,-2)
-    endif
-
-    deallocate(x1,iorder1,stat=err)
-    if (err /= 0) then
-      print *,  irp_here, ': Unable to deallocate arrays x1, iorder1'
-      stop
-    endif
-    return
-
-  else if (iradix == -2) then ! Positive
-
-    ! Find most significant bit
-
-    i0 = 0_$int_type
-    i4 = maxval(x)
-
-    iradix_new = max($integer_size-1-leadz(i4),1)
-    mask = ibset(0_$type,iradix_new)
-
-    allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err)
-    if (err /= 0) then
-      print *,  irp_here, ': Unable to allocate arrays'
-      stop
-    endif
-
-    i1=1_$int_type
-    i2=1_$int_type
-
-    do i=1_$int_type,isize
-      if (iand(mask,x(i)) == 0_$type) then
-        iorder1(i1) = iorder(i)
-        x1(i1) = x(i)
-        i1 = i1+1_$int_type
-      else
-        iorder2(i2) = iorder(i)
-        x2(i2) = x(i)
-        i2 = i2+1_$int_type
-      endif
-    enddo
-    i1=i1-1_$int_type
-    i2=i2-1_$int_type
-
-    do i=1_$int_type,i1
-      iorder(i0+i) = iorder1(i)
-      x(i0+i) = x1(i)
-    enddo
-    i0 = i0+i1
-    i3 = i0
-    deallocate(x1,iorder1,stat=err)
-    if (err /= 0) then
-      print *,  irp_here, ': Unable to deallocate arrays x1, iorder1'
-      stop
-    endif
-
-
-    do i=1_$int_type,i2
-      iorder(i0+i) = iorder2(i)
-      x(i0+i) = x2(i)
-    enddo
-    i0 = i0+i2
-    deallocate(x2,iorder2,stat=err)
-    if (err /= 0) then
-      print *,  irp_here, ': Unable to deallocate arrays x2, iorder2'
-      stop
-    endif
-
-
-    if (i3>1_$int_type) then
-      call $Xradix_sort$big(x,iorder,i3,iradix_new-1)
-    endif
-
-    if (isize-i3>1_$int_type) then
-      call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1)
-    endif
-
-    return
-  endif
-
-  ASSERT (iradix >= 0)
-
-  if (isize < 48) then
-    call insertion_$Xsort$big(x,iorder,isize)
-    return
-  endif
-
-
-  allocate(x2(isize),iorder2(isize),stat=err)
-  if (err /= 0) then
-    print *,  irp_here, ': Unable to allocate arrays x1, iorder1'
-    stop
-  endif
-
-
-  mask = ibset(0_$type,iradix)
-  i0=1_$int_type
-  i1=1_$int_type
-
-  do i=1_$int_type,isize
-    if (iand(mask,x(i)) == 0_$type) then
-      iorder(i0) = iorder(i)
-      x(i0) = x(i)
-      i0 = i0+1_$int_type
-    else
-      iorder2(i1) = iorder(i)
-      x2(i1) = x(i)
-      i1 = i1+1_$int_type
-    endif
-  enddo
-  i0=i0-1_$int_type
-  i1=i1-1_$int_type
-
-  do i=1_$int_type,i1
-    iorder(i0+i) = iorder2(i)
-    x(i0+i) = x2(i)
-  enddo
-
-  deallocate(x2,iorder2,stat=err)
-  if (err /= 0) then
-    print *,  irp_here, ': Unable to allocate arrays x2, iorder2'
-    stop
-  endif
-
-
-  if (iradix == 0) then
-    return
-  endif
-
-
-  if (i1>1_$int_type) then
-    call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1)
-  endif
-  if (i0>1) then
-    call $Xradix_sort$big(x,iorder,i0,iradix-1)
-  endif
-
- end
-
-SUBST [ X, type, integer_size, is_big, big, int_type ]
- i  ; 4 ; 32 ; .False. ;      ; 4 ;;
- i8 ; 8 ; 64 ; .False. ;      ; 4 ;;
- i2 ; 2 ; 16 ; .False. ;      ; 4 ;;
- i  ; 4 ; 32 ; .True.  ; _big ; 8 ;;
- i8 ; 8 ; 64 ; .True.  ; _big ; 8 ;;
-END_TEMPLATE
-
-
-
--- a/src/utils/units.irp.f
+++ b/src/utils/units.irp.f
@ -0,0 +1,22 @@
+BEGIN_PROVIDER [double precision, ha_to_ev]
+
+  implicit none
+  BEGIN_DOC
+  ! Converstion from Hartree to eV
+  END_DOC
+
+  ha_to_ev = 27.211396641308d0
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, au_to_D]
+
+  implicit none
+  BEGIN_DOC
+  ! Converstion from au to Debye
+  END_DOC
+
+  au_to_D = 2.5415802529d0
+
+END_PROVIDER
+
--- a/src/utils/util.irp.f
+++ b/src/utils/util.irp.f
@ -37,6 +37,10 @@ double precision function binom_func(i,j)
  else
    binom_func = dexp( logfact(i)-logfact(j)-logfact(i-j) )
  endif
+
+  ! To avoid .999999 numbers
+  binom_func = floor(binom_func + 0.5d0)
+
 end


--- a/tests/bats/common.bats.sh
+++ b/tests/bats/common.bats.sh
@ -46,7 +46,7 @@ function test_exe() {

 run_only_test() {
  if [[ "$BATS_TEST_DESCRIPTION" != "$1" ]] && [[ "$BATS_TEST_NUMBER" != "$1" ]]; then
-    if [[ -z $BATS_TEST_FILENAME ]] ; then
+    if [[ -z "$BATS_TEST_FILENAME" ]] ; then
      exit 0
    else
      skip
--- a/travis/compilation.sh
+++ b/travis/compilation.sh
@ -1,16 +0,0 @@
-#!/bin/bash
-# Stage 2
-
-# Extract cache from config stage
-cd ../
-tar -zxf $HOME/cache/config.tgz
-
-# Configure QP2
-cd qp2
-source ./quantum_package.rc
-ninja -j 1 -v || exit -1
-
-# Create cache
-cd ..
-tar -zcf $HOME/cache/compil.tgz qp2 && rm $HOME/cache/config.tgz
-
--- a/travis/configuration.sh
+++ b/travis/configuration.sh
@ -1,10 +0,0 @@
-#!/bin/bash
-# Stage 1
-
-# Configure QP2
-./configure --download all --install all --config ./config/travis.cfg || exit -1
-
-# Create cache
-cd ../
-tar -zcf $HOME/cache/config.tgz qp2
-
--- a/travis/testing.sh
+++ b/travis/testing.sh
@ -1,16 +0,0 @@
-#!/bin/bash
-# Stage 3
-
-# Extract cache from compile stage
-cd ../
-tar -zxf $HOME/cache/compil.tgz
-
-# Configure QP2
-cd qp2
-source ./quantum_package.rc
-exec qp_test -a && rm $HOME/cache/compil.tgz
-
-
-
-
-