Merge branch 'master' into feature/4idx

Conflicts:
	.travis.yml
	src/AO_Basis/aos.irp.f

.travis.yml

@@ -4,17 +4,21 @@
 #  - sudo apt-get install gfortran liblapack-dev gcc
 #  - sudo apt-get install graphviz
 
+os: linux
+
 dist: trusty
 
 sudo: false
 
+compiler: gfortran
+
 addons:
   apt:
     packages:
     - gfortran
     - gcc
-    - libblas-dev
+#    - liblapack-dev
-    - liblapack-dev
+#    - libblas-dev
     - graphviz
 
 cache:
@@ -28,6 +32,7 @@ python:
 script: 
   - ./configure  ./config/travis.cfg
   - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD_ZMQ mrcepa0 All_singles
+  - source ./quantum_package.rc ; cd ~ ; install_lapack.sh ; cd $QP_ROOT
   - source ./quantum_package.rc ; ninja
   - source ./quantum_package.rc ; cd ocaml ; make ; cd -
   - source ./quantum_package.rc ; cd tests ; ./run_tests.sh -v

README.md

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

config/gfortran.cfg

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

config/gfortran_debug.cfg

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

config/travis.cfg

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

install/scripts/install_lapack.sh

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

install/scripts/install_ocaml.sh

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

ocaml/Address.ml

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

ocaml/Atom.ml

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

ocaml/Basis.ml

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

ocaml/Bit.ml

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

ocaml/Bit.mli

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

ocaml/Bitlist.ml

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

ocaml/Charge.ml

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

ocaml/Charge.mli

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

ocaml/Determinant.ml

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

ocaml/Determinant.mli

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

ocaml/Element.ml

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

ocaml/Element.mli

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

ocaml/Excitation.ml

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

ocaml/Excitation.mli

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

ocaml/GaussianPrimitive.ml

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

ocaml/Generic_input_of_rst.ml

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

ocaml/Gto.ml

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

ocaml/Gto.mli

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

ocaml/Input.ml

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

ocaml/Input_ao_basis.ml

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

ocaml/Input_bi_integrals.ml

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

ocaml/Input_bitmasks.ml

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

ocaml/Input_determinants_by_hand.ml

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

ocaml/Input_electrons.ml

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

ocaml/Input_mo_basis.ml

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

ocaml/Input_nuclei_by_hand.ml

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

ocaml/Io_ext.ml

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

ocaml/Long_basis.ml

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

ocaml/Long_basis.mli

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

ocaml/MO_class.ml

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

ocaml/MO_class.mli

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

ocaml/MO_label.ml

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

ocaml/MO_label.mli

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

ocaml/Makefile

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

ocaml/Message.ml

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

ocaml/Molecule.ml

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

ocaml/Molecule.mli

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

ocaml/Multiplicity.ml

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

ocaml/Multiplicity.mli

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

ocaml/Point3d.ml

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

ocaml/Point3d.mli

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

ocaml/Primitive.mli

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

ocaml/Progress_bar.ml

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

ocaml/Pseudo.ml

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

ocaml/Qpackage.ml

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

ocaml/Qputils.ml

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

ocaml/README_qp_edit.rst

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

ocaml/Range.ml

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

ocaml/Range.mli

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

ocaml/String_ext.ml

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

ocaml/Symmetry.ml

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

ocaml/Symmetry.mli

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

ocaml/TaskServer.ml

@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 
@@ -63,7 +63,7 @@ let bind_socket ~socket_type ~socket ~port =
         ZMQ.Socket.bind socket @@ Printf.sprintf "tcp://*:%d" port;
         loop (-1)
       with
-      | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) )
+      | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_sec 1. ; loop (i-1) )
       | other_exception -> raise other_exception
   in loop 60
 

ocaml/To_md5.ml

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

ocaml/_tags

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

ocaml/create_git_sha1.sh

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

ocaml/qp_basis_clean.ml

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

ocaml/qp_create_ezfio_from_xyz.ml

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

ocaml/qp_create_guess.ml

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

ocaml/qp_find_pi_space.ml

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

ocaml/qp_overlap_of_wf.ml

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

ocaml/qp_print.ml

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

ocaml/qp_print_basis.ml

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

ocaml/qp_run.ml

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

ocaml/qp_set_mo_class.ml

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

ocaml/qptypes_generator.ml

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

ocaml/test_basis.ml

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

ocaml/test_gto.ml

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

ocaml/test_message.ml

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

ocaml/test_molecule.ml

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

ocaml/test_pseudo.ml

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

ocaml/test_queuing_system.ml

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

ocaml/test_symmetry.ml

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

ocaml/test_task_server.ml

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

plugins/CAS_SD_ZMQ/selection.irp.f

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

plugins/FOBOCI/diag_fock_inactiv_virt.irp.f

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

plugins/FOBOCI/new_approach.irp.f

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

plugins/FOBOCI/routines_foboci.irp.f

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

plugins/Full_CI_ZMQ/selection.irp.f

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

plugins/Hartree_Fock/DIIS.irp.f

@@ -94,7 +94,7 @@ END_PROVIDER
 
   do i=1,AO_num
     do j=1,AO_num
-      Xt(i,j) = X_Matrix_AO(j,i)
+      Xt(i,j) = S_half_inv(j,i)
     enddo
   enddo
 
@@ -103,7 +103,7 @@ END_PROVIDER
   call dgemm('N','N',AO_num,AO_num,AO_num,     &
        1.d0,                                   &
        Fock_matrix_AO,size(Fock_matrix_AO,1),  &
-       X_Matrix_AO,size(X_Matrix_AO,1),        &
+       S_half_inv,size(S_half_inv,1),        &
        0.d0,                                   &
        eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))       
 
@@ -130,67 +130,10 @@ END_PROVIDER
 
   call dgemm('N','N',AO_num,AO_num,AO_num,     &
        1.d0,                                   &
-       X_matrix_AO,size(X_matrix_AO,1),        &
+       S_half_inv,size(S_half_inv,1),        &
        scratch,size(scratch,1),                &
        0.d0,                                   &
        eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))       
    
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num,AO_num) ]
-
-  BEGIN_DOC
-!   Matrix X = S^{-1/2} obtained by SVD
-  END_DOC
-
-  implicit none
-  
-  integer                         :: num_linear_dependencies
-  integer                         :: LDA, LDC
-  double precision, allocatable   :: U(:,:),Vt(:,:), D(:)
-  integer                         :: info, i, j, k
- 
-  LDA = size(AO_overlap,1)
-  LDC = size(X_matrix_AO,1)
- 
-  allocate(         &
-    U(LDC,AO_num),  &
-    Vt(LDA,AO_num), &
-    D(AO_num))
-
-  call svd(            &
-       AO_overlap,LDA, &
-       U,LDC,          &
-       D,              &
-       Vt,LDA,         &
-       AO_num,AO_num)
-
-  num_linear_dependencies = 0
-  do i=1,AO_num
-    print*,D(i)
-    if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
-      D(i) = 0.d0
-      num_linear_dependencies += 1
-    else
-      ASSERT (D(i) > 0.d0)
-      D(i) = 1.d0/sqrt(D(i))
-    endif
-    do j=1,AO_num
-      X_matrix_AO(j,i) = 0.d0
-    enddo
-  enddo
-  write(*,*) 'linear dependencies',num_linear_dependencies
-!  stop
-
-  do k=1,AO_num
-    if(D(k) /= 0.d0) then
-      do j=1,AO_num
-        do i=1,AO_num
-          X_matrix_AO(i,j) = X_matrix_AO(i,j) + U(i,k)*D(k)*Vt(k,j)
-        enddo
-      enddo
-    endif
-  enddo
-  
-
-END_PROVIDER

plugins/Hartree_Fock/EZFIO.cfg

@@ -1,9 +1,3 @@
-[threshold_overlap_ao_eigenvalues]
-type: Threshold
-doc: Threshold on the magnitude of the smallest eigenvalues of the overlap matrix in the AO basis
-interface: ezfio,provider,ocaml
-default: 1.e-6
-
 [max_dim_diis]
 type: integer
 doc: Maximum size of the DIIS extrapolation procedure
@@ -32,7 +26,7 @@ default: 500
 type: Positive_float
 doc: Energy shift on the virtual MOs to improve SCF convergence
 interface: ezfio,provider,ocaml
-default: 0.0
+default: 0.1
 
 [scf_algorithm]
 type: character*(32)

plugins/Hartree_Fock/Fock_matrix.irp.f

@@ -263,17 +263,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num,mo_tot_num)
    BEGIN_DOC
    ! Fock matrix on the MO basis
    END_DOC
-   double precision, allocatable  :: T(:,:)
+   call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), &
-   allocate ( T(ao_num,mo_tot_num) )
+                 Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1))
-   call dgemm('N','N', ao_num, mo_tot_num, ao_num,                   &
-       1.d0, Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1),      &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0, T, size(T,1))
-   call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num,               &
-       1.d0, mo_coef,size(mo_coef,1),                                &
-       T, size(T,1),                                                 &
-       0.d0, Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1))
-   deallocate(T)
 END_PROVIDER
  
  
@@ -282,17 +273,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num,mo_tot_num)
    BEGIN_DOC
    ! Fock matrix on the MO basis
    END_DOC
-   double precision, allocatable  :: T(:,:)
+   call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), &
-   allocate ( T(ao_num,mo_tot_num) )
+                 Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1))
-   call dgemm('N','N', ao_num, mo_tot_num, ao_num,                   &
-       1.d0, Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1),        &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0, T, size(T,1))
-   call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num,               &
-       1.d0, mo_coef,size(mo_coef,1),                                &
-       T, size(T,1),                                                 &
-       0.d0, Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1))
-   deallocate(T)
 END_PROVIDER
  
 BEGIN_PROVIDER [ double precision, HF_energy ]
@@ -330,97 +312,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ]
      enddo
    enddo
  else
-   double precision, allocatable  :: T(:,:), M(:,:)
+   call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), &
-   integer                        :: ierr
+      Fock_matrix_ao,size(Fock_matrix_ao,1)) 
-   ! F_ao = S C F_mo C^t S
-   allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr)
-   if (ierr /=0 ) then
-      print *,  irp_here, ' : allocation failed'
-   endif
-
-!  ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
-!  -> M(ao_num,mo_tot_num)
-   call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0,               &
-       ao_overlap, size(ao_overlap,1),                               &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0,                                                         &
-       M, size(M,1))
-
-!  M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num)
-!  -> T(ao_num,mo_tot_num)
-   call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,           &
-       M, size(M,1),                                                 &
-       Fock_matrix_mo, size(Fock_matrix_mo,1),                       &
-       0.d0,                                                         &
-       T, size(T,1))
-
-!  T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
-!  -> M(ao_num,ao_num)
-   call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0,           &
-       T, size(T,1),                                                 &
-       mo_coef, size(mo_coef,1),                                     &
-       0.d0,                                                         &
-       M, size(M,1))
-
-!  M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
-!  -> Fock_matrix_ao(ao_num,ao_num)
-   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
-       M, size(M,1),                                                 &
-       ao_overlap, size(ao_overlap,1),                               &
-       0.d0,                                                         &
-       Fock_matrix_ao, size(Fock_matrix_ao,1))
-
-
-   deallocate(T)
  endif
 END_PROVIDER
 
-subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
-  implicit none
-  integer, intent(in)            :: LDFMO ! size(FMO,1)
-  integer, intent(in)            :: LDFAO ! size(FAO,1)
-  double precision, intent(in)   :: FMO(LDFMO,*)
-  double precision, intent(out)  :: FAO(LDFAO,*)
-  
-  double precision, allocatable  :: T(:,:), M(:,:)
-  integer                        :: ierr
-  ! F_ao = S C F_mo C^t S
-  allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr)
-  if (ierr /=0 ) then
-     print *,  irp_here, ' : allocation failed'
-  endif
-
-!  ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
-!  -> M(ao_num,mo_tot_num)
-  call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0,                    &
-      ao_overlap, size(ao_overlap,1),                                &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0,                                                          &
-      M, size(M,1))
-
-!  M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num)
-!  -> T(ao_num,mo_tot_num)
-  call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
-      M, size(M,1),                                                  &
-      FMO, size(FMO,1),                                              &
-      0.d0,                                                          &
-      T, size(T,1))
-
-!  T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
-!  -> M(ao_num,ao_num)
-  call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0,            &
-      T, size(T,1),                                                  &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0,                                                          &
-      M, size(M,1))
-
-!  M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
-!  -> Fock_matrix_ao(ao_num,ao_num)
-  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
-      M, size(M,1),                                                  &
-      ao_overlap, size(ao_overlap,1),                                &
-      0.d0,                                                          &
-      FAO, size(FAO,1))
-  deallocate(T,M)
-end
 

plugins/Hartree_Fock/HF_density_matrix_ao.irp.f

@@ -1,7 +1,7 @@
 BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! S^-1 x Alpha density matrix in the AO basis x S^-1
+   ! S^{-1}.P_alpha.S^{-1}
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
@@ -14,7 +14,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! S^-1 Beta density matrix in the AO basis x S^-1
+   ! S^{-1}.P_beta.S^{-1}
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
@@ -27,7 +27,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! S^-1 Density matrix in the AO basis S^-1
+   ! S^{-1}.P.S^{-1}  where P = C.C^t
    END_DOC
    ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
    if (elec_alpha_num== elec_beta_num) then

plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f

@@ -133,7 +133,7 @@ END_DOC
   write(output_hartree_fock,*)
   
   if(.not.no_oa_or_av_opt)then
-   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
+   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
   endif
 
   call write_double(output_hartree_fock, Energy_SCF, 'Hartree-Fock energy')

plugins/Hartree_Fock/SCF.irp.f

@@ -23,7 +23,7 @@ subroutine create_guess
       mo_coef = ao_ortho_lowdin_coef
       TOUCH mo_coef
       mo_label = 'Guess'
-      call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label)
+      call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.)
       SOFT_TOUCH mo_coef mo_label
     else if (mo_guess_type == "Huckel") then
       call huckel_guess

plugins/Hartree_Fock/damping_SCF.irp.f

@@ -119,7 +119,7 @@ subroutine damping_SCF
   write(output_hartree_fock,*)
   
   if(.not.no_oa_or_av_opt)then
-   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
+   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
   endif
 
   call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')

plugins/Hartree_Fock/huckel.irp.f

@@ -7,25 +7,28 @@ subroutine huckel_guess
   double precision               :: accu
   double precision               :: c
   character*(64)                 :: label
-
+  double precision, allocatable  :: A(:,:)
   label = "Guess"
-  call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
-                                     size(mo_mono_elec_integral,1),  &
-                                     size(mo_mono_elec_integral,2),label,1)
-  TOUCH mo_coef
-
   c = 0.5d0 * 1.75d0
 
+  allocate (A(ao_num, ao_num))
+  A = 0.d0
   do j=1,ao_num
     do i=1,ao_num
-      Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + &
+      A(i,j) = c * ao_overlap(i,j) * (ao_mono_elec_integral_diag(i) + ao_mono_elec_integral_diag(j))
-                                                 ao_mono_elec_integral_diag(j))
     enddo
-    Fock_matrix_ao(j,j) = Fock_matrix_ao_alpha(j,j)
+    A(j,j) = ao_mono_elec_integral_diag(j) + ao_bi_elec_integral_alpha(j,j)
   enddo
-  TOUCH Fock_matrix_ao
+
+  Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
+  Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
+
+!  TOUCH mo_coef
+
+  TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta
   mo_coef = eigenvectors_fock_matrix_mo
   SOFT_TOUCH mo_coef
   call save_mos
+  deallocate(A)
 
 end

plugins/MRPT/jm_mrpt2.irp.f

@@ -11,7 +11,7 @@ end
 subroutine routine_3
  implicit none
 !provide fock_virt_total_spin_trace
- provide delta_ij 
+ provide delta_ij_mrpt 
  
  print *,  'N_det    = ', N_det
  print *,  'N_states = ', N_states

plugins/MRPT_Utils/mrpt_utils.irp.f

@@ -1,5 +1,5 @@
 
- BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
+ BEGIN_PROVIDER [ double precision, delta_ij_mrpt, (N_det,N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, second_order_pt_new, (N_states) ]
 &BEGIN_PROVIDER [ double precision, second_order_pt_new_1h, (N_states) ]
 &BEGIN_PROVIDER [ double precision, second_order_pt_new_1p, (N_states) ]
@@ -19,7 +19,7 @@
  double precision, allocatable :: delta_ij_tmp(:,:,:)
 
 
- delta_ij = 0.d0
+ delta_ij_mrpt = 0.d0
 
  allocate (delta_ij_tmp(N_det,N_det,N_states))
 
@@ -32,7 +32,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
    second_order_pt_new_1h(i_state) = accu(i_state) 
@@ -47,7 +47,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_1p(i_state) = accu(i_state) 
@@ -63,7 +63,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_1h1p(i_state) = accu(i_state) 
@@ -79,7 +79,7 @@
   do i = 1, N_det
    do j = 1, N_det
     accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-    delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+    delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
    enddo
   enddo
   second_order_pt_new_1h1p(i_state) = accu(i_state) 
@@ -95,7 +95,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_2h(i_state) = accu(i_state) 
@@ -110,7 +110,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_2p(i_state) = accu(i_state) 
@@ -126,7 +126,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_1h2p(i_state) = accu(i_state) 
@@ -142,7 +142,7 @@
  do i = 1, N_det
   do j = 1, N_det
    accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
   enddo
  enddo
  second_order_pt_new_2h1p(i_state) = accu(i_state) 
@@ -157,7 +157,7 @@
 !do i = 1, N_det
 ! do j = 1, N_det
 !  accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-!  delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+!  delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
 ! enddo
 !enddo
 !second_order_pt_new_2h2p(i_state) = accu(i_state) 
@@ -168,7 +168,7 @@
  call give_2h2p(contrib_2h2p)
  do i_state = 1, N_states
  do i = 1, N_det
-   delta_ij(i,i,i_state) += contrib_2h2p(i_state)
+   delta_ij_mrpt(i,i,i_state) += contrib_2h2p(i_state)
  enddo
  second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) 
  enddo
@@ -179,9 +179,9 @@
  accu = 0.d0
  do i_state = 1, N_states
  do i = 1, N_det
-! write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
+! write(*,'(1000(F16.10,x))')delta_ij_mrpt(i,:,:)
   do j = i_state, N_det
-   accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
+   accu(i_state) += delta_ij_mrpt(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
   enddo
  enddo
  second_order_pt_new(i_state) = accu(i_state) 
@@ -199,7 +199,7 @@ END_PROVIDER
  do i_state = 1, N_states
   do i = 1,N_det
    do j = 1,N_det
-    Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij(j,i,i_state)
+    Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij_mrpt(j,i,i_state)
    enddo
   enddo
  enddo
@@ -214,7 +214,7 @@ END_PROVIDER
   do i = 1,N_det
    do j = i,N_det
     Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) =  H_matrix_all_dets(j,i) & 
-                                            + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) )
+                                            + 0.5d0 * ( delta_ij_mrpt(j,i,i_state) + delta_ij_mrpt(i,j,i_state) )
     Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) =  Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) 
    enddo
   enddo

plugins/analyze_wf/attachment.irp.f

@@ -0,0 +1,64 @@
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_detachment, (mo_tot_num,mo_tot_num,2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_attachment, (mo_tot_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Detachment and attachment density matrices in MO basis
+ END_DOC
+ integer :: i,j, k, istate
+ double precision :: km(mo_tot_num), kp(mo_tot_num)
+
+ one_body_dm_mo_detachment = 0.d0
+ one_body_dm_mo_attachment = 0.d0
+
+ do istate=2,N_states
+
+   km(:) = 0.d0
+   kp(:) = 0.d0
+   do i=1,mo_tot_num
+      if (one_body_dm_mo_diff_eigvalues(i,istate) < 0) then
+        km(i) = -one_body_dm_mo_diff_eigvalues(i,istate)
+      else
+        kp(i) =  one_body_dm_mo_diff_eigvalues(i,istate)
+      endif
+   enddo
+
+   ! Attachment
+  do k=1,mo_tot_num
+     do j=1,mo_tot_num
+       do i=1,mo_tot_num
+         one_body_dm_mo_detachment(i,j,istate) = one_body_dm_mo_detachment(i,j,istate) + &
+           one_body_dm_mo_diff_eigvectors(i,k,istate)*km(k)* &
+           one_body_dm_mo_diff_eigvectors(j,k,istate)
+         one_body_dm_mo_attachment(i,j,istate) = one_body_dm_mo_attachment(i,j,istate) + &
+           one_body_dm_mo_diff_eigvectors(i,k,istate)*kp(k)* &
+           one_body_dm_mo_diff_eigvectors(j,k,istate)
+       enddo
+    enddo
+  enddo
+
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, one_body_dm_ao_detachment, (ao_num,ao_num,2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_ao_attachment, (ao_num,ao_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Detachment and attachment density matrices in AO basis
+ END_DOC
+ integer :: istate
+ do istate=2,N_states
+   call mo_to_ao_no_overlap(                                         &
+       one_body_dm_mo_attachment(1,1,istate),                        &
+       size(one_body_dm_mo_attachment,1),                            &
+       one_body_dm_ao_attachment(1,1,istate),                        &
+       size(one_body_dm_ao_attachment,1) )
+   call mo_to_ao_no_overlap(                                         &
+       one_body_dm_mo_detachment(1,1,istate),                        &
+       size(one_body_dm_mo_detachment,1),                            &
+       one_body_dm_ao_detachment(1,1,istate),                        &
+       size(one_body_dm_ao_detachment,1) )
+ enddo
+
+END_PROVIDER
+

plugins/analyze_wf/dump_nto.irp.f

@@ -0,0 +1,19 @@
+program dump_nto
+  implicit none
+  integer :: i,j, istate
+
+  print *,  'Phi_S'
+  do i=2,N_states
+    print *,  i, Phi_S(i)
+  enddo
+
+  do istate=2,N_states
+    do j=1,mo_tot_num
+      print *,  'MO: ', j, 'State:', istate, 'Eig:', one_body_dm_mo_diff_eigvalues(j,istate)
+      do i=1,ao_num
+        print *,  i, transition_natorb(i,j,istate)
+      enddo
+    enddo
+  enddo
+
+end

plugins/analyze_wf/dump_one_body_mos.irp.f

@@ -0,0 +1,38 @@
+program dump_one_body_mos
+  implicit none
+  BEGIN_DOC
+! Output density matrices of all the states
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run()
+end
+
+subroutine run
+  implicit none
+  integer                        :: istate
+  integer, parameter             :: iunit = 66
+  character*(64)                 :: filename, fmt
+  integer                        :: i,j,k
+
+  write(fmt,'(''('',I4.4,''(X,E20.14))'')') mo_tot_num
+  do istate=1,N_states
+    write(filename,'(''state.'',I2.2)') istate
+    open(unit=iunit, form='formatted', file=filename)
+    write(iunit,*) mo_tot_num
+    do j=1,mo_tot_num
+      write(iunit,fmt) one_body_dm_mo_alpha(1:mo_tot_num,j,istate) + one_body_dm_mo_beta(1:mo_tot_num,j,istate)
+    enddo
+  enddo
+
+
+  call run2()
+end
+subroutine run2
+  integer :: i,j, istate
+  print *,  'Phi_S'
+  do i=2,N_states
+    print *,  i, Phi_S(i)
+  enddo
+end
+

plugins/analyze_wf/phi_s.irp.f

@@ -0,0 +1,128 @@
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvalues,  (mo_tot_num, 2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvectors, (mo_tot_num, mo_tot_num, 2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Eigenvalues and eigenvectors of one_body_dm_mo_diff
+ END_DOC
+ integer                        :: i,j,istate
+ integer                        :: liwork, lwork, n, info
+ integer, allocatable           :: iwork(:)
+ double precision, allocatable  :: work(:)
+ 
+ 
+ one_body_dm_mo_diff_eigvectors(1:mo_tot_num, 1:mo_tot_num, 2:N_states) =&
+     one_body_dm_mo_diff(1:mo_tot_num, 1:mo_tot_num, 2:N_states)
+ 
+ n = mo_tot_num
+ lwork = 1+6*n + 2*n*n
+ liwork = 3 + 5*n
+ 
+ allocate(work(lwork), iwork(liwork))
+ 
+ lwork=-1
+ liwork=-1
+ istate=2
+ 
+ call dsyevd( 'V', 'U', mo_tot_num,                                  &
+     one_body_dm_mo_diff_eigvectors(1,1,istate),                     &
+     size(one_body_dm_mo_diff_eigvectors,1),                         &
+     one_body_dm_mo_diff_eigvalues(1,istate),                        &
+     work, lwork, iwork, liwork, info)
+ 
+ 
+ if (info /= 0) then
+   print *,  irp_here//' DSYEVD failed : ', info
+   stop 1
+ endif
+ lwork = int(work(1))
+ liwork = iwork(1)
+ deallocate(iwork,work)
+ 
+ allocate(work(lwork), iwork(liwork))
+ 
+ do istate=2,N_states
+   call dsyevd( 'V', 'U', mo_tot_num,                                &
+       one_body_dm_mo_diff_eigvectors(1,1,istate),                   &
+       size(one_body_dm_mo_diff_eigvectors,1),                       &
+       one_body_dm_mo_diff_eigvalues(1,istate),                      &
+       work, lwork, iwork, liwork, info)
+
+   if (info /= 0) then
+       print *,  irp_here//' DSYEVD failed : ', info
+       stop 1
+   endif
+
+ enddo
+
+ deallocate(iwork,work)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, transition_natorb, (ao_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Natural transition molecular orbitals
+ END_DOC
+
+ integer :: istate
+
+ do istate=2,N_states
+   call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
+       mo_coef, size(mo_coef,1),                                     &
+       one_body_dm_mo_diff_eigvectors(1,1,istate),                   &
+       size(one_body_dm_mo_diff_eigvectors,1), 0.d0,                 &
+       transition_natorb(1,1,istate), size(transition_natorb,1))
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, phi_s, (2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! 
+ END_DOC
+
+ integer :: i,istate
+ double precision, allocatable :: T(:,:), A(:,:), D(:,:)
+ double precision :: trace, norm
+ allocate(T(ao_num,ao_num), A(ao_num,ao_num), D(ao_num,ao_num))
+
+  do istate=2,N_states
+
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        S_half, size(S_half,1),                                      &
+        one_body_dm_ao_attachment(1,1,istate), size(one_body_dm_ao_attachment,1), 0.d0,&
+        T, size(T,1))
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        T, size(T,1),                                                &
+        S_half, size(S_half,1), 0.d0,                                &
+        A, size(A,1))
+!
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        S_half, size(S_half,1),                                      &
+        one_body_dm_ao_detachment(1,1,istate), size(one_body_dm_ao_detachment,1), 0.d0,&
+        T, size(T,1))
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        T, size(T,1),                                                &
+        S_half, size(S_half,1), 0.d0,                                &
+        D, size(D,1))
+
+    trace = 0.d0
+    do i=1,ao_num
+      trace = trace + A(i,i)
+    enddo
+    norm = 0.d0
+    do i=1,ao_num
+      norm = norm + D(i,i)
+    enddo
+    norm = 0.5d0*(norm + trace)
+
+    trace = 0.d0
+    do i=1,mo_tot_num
+      trace = trace + dsqrt(A(i,i)*D(i,i))
+    enddo
+    phi_s(istate) = trace/norm 
+  enddo
+
+END_PROVIDER
+