Merge pull request #1 from LCPQ/master

merge

Makefile

@@ -32,6 +32,7 @@ bin/irpf90:
 	  (echo Unable to download IRPF90 : $(WWW_SERVER)/$(IRPF90_TGZ) ; exit 1)
 	tar -zxf $(IRPF90_TGZ) && rm $(IRPF90_TGZ)
 	$(MAKE) -C irpf90 | tee install_irpf90.log
+	rm -rf -- $$PWD/bin/irpf90 $$PWD/bin/irpman
 	ln -s $$PWD/irpf90/bin/irpf90 $$PWD/bin/irpf90
 	ln -s $$PWD/irpf90/bin/irpman $$PWD/bin/irpman
 
@@ -53,7 +54,9 @@ bin/m4:
 	QPACKAGE_ROOT=$$PWD ./scripts/install_m4.sh | tee install_m4.log
 
 
-ocaml: ocaml/Qptypes.ml curl m4
+ocaml: curl m4
+	- rm -f -- ocaml/Qptypes.ml
+	$(MAKE) ocaml/Qptypes.ml
 
 ocaml/Qptypes.ml: 
 	$(info $(BLUE)===== Installing ocaml =====$(BLACK))

README.md

@@ -0,0 +1,8 @@
+Quantum package
+===============
+
+[![Gitter](https://badges.gitter.im/Join Chat.svg)](https://gitter.im/LCPQ/quantum_package?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)
+
+Set of quantum chemistry programs and libraries.
+
+For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>)

README.rst

@@ -1,9 +0,0 @@
-===============
-Quantum package
-===============
-
-Set of quantum chemistry programs and libraries.
-
-For more information, you can visit the
-`wiki of the project <http://github.com/LCPQ/quantum_package/wiki>`_
-

ocaml/Atom.ml

@@ -6,7 +6,7 @@ type t =
 { element : Element.t ;
   charge  : Charge.t ;
   coord   : Point3d.t ;
-}
+} with sexp
 
 (** Read xyz coordinates of the atom with unit u *)
 let of_string u s =
@@ -23,7 +23,7 @@ let of_string u s =
   | [ name; x; y; z ] ->
     let e = Element.of_string name in
     { element = e ;
-      charge  = Charge.of_int (Element.to_charge e);
+      charge  = Element.to_charge e;
       coord   = Point3d.of_string u (String.concat [x; y; z] ~sep:" ")
     }
   | _ -> raise (AtomError s)
@@ -33,6 +33,6 @@ let to_string u a =
   [ Element.to_string a.element ;
     Charge.to_string  a.charge ;
     Point3d.to_string u a.coord ]
-  |> String.concat ?sep:(Some "   ")
+  |> String.concat ~sep:"   "
 ;;
 

ocaml/Atom.mli

@@ -0,0 +1,9 @@
+exception AtomError of string
+
+type t = { element : Element.t; charge : Charge.t; coord : Point3d.t; } 
+
+val t_of_sexp : Sexplib.Sexp.t -> t
+val sexp_of_t : t -> Sexplib.Sexp.t
+
+val of_string : Units.units -> string -> t
+val to_string : Units.units -> t -> string

ocaml/Basis.ml

@@ -1,7 +1,7 @@
 open Core.Std;;
 open Qptypes;;
 
-type t = (Gto.t * Nucl_number.t) list;;
+type t = (Gto.t * Nucl_number.t) list with sexp;;
 
 (** Read all the basis functions of an element *)
 let read in_channel at_number =
@@ -36,8 +36,28 @@ let read_element in_channel at_number element =
 ;;
 
 let to_string b =
-  List.map ~f:(fun (g,n) ->
+  let new_nucleus n = 
-     let n = Nucl_number.to_int n in
+    Printf.sprintf "Atom %d" n
-     (Int.to_string n)^":"^(Gto.to_string g)) b
+  in
+  
+  let rec do_work accu current_nucleus = function
+  | [] -> List.rev accu
+  | (g,n)::tail -> 
+    let n = Nucl_number.to_int n
+    in
+    let accu = 
+       if (n <> current_nucleus) then
+         (new_nucleus n)::""::accu
+       else
+         accu
+    in
+    do_work ((Gto.to_string g)::accu) n tail
+  in
+  do_work [new_nucleus 1] 1 b
   |> String.concat ~sep:"\n"
 ;;
+
+include To_md5;;
+let to_md5 = to_md5 sexp_of_t
+;;
+

ocaml/Basis.mli

@@ -15,3 +15,6 @@ val read_element :
 
 (** Convert the basis to a string *)
 val to_string :  (Gto.t * Nucl_number.t) list -> string
+
+(** Convert the basis to an MD5 hash *)
+val to_md5 :  (Gto.t * Nucl_number.t) list -> MD5.t

ocaml/Bit.ml

@@ -1,3 +1,4 @@
+open Core.Std;;
 
 (*
 Type for bits 
@@ -7,9 +8,10 @@ Zero | One
 
 *)
 
-type bit      = 
+type t      = 
-   | One 
+| One 
-   | Zero 
+| Zero 
+with sexp
 
 let to_string = function
   | Zero -> "0"

ocaml/Bit.mli

@@ -0,0 +1,10 @@
+type t = One | Zero with sexp
+
+(** String conversions for printing *)
+val to_string : t -> string
+
+(** Logical operations *)
+val and_operator : t -> t -> t
+val or_operator  : t -> t -> t
+val xor_operator : t -> t -> t
+val not_operator : t -> t

ocaml/Bitlist.ml

@@ -1,4 +1,5 @@
 open Qptypes;;
+open Core.Std;;
 
 (*
 Type for bits strings
@@ -7,7 +8,7 @@ Type for bits strings
 list of Bits
 *)
 
-type bit_list = Bit.bit list
+type t = Bit.t list
 
 (* String representation *)
 let to_string b = 
@@ -20,6 +21,13 @@ let to_string b =
   do_work "" b
 ;;
 
+let of_string ?(zero='0') ?(one='1') s =
+  String.to_list s
+  |> List.rev_map ~f:( fun c ->
+     if (c = zero) then Bit.Zero
+     else if (c = one) then Bit.One
+     else (failwith ("Error in string "^s) ) )
+;;
 
 (* Create a bit list from an int64 *)
 let of_int64 i = 
@@ -27,7 +35,7 @@ let of_int64 i =
   | 0L -> [ Bit.Zero ]
   | 1L -> [ Bit.One ]
   | i -> let b =
-    match (Int64.logand i 1L ) with
+    match (Int64.bit_and i 1L ) with
     | 0L -> Bit.Zero
     | 1L -> Bit.One 
     | _ -> raise (Failure "i land 1 not in (0,1)")  
@@ -51,15 +59,14 @@ let to_int64 l =
   let rec do_work accu = function
     | [] -> accu
     | Bit.Zero::tail -> do_work Int64.(shift_left accu 1) tail
-    | Bit.One::tail  -> do_work Int64.(logor one (shift_left accu 1)) tail
+    | Bit.One::tail  -> do_work Int64.(bit_or one (shift_left accu 1)) tail
   in do_work Int64.zero (List.rev l)
 ;;
 
 (* Create a bit list from a list of int64 *)
 let of_int64_list l = 
-  let list_of_lists = List.map of_int64 l in
+  List.map ~f:of_int64 l 
-  let result = List.rev list_of_lists in
+  |> List.concat
-  List.flatten result
 ;;
 
 (* Compute n_int *)
@@ -92,27 +99,28 @@ let to_int64_list l =
   in
   let l = do_work [] [] 1 l
   in
-  List.rev_map to_int64 l
+  List.rev_map ~f:to_int64 l
 ;;
 
 (* Create a bit list from a list of MO indices *)
 let of_mo_number_list n_int l = 
   let n_int = N_int_number.to_int n_int in
   let length = n_int*64 in
-  let a = Array.make length (Bit.Zero) in
+  let a = Array.create length (Bit.Zero) in
-  List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l;
+  List.iter ~f:(fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l;
   Array.to_list a
 ;;
 
 let to_mo_number_list l =
   let a = Array.of_list l in
+  let mo_tot_num = MO_number.get_max () in
   let rec do_work accu = function
   | 0 -> accu
   | i ->
       begin
         let new_accu = 
         match a.(i-1) with
-        | Bit.One  -> (MO_number.of_int i)::accu 
+        | Bit.One  -> (MO_number.of_int ~max:mo_tot_num i)::accu 
         | Bit.Zero -> accu 
         in
         do_work new_accu (i-1)
@@ -161,30 +169,4 @@ let popcnt b =
   in popcnt 0 b
 ;;
 
-let test_module () = 
-  let test = of_int64_list ([-1231L;255L]) in
-  print_string (to_string test);
-  print_newline ();
-  print_string (string_of_int (String.length (to_string test)));
-  print_newline ();
-  print_string ( Bit.to_string Bit.One );
-
-  let a = of_int64_list ([-1L;0L]) 
-  and b = of_int64_list ([128L;127L])
-  in begin
-   print_newline ();
-   print_newline ();
-   print_string (to_string a);
-   print_newline ();
-   print_string (to_string b);
-   print_newline ();
-   print_string (to_string (and_operator a b));
-   print_newline ();
-   print_string (to_string (or_operator a b));
-   print_newline ();
-   print_string (to_string (xor_operator a b));
-   print_string (to_string a);
-   print_int (popcnt a);
-  end
-;;
 

ocaml/Bitlist.mli

@@ -0,0 +1,35 @@
+type t = Bit.t list
+
+(** The zero bit list *)
+val zero : Qptypes.N_int_number.t -> t
+
+(** Convert to a string for printing *)
+val to_string : t -> string
+
+(** Convert to a string for printing *)
+val of_string : ?zero:char -> ?one:char -> string -> t
+
+(** int64 conversion functions *)
+
+val of_int64 : int64 -> t
+val to_int64 : t -> int64
+
+val of_int64_list : int64 list -> t
+val to_int64_list : t -> int64 list
+
+(** Get the number of needed int64 elements to encode the bit list *)
+val n_int_of_mo_tot_num : int -> Qptypes.N_int_number.t
+
+(** Conversion to MO numbers *)
+val to_mo_number_list : t -> Qptypes.MO_number.t list
+val of_mo_number_list :
+  Qptypes.N_int_number.t -> Qptypes.MO_number.t list -> t
+
+(** Logical operators *)
+val and_operator : t -> t -> t
+val xor_operator : t -> t -> t
+val or_operator  : t -> t -> t
+val not_operator : t -> t
+
+(** Count the number of bits set to one *)
+val popcnt : t -> int

ocaml/Charge.ml

@@ -1,11 +1,12 @@
 open Core.Std;;
 
-type t = float
+type t = float with sexp;;
 
 let of_float x = x
 let of_int   i = Float.of_int i
 let of_string s = Float.of_string s
 
+
 let to_float x = x
 let to_int   x = Float.to_int x
 let to_string x = 

ocaml/Charge.mli

@@ -1,8 +1,13 @@
-type t = float
+type t = float with sexp
 
+(** Float conversion functions *)
 val to_float : t -> float
-val to_int   : t -> int
-val to_string: t -> string 
 val of_float : float -> t
+
+(** Int conversion functions *)
+val to_int   : t -> int
 val of_int   : int -> t
+
+(** String conversion functions *)
+val to_string: t -> string 
 val of_string: string -> t

ocaml/Determinant.ml

@@ -0,0 +1,73 @@
+open Core.Std;;
+open Qptypes;;
+
+type t = int64 array with sexp
+
+let to_int64_array (x:t) = (x:int64 array)
+;;
+
+let to_alpha_beta x = 
+  let x = to_int64_array x in
+  let n_int = (Array.length x)/2 in
+  ( Array.init n_int ~f:(fun i -> x.(i)) ,
+    Array.init n_int ~f:(fun i -> x.(i+n_int)) )
+;;
+
+let to_bitlist_couple x =
+  let (xa,xb) = to_alpha_beta x in
+  let xa = to_int64_array xa
+  |> Array.to_list
+  |> Bitlist.of_int64_list
+  and xb = to_int64_array xb
+  |> Array.to_list
+  |> Bitlist.of_int64_list
+  in (xa,xb)
+;;
+
+let bitlist_to_string ~mo_tot_num x =
+  List.map x ~f:(fun i -> match i with
+    | Bit.Zero -> "-"
+    | Bit.One  -> "+" )
+  |> String.concat
+  |> String.sub ~pos:0 ~len:(MO_number.to_int mo_tot_num) 
+;;
+
+
+let of_int64_array ~n_int ~alpha ~beta x =
+   assert ((Array.length x) = (N_int_number.to_int n_int)*2) ;
+   let (a,b) = to_bitlist_couple x 
+   and alpha = Elec_alpha_number.to_int alpha
+   and beta  = Elec_beta_number.to_int beta
+   in
+   if ( (Bitlist.popcnt a) <> alpha) then
+     begin
+       let mo_tot_num = MO_number.get_max () in
+       let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num  in
+       failwith (Printf.sprintf "Expected %d electrons in alpha determinant
+%s" alpha (bitlist_to_string ~mo_tot_num:mo_tot_num a) )
+     end;
+   if ( (Bitlist.popcnt b) <> beta ) then
+     begin
+       let mo_tot_num = MO_number.get_max () in
+       let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num  in
+       failwith (Printf.sprintf "Expected %d electrons in beta determinant
+%s" beta (bitlist_to_string ~mo_tot_num:mo_tot_num b) )
+     end;
+   x
+;;
+
+let of_bitlist_couple ~alpha ~beta (xa,xb) =
+  let ba = Bitlist.to_int64_list xa in
+  let bb = Bitlist.to_int64_list xb in
+  let n_int = Bitlist.n_int_of_mo_tot_num (List.length xa) in
+  of_int64_array ~n_int:n_int ~alpha:alpha ~beta:beta (Array.of_list (ba@bb))
+;;
+
+let to_string ~mo_tot_num x = 
+  let (xa,xb) = to_bitlist_couple x in
+  [ bitlist_to_string ~mo_tot_num:mo_tot_num xa ;
+    bitlist_to_string ~mo_tot_num:mo_tot_num xb ]
+  |> String.concat ~sep:"\n"
+;;
+
+

ocaml/Determinant.mli

@@ -0,0 +1,32 @@
+(** Determinants are stored as follows :
+  *    <-------- N_int ---------->
+  * [|  i1_alpha ; i2_alpha ; ... ;
+  *     i1_beta  ; i2_beta  ; ... ; |]
+  * 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
+
+(** Transform to an int64 array *)
+val to_int64_array : t -> int64 array
+
+(** Create from an int64 array, checking the number of alpha 
+  * and beta electrons *)
+val of_int64_array : n_int:Qptypes.N_int_number.t -> 
+  alpha:Qptypes.Elec_alpha_number.t ->
+  beta:Qptypes.Elec_beta_number.t -> 
+  int64 array -> t
+
+(** Split into an alpha-only and a beta-only determinant *)
+val to_alpha_beta : t -> (int64 array)*(int64 array)
+
+(** Transform to a bit list *)
+val to_bitlist_couple : t -> Bitlist.t * Bitlist.t
+
+(** Create from a bit list *)
+val of_bitlist_couple : alpha:Qptypes.Elec_alpha_number.t ->
+  beta:Qptypes.Elec_beta_number.t -> 
+  Bitlist.t * Bitlist.t -> t
+
+(** String representation *)
+val to_string : mo_tot_num:Qptypes.MO_number.t -> t -> string

ocaml/Element.ml

@@ -8,7 +8,7 @@ type t =
 |Li|Be                              |B |C |N |O |F |Ne
 |Na|Mg                              |Al|Si|P |S |Cl|Ar
 |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr
-;;
+with sexp;;
 
 let of_string x = 
   match (String.capitalize (String.lowercase x)) with
@@ -132,47 +132,49 @@ let to_long_string = function
 | Kr  -> "Krypton"
 ;;
 
-let to_charge = function
+let to_charge c = 
-| X   -> 0
+  let result = match c with
-| H   -> 1
+  | X   -> 0
-| He  -> 2
+  | H   -> 1
-| Li  -> 3
+  | He  -> 2
-| Be  -> 4
+  | Li  -> 3
-| B   -> 5
+  | Be  -> 4
-| C   -> 6
+  | B   -> 5
-| N   -> 7
+  | C   -> 6
-| O   -> 8
+  | N   -> 7
-| F   -> 9
+  | O   -> 8
-| Ne  -> 10
+  | F   -> 9
-| Na  -> 11
+  | Ne  -> 10
-| Mg  -> 12
+  | Na  -> 11
-| Al  -> 13
+  | Mg  -> 12
-| Si  -> 14
+  | Al  -> 13
-| P   -> 15
+  | Si  -> 14
-| S   -> 16
+  | P   -> 15
-| Cl  -> 17
+  | S   -> 16
-| Ar  -> 18
+  | Cl  -> 17
-| K   -> 19
+  | Ar  -> 18
-| Ca  -> 20
+  | K   -> 19
-| Sc  -> 21
+  | Ca  -> 20
-| Ti  -> 22
+  | Sc  -> 21
-| V   -> 23
+  | Ti  -> 22
-| Cr  -> 24
+  | V   -> 23
-| Mn  -> 25
+  | Cr  -> 24
-| Fe  -> 26
+  | Mn  -> 25
-| Co  -> 27
+  | Fe  -> 26
-| Ni  -> 28
+  | Co  -> 27
-| Cu  -> 29
+  | Ni  -> 28
-| Zn  -> 30
+  | Cu  -> 29
-| Ga  -> 31
+  | Zn  -> 30
-| Ge  -> 32
+  | Ga  -> 31
-| As  -> 33
+  | Ge  -> 32
-| Se  -> 34
+  | As  -> 33
-| Br  -> 35
+  | Se  -> 34
-| Kr  -> 36
+  | Br  -> 35
+  | Kr  -> 36
+  in Charge.of_int result
 ;;
 
-let of_charge = function
+let of_charge c = match (Charge.to_int c) with
 |  0   ->  X   
 |  1   ->  H   
 |  2   ->  He  

ocaml/Element.mli

@@ -0,0 +1,18 @@
+exception ElementError of string
+
+type t =
+|X
+|H                                                 |He
+|Li|Be                              |B |C |N |O |F |Ne
+|Na|Mg                              |Al|Si|P |S |Cl|Ar
+|K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr
+with sexp
+
+(** String conversion functions *)
+val of_string : string -> t
+val to_string : t -> string
+val to_long_string : t -> string
+
+(** get the positive charge *)
+val to_charge : t -> Charge.t
+val of_charge : Charge.t -> t

ocaml/Excitation.ml

@@ -1,22 +1,14 @@
 open Core.Std;;
 open Qptypes;;
 
-module Hole : sig
+module Hole = struct
-    type t
+    type t = MO_class.t with sexp
-    val to_mo_class  : t -> MO_class.t
-    val of_mo_class  : MO_class.t -> t
-end = struct
-    type t = MO_class.t
     let of_mo_class x = x
     let to_mo_class x = x
 end
 
-module Particle : sig
+module Particle = struct
-    type t
+    type t = MO_class.t with sexp
-    val to_mo_class  : t -> MO_class.t
-    val of_mo_class  : MO_class.t -> t
-end = struct
-    type t = MO_class.t
     let of_mo_class x = x
     let to_mo_class x = x
 end
@@ -24,10 +16,7 @@ end
 type t =
 | Single of Hole.t*Particle.t
 | Double of Hole.t*Particle.t*Hole.t*Particle.t
-;;
+with sexp;;
-
-let failwith s = raise (Failure s)
-;;
 
 let create_single ~hole ~particle =
   MO_class.(

ocaml/Excitation.mli

@@ -0,0 +1,30 @@
+module Hole :
+  sig
+    type t
+    val to_mo_class : t -> MO_class.t
+    val of_mo_class : MO_class.t -> t
+    val t_of_sexp : Sexplib.Sexp.t -> t
+    val sexp_of_t : t -> Sexplib.Sexp.t
+  end
+module Particle :
+  sig
+    type t
+    val to_mo_class : t -> MO_class.t
+    val of_mo_class : MO_class.t -> t
+    val t_of_sexp : Sexplib.Sexp.t -> t
+    val sexp_of_t : t -> Sexplib.Sexp.t
+  end
+
+type t =
+  | Single of Hole.t * Particle.t
+  | Double of Hole.t * Particle.t * Hole.t * Particle.t
+with sexp
+
+val create_single : hole:MO_class.t -> particle:MO_class.t -> t
+
+val double_of_singles : t -> t -> t
+
+val create_double : hole1:MO_class.t -> particle1:MO_class.t ->
+                    hole2:MO_class.t -> particle2:MO_class.t -> t
+
+val to_string : t -> string

ocaml/Gto.ml

@@ -7,7 +7,7 @@ exception End_Of_Basis
 type t =
 { sym  : Symmetry.t ;
   lc   : ((Primitive.t * AO_coef.t) list)
-}
+} with sexp
 ;;
 
 let of_prim_coef_list pc =
@@ -52,24 +52,37 @@ let read_one in_channel =
       match buffer with
       | [ j ; expo ; coef ] ->
         begin
-          let p = { Primitive.sym = sym ;
+          let p =
-                    Primitive.expo = AO_expo.of_float
+            Primitive.of_sym_expo sym 
-                        (Float.of_string expo)
+              (AO_expo.of_float (Float.of_string expo) )
-                  }
           and c = AO_coef.of_float (Float.of_string coef) in
           read_lines ( (p,c)::result) (i-1)
         end
       | _ -> raise (GTO_Read_Failure line_buffer)
     end
   in read_lines [] n
+  |> List.rev
   |> of_prim_coef_list
 ;;
 
 
 (** Transform the gto to a string *)
 let to_string { sym = sym ; lc = lc } =
-  let f (p,c) = Printf.sprintf "( %s, %f )" (Primitive.to_string p) (AO_coef.to_float c)
+  let result = 
+    Printf.sprintf "%s %3d" (Symmetry.to_string sym) (List.length lc)
   in
-  Printf.sprintf "[ %s : %s ]" (Symmetry.to_string sym)
+  let rec do_work accu i = function
-    (String.concat (List.map ~f:f lc) ~sep:", ")
+  | [] -> List.rev accu
+  | (p,c)::tail -> 
+    let p = AO_expo.to_float p.Primitive.expo
+    and c = AO_coef.to_float c
+    in
+    let result = 
+      Printf.sprintf "%3d %16f  %16f" i p c
+    in
+    do_work (result::accu) (i+1) tail
+  in
+  (do_work [result] 1 lc)
+  |> String.concat ~sep:"\n"
 ;;
+

ocaml/Gto.mli

@@ -0,0 +1,16 @@
+exception GTO_Read_Failure of string
+exception End_Of_Basis
+type t =
+  { sym : Symmetry.t ;
+    lc  : (Primitive.t * Qptypes.AO_coef.t) list;
+  } with sexp
+
+(** Create from a list of Primitive.t * Qptypes.AO_coef.t *)
+val of_prim_coef_list :
+  (Primitive.t * Qptypes.AO_coef.t) list -> t
+
+(** Read from a file *)
+val read_one : in_channel -> t
+
+(** Convert to string for printing *)
+val to_string : t -> string

ocaml/Input_ao_basis.ml

@@ -12,10 +12,12 @@ module Ao_basis : sig
       ao_power        : Symmetry.Xyz.t array;
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
-    }
+    } with sexp
   ;;
   val read : unit -> t
   val to_string : t -> string
+  val to_md5 : t -> MD5.t
+  val to_rst : t -> Rst_string.t
 end = struct
   type t = 
     { ao_basis        : string ;
@@ -26,7 +28,7 @@ end = struct
       ao_power        : Symmetry.Xyz.t array;
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "ao_basis";;
@@ -44,28 +46,29 @@ end = struct
   ;;
 
   let read_ao_prim_num () =
-    (Ezfio.get_ao_basis_ao_prim_num () ).Ezfio.data
+    Ezfio.get_ao_basis_ao_prim_num () 
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.map ~f:AO_prim_number.of_int
   ;;
 
   let read_ao_prim_num_max () =
-    (Ezfio.get_ao_basis_ao_prim_num () ).Ezfio.data
+    Ezfio.get_ao_basis_ao_prim_num ()
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.fold ~f:(fun x y -> if x>y then x else y) ~init:0
     |> AO_prim_number.of_int
   ;;
 
   let read_ao_nucl () =
-    (Ezfio.get_ao_basis_ao_nucl () ).Ezfio.data
+    let nmax = Nucl_number.get_max () in
-    |> Ezfio.flattened_ezfio_data
+    Ezfio.get_ao_basis_ao_nucl () 
-    |> Array.map ~f:Nucl_number.of_int
+    |> Ezfio.flattened_ezfio
+    |> Array.map ~f:(fun x-> Nucl_number.of_int ~max:nmax x)
   ;;
 
   let read_ao_power () =
     let x = Ezfio.get_ao_basis_ao_power () in
     let dim = x.Ezfio.dim.(0) in
-    let data = Ezfio.flattened_ezfio_data x.Ezfio.data in
+    let data = Ezfio.flattened_ezfio x in
     let result = Array.init dim ~f:(fun x -> "") in
     for i=1 to dim
     do
@@ -80,14 +83,14 @@ end = struct
   ;;
 
   let read_ao_coef () =
-    (Ezfio.get_ao_basis_ao_coef () ).Ezfio.data
+    Ezfio.get_ao_basis_ao_coef () 
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.map ~f:AO_coef.of_float
   ;;
 
   let read_ao_expo () =
-    (Ezfio.get_ao_basis_ao_expo () ).Ezfio.data
+    Ezfio.get_ao_basis_ao_expo ()
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.map ~f:AO_expo.of_float
   ;;
 
@@ -103,9 +106,104 @@ end = struct
     }
   ;;
     
+  let to_long_basis b =
+    let ao_num = AO_number.to_int b.ao_num in
+    let gto_array = Array.init (AO_number.to_int b.ao_num)
+      ~f:(fun i ->
+        let s = Symmetry.Xyz.to_symmetry b.ao_power.(i) in
+        let ao_prim_num = AO_prim_number.to_int b.ao_prim_num.(i) in
+        let prims = List.init ao_prim_num ~f:(fun j ->
+          let prim = { Primitive.sym  = s ;
+                       Primitive.expo = b.ao_expo.(ao_num*j+i)
+                     }
+          in
+          let coef = b.ao_coef.(ao_num*j+i) in
+          (prim,coef)
+        ) in
+        Gto.of_prim_coef_list prims
+      )
+    in
+    let rec do_work accu sym gto nucl =
+      match (sym, gto, nucl) with
+        | (s::srest, g::grest, n::nrest) -> 
+          do_work ((s,g,n)::accu) srest grest nrest
+        | ([],[],[]) -> List.rev accu
+        | _ -> assert false
+    in
+    do_work [] 
+      (Array.to_list b.ao_power)
+      (Array.to_list gto_array)
+      (Array.to_list b.ao_nucl)
+  ;;
+
+  let to_basis b =
+    to_long_basis b
+    |> Long_basis.to_basis
+  ;;
+      
+  let to_rst b =
+    let print_sym = 
+      let l = List.init (Array.length b.ao_power) ~f:(
+         fun i -> ( (i+1),b.ao_nucl.(i),b.ao_power.(i) ) ) in
+      let rec do_work = function
+      | [] -> []
+      | (i,n,x)::tail  -> 
+          (Printf.sprintf " %5d  %6d     %-8s\n" i (Nucl_number.to_int n) (Symmetry.Xyz.to_string x))::
+          (do_work tail)
+      in do_work l
+      |> String.concat
+    in
+
+    let short_basis = to_basis b in
+    Printf.sprintf "
+Name of the AO basis ::
+
+  ao_basis = %s
+
+Basis set ::
+  
+%s
+
+
+======= ========= ===========
+ Basis   Nucleus   Symmetries
+======= ========= ===========
+%s
+======= ========= ===========
+
+" b.ao_basis
+    (Basis.to_string short_basis 
+     |> String.split ~on:'\n'
+     |> List.map ~f:(fun x-> "  "^x)
+     |> String.concat ~sep:"\n"
+    ) print_sym
+  |> Rst_string.of_string
+  ;;
+
+  let read_rst s = 
+    let s = Rst_string.to_string s
+    |> String.split ~on:'\n'
+    in
+    let rec extract_basis = function
+    | [] -> failwith "Error in basis set"
+    | line :: tail ->
+      let line = String.strip line in
+      if line = "Basis set ::" then
+        String.concat tail ~sep:"\n"
+      else
+        extract_basis tail
+    in
+    extract_basis s
+  ;;
+
+  let to_md5 b =
+    let short_basis = to_basis b in
+    Basis.to_md5 short_basis
+  ;;
+    
   let to_string b =
     Printf.sprintf "
-ao_basis        = \"%s\"
+ao_basis        = %s
 ao_num          = %s
 ao_prim_num     = %s
 ao_prim_num_max = %s
@@ -113,6 +211,7 @@ ao_nucl         = %s
 ao_power        = %s
 ao_coef         = %s
 ao_expo         = %s
+md5             = %s
 "
     b.ao_basis
     (AO_number.to_string b.ao_num)
@@ -127,7 +226,8 @@ ao_expo         = %s
       |> String.concat ~sep:", ")
     (b.ao_expo  |> Array.to_list |> List.map ~f:AO_expo.to_string
       |> String.concat ~sep:", ")
+    (to_md5 b |> MD5.to_string )
 
+  ;;
 end
 
-

ocaml/Input_bi_integrals.ml

@@ -11,10 +11,13 @@ module Bielec_integrals : sig
       threshold_ao       : Threshold.t;
       threshold_mo       : Threshold.t;
       direct             : bool;
-    }
+    } with sexp
   ;;
-  val read : unit -> t
+  val read  : unit -> t
+  val write : t -> unit
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
+  val of_rst : Rst_string.t -> t 
 end = struct
   type t = 
     { read_ao_integrals  : bool;
@@ -24,7 +27,7 @@ end = struct
       threshold_ao       : Threshold.t;
       threshold_mo       : Threshold.t;
       direct             : bool;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "bielec_integrals";;
@@ -38,6 +41,11 @@ end = struct
     Ezfio.get_bielec_integrals_read_ao_integrals ()
   ;;
 
+  let write_read_ao_integrals =
+    Ezfio.set_bielec_integrals_read_ao_integrals
+  ;;
+
+
   let read_read_mo_integrals () = 
     if not (Ezfio.has_bielec_integrals_read_mo_integrals ()) then
        get_default "read_mo_integrals"
@@ -47,6 +55,11 @@ end = struct
     Ezfio.get_bielec_integrals_read_mo_integrals ()
   ;;
 
+  let write_read_mo_integrals =
+    Ezfio.set_bielec_integrals_read_mo_integrals
+  ;;
+
+
   let read_write_ao_integrals () = 
     if not (Ezfio.has_bielec_integrals_write_ao_integrals ()) then
        get_default "write_ao_integrals"
@@ -56,6 +69,11 @@ end = struct
     Ezfio.get_bielec_integrals_write_ao_integrals ()
   ;;
 
+  let write_write_ao_integrals =
+    Ezfio.set_bielec_integrals_write_ao_integrals
+  ;;
+
+
   let read_write_mo_integrals () = 
     if not (Ezfio.has_bielec_integrals_write_mo_integrals ()) then
        get_default "write_mo_integrals"
@@ -65,6 +83,11 @@ end = struct
     Ezfio.get_bielec_integrals_write_mo_integrals ()
   ;;
 
+  let write_write_mo_integrals =
+    Ezfio.set_bielec_integrals_write_mo_integrals
+  ;;
+
+
   let read_direct () = 
     if not (Ezfio.has_bielec_integrals_direct ()) then
        get_default "direct"
@@ -74,6 +97,11 @@ end = struct
     Ezfio.get_bielec_integrals_direct ()
   ;;
 
+  let write_direct =
+    Ezfio.set_bielec_integrals_direct
+  ;;
+
+
   let read_threshold_ao () = 
     if not (Ezfio.has_bielec_integrals_threshold_ao ()) then
        get_default "threshold_ao"
@@ -84,6 +112,12 @@ end = struct
     |> Threshold.of_float
   ;;
 
+  let write_threshold_ao t =
+    Threshold.to_float t
+    |> Ezfio.set_bielec_integrals_threshold_ao
+  ;;
+
+
   let read_threshold_mo () = 
     if not (Ezfio.has_bielec_integrals_threshold_mo ()) then
        get_default "threshold_mo"
@@ -94,6 +128,12 @@ end = struct
     |> Threshold.of_float
   ;;
 
+  let write_threshold_mo t =
+    Threshold.to_float t
+    |> Ezfio.set_bielec_integrals_threshold_mo
+  ;;
+
+
   let read ()= 
     let result = 
     { read_ao_integrals  = read_read_ao_integrals();
@@ -113,15 +153,31 @@ end = struct
     result
   ;;
 
+  let write b =
+    if (b.read_ao_integrals &&
+        b.write_ao_integrals) then
+          failwith "Read and Write AO integrals are both true.";
+    if (b.read_mo_integrals &&
+        b.write_mo_integrals) then
+          failwith "Read and Write MO integrals are both true.";
+    write_read_ao_integrals  b.read_ao_integrals;
+    write_read_mo_integrals  b.read_mo_integrals;
+    write_write_ao_integrals b.write_ao_integrals ;
+    write_write_mo_integrals b.write_mo_integrals ;
+    write_threshold_ao       b.threshold_ao;
+    write_threshold_mo       b.threshold_mo;
+    write_direct             b.direct;
+  ;;
+
   let to_string b =
     Printf.sprintf "
-read_ao_integrals  = %s
+read_ao_integrals = %s
-read_mo_integrals  = %s
+read_mo_integrals = %s
 write_ao_integrals = %s
 write_mo_integrals = %s
-threshold_ao       = %s
+threshold_ao = %s
-threshold_mo       = %s
+threshold_mo = %s
-direct             = %s
+direct = %s
 "
         (Bool.to_string b.read_ao_integrals)
         (Bool.to_string b.read_mo_integrals)
@@ -130,6 +186,55 @@ direct             = %s
         (Threshold.to_string b.threshold_ao)
         (Threshold.to_string b.threshold_mo)
         (Bool.to_string b.direct)
+  ;;
+
+  let to_rst b =
+    Printf.sprintf "
+Read AO/MO integrals from disk ::
+
+  read_ao_integrals = %s
+  read_mo_integrals = %s
+
+Write AO/MO integrals to disk ::
+
+  write_ao_integrals = %s
+  write_mo_integrals = %s
+
+Thresholds on integrals ::
+
+  threshold_ao = %s
+  threshold_mo = %s
+
+Direct calculation of integrals ::
+
+  direct = %s
+
+"
+        (Bool.to_string b.read_ao_integrals)
+        (Bool.to_string b.read_mo_integrals)
+        (Bool.to_string b.write_ao_integrals)
+        (Bool.to_string b.write_mo_integrals)
+        (Threshold.to_string b.threshold_ao)
+        (Threshold.to_string b.threshold_mo)
+        (Bool.to_string b.direct)
+  |> Rst_string.of_string
+  ;;
+
+  let of_rst s =
+    let s = Rst_string.to_string s
+    |> String.split ~on:'\n'
+    |> List.filter ~f:(fun line ->
+        String.contains line '=')
+    |> List.map ~f:(fun line ->
+        "("^(
+        String.tr line ~target:'=' ~replacement:' '
+        )^")" )
+    |> String.concat
+    in
+    Sexp.of_string ("("^s^")")
+    |> t_of_sexp
+  ;;
+
 end
 
 

ocaml/Input_bitmasks.ml

@@ -8,7 +8,7 @@ module Bitmasks : sig
       bit_kind           : Bit_kind.t;
       n_mask_gen         : Bitmask_number.t;
       generators         : int64 array;
-    }
+    } with sexp
   ;;
   val read : unit -> t
   val to_string : t -> string
@@ -18,7 +18,7 @@ end = struct
       bit_kind           : Bit_kind.t;
       n_mask_gen         : Bitmask_number.t;
       generators         : int64 array;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "bitmasks";;
@@ -72,8 +72,8 @@ end = struct
       in
       Ezfio.set_bitmasks_generators generators
       end;
-    (Ezfio.get_bitmasks_generators ()).Ezfio.data
+    Ezfio.get_bitmasks_generators ()
-    |> Ezfio.flattened_ezfio_data 
+    |> Ezfio.flattened_ezfio
   ;;
 
   let read () = 

ocaml/Input_cis.ml

@@ -10,10 +10,11 @@ module Cis_dressed : sig
       mp2_dressing       : bool;
       standard_doubles   : bool;
       en_2_2             : bool;
-    }
+    } with sexp
   ;;
   val read : unit -> t
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
 end = struct
   type t = 
     { n_state_cis        : States_number.t;
@@ -22,7 +23,7 @@ end = struct
       mp2_dressing       : bool;
       standard_doubles   : bool;
       en_2_2             : bool;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "cis_dressed";;
@@ -95,9 +96,9 @@ end = struct
 
   let to_string b =
     Printf.sprintf "
-n_state_cis        = %s
+n_state_cis = %s
-n_core_cis         = %s
+n_core_cis = %s
-n_act_cis          = %s
+n_act_cis  = %s
 mp2_dressing       = %s
 standard_doubles   = %s
 en_2_2             = %s
@@ -108,6 +109,42 @@ en_2_2             = %s
         (Bool.to_string b.mp2_dressing)
         (Bool.to_string b.standard_doubles)
         (Bool.to_string b.en_2_2)
+  ;;
+
+
+  let to_rst b =
+    Printf.sprintf "
+Number of states ::
+
+  n_state_cis = %s
+
+Core and active MOs ::
+
+  n_core_cis = %s
+  n_act_cis  = %s
+
+Dress with MP2 perturbation ::
+
+  mp2_dressing       = %s
+
+Use standard double-excitations ::
+
+  standard_doubles   = %s
+
+Epstein-Nesbet 2x2 diagonalization ::
+
+  en_2_2             = %s
+
+"
+        (States_number.to_string b.n_state_cis)
+        (Positive_int.to_string b.n_core_cis)
+        (Positive_int.to_string b.n_act_cis)
+        (Bool.to_string b.mp2_dressing)
+        (Bool.to_string b.standard_doubles)
+        (Bool.to_string b.en_2_2)
+
+   |> Rst_string.of_string
+  ;;
 end
 
 

ocaml/Input_cisd_sc2.ml

@@ -4,19 +4,22 @@ open Core.Std;;
 
 module Cisd_sc2 : sig
   type t = 
-    { n_det_max_cisd_sc2 : Det_number.t;
+    { n_det_max_cisd_sc2 : Det_number_max.t;
       pt2_max            : PT2_energy.t;
       do_pt2_end         : bool;
-    }
+    } with sexp
   ;;
-  val read : unit -> t
+  val read  : unit -> t
+  val write : t -> unit
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
+  val of_rst : Rst_string.t -> t
 end = struct
   type t = 
-    { n_det_max_cisd_sc2 : Det_number.t;
+    { n_det_max_cisd_sc2 : Det_number_max.t;
       pt2_max            : PT2_energy.t;
       do_pt2_end         : bool;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "cisd_sc2_selected";;
@@ -28,7 +31,12 @@ end = struct
        |> Ezfio.set_cisd_sc2_selected_n_det_max_cisd_sc2 
     ;
     Ezfio.get_cisd_sc2_selected_n_det_max_cisd_sc2 ()
-    |> Det_number.of_int
+    |> Det_number_max.of_int
+  ;;
+
+  let write_n_det_max_cisd_sc2 n =
+    Det_number_max.to_int n
+    |> Ezfio.set_cisd_sc2_selected_n_det_max_cisd_sc2
   ;;
 
 
@@ -42,6 +50,11 @@ end = struct
     |> PT2_energy.of_float
   ;;
 
+  let write_pt2_max p =
+    PT2_energy.to_float p
+    |> Ezfio.set_cisd_sc2_selected_pt2_max
+  ;;
+
 
   let read_do_pt2_end () = 
     if not (Ezfio.has_cisd_sc2_selected_do_pt2_end ()) then
@@ -52,6 +65,10 @@ end = struct
     Ezfio.get_cisd_sc2_selected_do_pt2_end ()
   ;;
 
+  let write_do_pt2_end =
+    Ezfio.set_cisd_sc2_selected_do_pt2_end
+  ;;
+
 
   let read () = 
     { n_det_max_cisd_sc2  = read_n_det_max_cisd_sc2 ();
@@ -60,15 +77,64 @@ end = struct
     }
   ;;
 
+  let write { n_det_max_cisd_sc2 ;
+              pt2_max            ;
+              do_pt2_end         ;
+            } =
+     write_n_det_max_cisd_sc2 n_det_max_cisd_sc2;
+     write_pt2_max            pt2_max;
+     write_do_pt2_end         do_pt2_end;
+  ;;
+
   let to_string b =
     Printf.sprintf "
 n_det_max_cisd_sc2 = %s
-pt2_max            = %s
+pt2_max = %s
-do_pt2_end         = %s
+do_pt2_end = %s
 "
-        (Det_number.to_string b.n_det_max_cisd_sc2)
+        (Det_number_max.to_string b.n_det_max_cisd_sc2)
         (PT2_energy.to_string b.pt2_max)
         (Bool.to_string b.do_pt2_end)
+  ;;
+
+  let to_rst b =
+    Printf.sprintf "
+Stop when the `n_det` > `n_det_max_cisd_sc2` ::
+
+  n_det_max_cisd_sc2 = %s
+
+Stop when -E(PT2) < `pt2_max` ::
+
+  pt2_max = %s
+
+Compute E(PT2) at the end ::
+
+  do_pt2_end = %s
+
+"
+        (Det_number_max.to_string b.n_det_max_cisd_sc2)
+        (PT2_energy.to_string b.pt2_max)
+        (Bool.to_string b.do_pt2_end)
+   |> Rst_string.of_string
+  ;;
+
+  let of_rst s =
+    let s = Rst_string.to_string s
+    |> String.split ~on:'\n'
+    |> List.filter ~f:(fun line ->
+        String.contains line '=')
+    |> List.map ~f:(fun line ->
+        "("^(
+        String.tr line ~target:'=' ~replacement:' '
+        )^")" )
+    |> String.concat
+    in
+    Sexp.of_string ("("^s^")")
+    |> t_of_sexp
+  ;;
+
+
+
 end
 
 

ocaml/Input_determinants.ml

@@ -6,11 +6,11 @@ module Determinants : sig
   type t = 
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
-      mo_label               : Non_empty_string.t;
+      mo_label               : MO_label.t;
       n_det                  : Det_number.t;
       n_states               : States_number.t;
       n_states_diag          : States_number.t;
-      n_det_max_jacobi       : Det_number.t;
+      n_det_max_jacobi       : Strictly_positive_int.t;
       threshold_generators   : Threshold.t;
       threshold_selectors    : Threshold.t; 
       read_wf                : bool;
@@ -18,19 +18,21 @@ module Determinants : sig
       s2_eig                 : bool;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
-    }
+    } with sexp
-  ;;
+  val read  : unit -> t
-  val read : unit -> t
+  val write : t -> unit
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
+  val of_rst : Rst_string.t -> t 
 end = struct
   type t = 
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
-      mo_label               : Non_empty_string.t;
+      mo_label               : MO_label.t;
       n_det                  : Det_number.t;
       n_states               : States_number.t;
       n_states_diag          : States_number.t;
-      n_det_max_jacobi       : Det_number.t;
+      n_det_max_jacobi       : Strictly_positive_int.t;
       threshold_generators   : Threshold.t;
       threshold_selectors    : Threshold.t; 
       read_wf                : bool;
@@ -38,7 +40,7 @@ end = struct
       s2_eig                 : bool;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "determinants";;
@@ -54,6 +56,12 @@ end = struct
     |> N_int_number.of_int
   ;;
 
+  let write_n_int n =
+    N_int_number.to_int n
+    |> Ezfio.set_determinants_n_int
+  ;;
+
+
   let read_bit_kind () =
     if not (Ezfio.has_determinants_bit_kind ()) then
       Lazy.force Qpackage.bit_kind
@@ -64,15 +72,27 @@ end = struct
     |> Bit_kind.of_int
   ;;
 
+  let write_bit_kind b =
+    Bit_kind.to_int b
+    |> Ezfio.set_determinants_bit_kind
+  ;;
+
+
   let read_mo_label () =
     if (not (Ezfio.has_determinants_mo_label ())) then
       Ezfio.get_mo_basis_mo_label ()
       |> Ezfio.set_determinants_mo_label 
       ;
     Ezfio.get_determinants_mo_label ()
-    |> Non_empty_string.of_string
+    |> MO_label.of_string
   ;;
 
+  let write_mo_label l =
+    MO_label.to_string l
+    |> Ezfio.set_determinants_mo_label
+  ;;
+
+
   let read_n_det () =
     if not (Ezfio.has_determinants_n_det ()) then
       Ezfio.set_determinants_n_det 1
@@ -81,6 +101,12 @@ end = struct
     |> Det_number.of_int
   ;;
 
+  let write_n_det n =
+    Det_number.to_int n
+    |> Ezfio.set_determinants_n_det
+  ;;
+
+
   let read_n_states () =
     if not (Ezfio.has_determinants_n_states ()) then
       Ezfio.set_determinants_n_states 1
@@ -89,6 +115,12 @@ end = struct
     |> States_number.of_int
   ;;
 
+  let write_n_states n =
+    States_number.to_int n
+    |> Ezfio.set_determinants_n_states
+  ;;
+
+
   let read_n_states_diag () =
     if not (Ezfio.has_determinants_n_states_diag ()) then
       read_n_states ()
@@ -99,6 +131,12 @@ end = struct
     |> States_number.of_int
   ;;
 
+  let write_n_states_diag n =
+    States_number.to_int n
+    |> Ezfio.set_determinants_n_states_diag
+  ;;
+
+
   let read_n_det_max_jacobi () =
     if not (Ezfio.has_determinants_n_det_max_jacobi ()) then
       get_default "n_det_max_jacobi"
@@ -106,9 +144,15 @@ end = struct
       |> Ezfio.set_determinants_n_det_max_jacobi 
     ;
     Ezfio.get_determinants_n_det_max_jacobi ()
-    |> Det_number.of_int
+    |> Strictly_positive_int.of_int
   ;;
 
+  let write_n_det_max_jacobi n =
+    Strictly_positive_int.to_int n
+    |> Ezfio.set_determinants_n_det_max_jacobi
+  ;;
+
+
   let read_threshold_generators () =
     if not (Ezfio.has_determinants_threshold_generators ()) then
       get_default "threshold_generators"
@@ -119,6 +163,12 @@ end = struct
     |> Threshold.of_float
   ;;
 
+  let write_threshold_generators t =
+    Threshold.to_float t
+    |> Ezfio.set_determinants_threshold_generators
+  ;;
+
+
   let read_threshold_selectors () =
     if not (Ezfio.has_determinants_threshold_selectors ()) then
       get_default "threshold_selectors"
@@ -129,6 +179,12 @@ end = struct
     |> Threshold.of_float
   ;;
 
+  let write_threshold_selectors t =
+    Threshold.to_float t
+    |> Ezfio.set_determinants_threshold_selectors
+  ;;
+
+
   let read_read_wf () =
     if not (Ezfio.has_determinants_read_wf ()) then
       get_default "read_wf"
@@ -138,6 +194,9 @@ end = struct
     Ezfio.get_determinants_read_wf ()
   ;;
 
+  let write_read_wf = Ezfio.set_determinants_read_wf ;;
+    
+
   let read_expected_s2 () =
     if not (Ezfio.has_determinants_expected_s2 ()) then
       begin
@@ -153,6 +212,12 @@ end = struct
     |> Positive_float.of_float
   ;;
 
+  let write_expected_s2 s2 =
+    Positive_float.to_float s2
+    |> Ezfio.set_determinants_expected_s2 
+  ;;
+
+
   let read_s2_eig () =
     if not (Ezfio.has_determinants_s2_eig ()) then
       get_default "s2_eig"
@@ -162,59 +227,80 @@ end = struct
     Ezfio.get_determinants_s2_eig ()
   ;;
 
+  let write_s2_eig = Ezfio.set_determinants_s2_eig ;;
+
+
   let read_psi_coef () =
     if not (Ezfio.has_determinants_psi_coef ()) then
         Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| 1 |] ~data:[1.]
         |> Ezfio.set_determinants_psi_coef 
       ;
-    (Ezfio.get_determinants_psi_coef ()).Ezfio.data
+    Ezfio.get_determinants_psi_coef ()
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.map ~f:Det_coef.of_float
   ;;
 
+  let write_psi_coef ~n_det c =
+    let n_det = Det_number.to_int n_det
+    and c = Array.to_list c
+            |> List.map ~f:Det_coef.to_float
+    in
+    Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| n_det |] ~data:c
+    |> Ezfio.set_determinants_psi_coef 
+  ;;
+
+
   let read_psi_det () =
-    let n_int = read_n_int () in
+    let n_int = read_n_int () 
+    and n_alpha = Ezfio.get_electrons_elec_alpha_num ()
+        |> Elec_alpha_number.of_int 
+    and n_beta = Ezfio.get_electrons_elec_beta_num ()
+        |> Elec_beta_number.of_int 
+    in
     if not (Ezfio.has_determinants_psi_det ()) then
       begin
+        let mo_tot_num = MO_number.get_max () in
         let rec build_data accu =  function
           | 0 -> accu
-          | n -> build_data ((MO_number.of_int n)::accu) (n-1)
+          | n -> build_data ((MO_number.of_int ~max:mo_tot_num n)::accu) (n-1)
         in
-        let det_a = build_data [] (Ezfio.get_electrons_elec_alpha_num ())
+        let det_a = build_data [] (Elec_alpha_number.to_int n_alpha)
           |> Bitlist.of_mo_number_list n_int
-        and det_b = build_data [] (Ezfio.get_electrons_elec_beta_num  ())
+        and det_b = build_data [] (Elec_beta_number.to_int n_beta)
           |> Bitlist.of_mo_number_list n_int
         in
         let data = ( (Bitlist.to_int64_list det_a) @ 
           (Bitlist.to_int64_list det_b) ) 
         in
         Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; 1 |] ~data:data
-        |> Ezfio.set_determinants_psi_det 
+          |> Ezfio.set_determinants_psi_det ;
       end  ;
     let n_int = N_int_number.to_int n_int in
-    let rec transform accu1 accu2 n_rest = function 
+    let psi_det_array = Ezfio.get_determinants_psi_det () in
-      | [] ->
+    let dim = psi_det_array.Ezfio.dim
-          let accu1 = List.rev accu1
+    and data =  Ezfio.flattened_ezfio psi_det_array
-          |> Array.of_list 
-          |> Determinant.of_int64_array
-          in
-          List.rev (accu1::accu2) |> Array.of_list
-      | i::rest ->
-          if (n_rest > 0) then
-            transform (i::accu1) accu2 (n_rest-1) rest
-          else
-            let accu1 = List.rev accu1
-            |> Array.of_list 
-            |> Determinant.of_int64_array
-            in
-            transform [] (accu1::accu2) (2*n_int) rest
     in
-    (Ezfio.get_determinants_psi_det ()).Ezfio.data
+    assert (n_int = dim.(0));
-    |> Ezfio.flattened_ezfio_data
+    assert (dim.(1) = 2);
-    |> Array.to_list
+    assert (dim.(2) = (Det_number.to_int (read_n_det ())));
-    |> transform [] [] (2*n_int)
+    List.init dim.(2) ~f:(fun i ->
+      Array.sub ~pos:(2*n_int*i) ~len:(2*n_int) data)
+    |> List.map ~f:(Determinant.of_int64_array
+      ~n_int:(N_int_number.of_int n_int)
+      ~alpha:n_alpha ~beta:n_beta )
+    |> Array.of_list
   ;;
 
+  let write_psi_det ~n_int ~n_det d =
+    let data = Array.to_list d
+      |> Array.concat 
+      |> Array.to_list
+    in
+    Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data
+    |> Ezfio.set_determinants_psi_det
+  ;;
+
+
   let read () =
     { n_int                  = read_n_int ()                ;
       bit_kind               = read_bit_kind ()             ;
@@ -233,7 +319,108 @@ end = struct
     }
   ;;
 
+  let write { n_int                ;
+              bit_kind             ;
+              mo_label             ;
+              n_det                ;
+              n_states             ;
+              n_states_diag        ;
+              n_det_max_jacobi     ;
+              threshold_generators ;
+              threshold_selectors  ;
+              read_wf              ;
+              expected_s2          ;
+              s2_eig               ;
+              psi_coef             ;
+              psi_det              ;
+            } =
+     write_n_int n_int ;
+     write_bit_kind bit_kind;
+     write_mo_label mo_label;
+     write_n_det n_det;
+     write_n_states n_states;
+     write_n_states_diag n_states_diag;
+     write_n_det_max_jacobi n_det_max_jacobi;
+     write_threshold_generators threshold_generators;
+     write_threshold_selectors threshold_selectors;
+     write_read_wf read_wf;
+     write_expected_s2 expected_s2;
+     write_s2_eig s2_eig;
+     write_psi_coef ~n_det:n_det psi_coef;
+     write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
+  ;;
+
+
+  let to_rst b =
+    let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in
+    let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in
+    let det_text = 
+      List.map2_exn ~f:(fun coef det ->
+        Printf.sprintf "  %F\n%s\n"
+        (Det_coef.to_float coef)
+        (Determinant.to_string ~mo_tot_num:mo_tot_num det 
+         |> String.split ~on:'\n'
+         |> List.map ~f:(fun x -> "  "^x)
+         |> String.concat ~sep:"\n"
+        )
+      ) (Array.to_list b.psi_coef) (Array.to_list b.psi_det)
+      |> String.concat ~sep:"\n"
+    in
+    Printf.sprintf "
+Read the current wave function ::
+
+  read_wf = %s
+
+Label of the MOs on which the determinants were computed ::
+
+  mo_label = %s
+
+Force the selected wave function to be an eigenfunction of S^2.
+If true, input the expected value of S^2 ::
+
+  s2_eig      = %s
+  expected_s2 = %s
+
+Thresholds on generators and selectors (fraction of the norm) ::
+
+  threshold_generators = %s
+  threshold_selectors  = %s
+
+Number of requested states, and number of states used for the
+Davidson diagonalization ::
+
+  n_states      = %s
+  n_states_diag = %s
+
+Maximum size of the Hamiltonian matrix that will be fully diagonalized ::
+
+  n_det_max_jacobi = %s
+
+Number of determinants ::
+
+  n_det = %s
+
+Determinants ::
+
+%s
+"
+     (b.read_wf       |> Bool.to_string)
+     (b.mo_label      |> MO_label.to_string)
+     (b.s2_eig        |> Bool.to_string)
+     (b.expected_s2   |> Positive_float.to_string)
+     (b.threshold_generators |> Threshold.to_string)
+     (b.threshold_selectors |> Threshold.to_string)
+     (b.n_states      |> States_number.to_string)
+     (b.n_states_diag |> States_number.to_string)
+     (b.n_det_max_jacobi |> Strictly_positive_int.to_string)
+     (b.n_det         |> Det_number.to_string)
+     det_text
+     |> Rst_string.of_string
+  ;;
+
   let to_string b =
+    let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in
+    let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in
     Printf.sprintf "
 n_int                  = %s
 bit_kind               = %s
@@ -252,11 +439,11 @@ psi_det                = %s
 "
      (b.n_int         |> N_int_number.to_string)
      (b.bit_kind      |> Bit_kind.to_string)
-     (b.mo_label      |> Non_empty_string.to_string)
+     (b.mo_label      |> MO_label.to_string)
      (b.n_det         |> Det_number.to_string)
      (b.n_states      |> States_number.to_string)
      (b.n_states_diag |> States_number.to_string)
-     (b.n_det_max_jacobi |> Det_number.to_string)
+     (b.n_det_max_jacobi |> Strictly_positive_int.to_string)
      (b.threshold_generators |> Threshold.to_string)
      (b.threshold_selectors |> Threshold.to_string)
      (b.read_wf       |> Bool.to_string)
@@ -264,13 +451,96 @@ psi_det                = %s
      (b.s2_eig        |> Bool.to_string)
      (b.psi_coef  |> Array.to_list |> List.map ~f:Det_coef.to_string
       |> String.concat ~sep:", ")
-     (b.psi_det   |> Array.map ~f:(fun x -> Determinant.to_int64_array x
+     (b.psi_det   |> Array.to_list |> List.map ~f:(Determinant.to_string
-      |> Array.map ~f:(fun x-> 
+       ~mo_tot_num:mo_tot_num) |> String.concat ~sep:"\n\n")
-          Int64.to_string x )|> Array.to_list |>
+  ;;
-       String.concat ~sep:", ") |> Array.to_list
+
-      |> String.concat ~sep:" | ")
+  let of_rst r =
-     ;
+    let r = Rst_string.to_string r
-;;
+    in
+
+    (* Split into header and determinants data *)
+    let idx = String.substr_index_exn r ~pos:0 ~pattern:"\nDeterminants"
+    in
+    let (header, dets) = 
+       (String.prefix r idx, String.suffix r ((String.length r)-idx) )
+    in
+
+    (* Handle header *)
+    let header = r
+    |> String.split ~on:'\n'
+    |> List.filter ~f:(fun line ->
+        if (line = "") then
+          false
+        else
+          ( (String.contains line '=') && (line.[0] = ' ') )
+       )
+    |> List.map ~f:(fun line ->
+        "("^(
+        String.tr line ~target:'=' ~replacement:' '
+        |> String.strip
+        )^")" )
+    |> String.concat 
+    in
+
+    (* Handle determinant coefs *)
+    let dets = match ( dets
+    |> String.split ~on:'\n'
+    |> List.map ~f:(String.strip)
+    ) with 
+    | _::lines -> lines 
+    | _ -> failwith "Error in determinants"
+    in
+
+    let psi_coef = 
+      let rec read_coefs accu = function
+      | [] -> List.rev accu
+      | ""::c::tail -> 
+          read_coefs (c::accu) tail
+      | _::tail -> read_coefs accu tail
+      in
+      let a = read_coefs [] dets
+      |> String.concat ~sep:" "
+      in
+      "(psi_coef ("^a^"))"
+    in
+
+    (* Handle determinants *)
+    let psi_det = 
+      let n_alpha = Ezfio.get_electrons_elec_alpha_num ()
+        |> Elec_alpha_number.of_int 
+      and n_beta = Ezfio.get_electrons_elec_beta_num ()
+        |> Elec_beta_number.of_int 
+      in
+      let rec read_dets accu = function
+      | [] -> List.rev accu
+      | ""::c::alpha::beta::tail -> 
+          begin
+            let alpha = String.rev alpha |> Bitlist.of_string ~zero:'-' ~one:'+'
+            and beta  = String.rev beta  |> Bitlist.of_string ~zero:'-' ~one:'+'
+            in
+            let newdet = Determinant.of_bitlist_couple 
+              ~alpha:n_alpha ~beta:n_beta (alpha,beta) 
+            |> Determinant.sexp_of_t |> Sexplib.Sexp.to_string
+            in
+            read_dets (newdet::accu) tail
+          end
+      | _::tail -> read_dets accu tail
+      in
+      let a = read_dets [] dets
+      |> String.concat 
+      in
+      "(psi_det ("^a^"))"
+    in
+
+    let bitkind = Printf.sprintf "(bit_kind %d)" (Lazy.force Qpackage.bit_kind
+      |> Bit_kind.to_int)
+    and n_int = Printf.sprintf "(n_int %d)" (N_int_number.get_max ()) in
+    let s = String.concat [ header ; bitkind ; n_int ; psi_coef ; psi_det]
+    in
+    Sexp.of_string ("("^s^")")
+    |> t_of_sexp
+  ;;
 
 end
 

ocaml/Input_electrons.ml

@@ -6,17 +6,19 @@ module Electrons : sig
   type t = 
     { elec_alpha_num     : Elec_alpha_number.t;
       elec_beta_num      : Elec_beta_number.t;
-      elec_num           : Elec_number.t;
+    } with sexp
-    }
   ;;
-  val read : unit -> t
+  val read  : unit -> t
+  val write : t -> unit
+  val read_elec_num : unit -> Elec_number.t
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
+  val of_rst : Rst_string.t -> t
 end = struct
   type t = 
     { elec_alpha_num     : Elec_alpha_number.t;
       elec_beta_num      : Elec_beta_number.t;
-      elec_num           : Elec_number.t;
+    } with sexp
-    }
   ;;
 
   let get_default = Qpackage.get_ezfio_default "electrons";;
@@ -26,11 +28,22 @@ end = struct
     |> Elec_alpha_number.of_int
   ;;
 
+  let write_elec_alpha_num n = 
+    Elec_alpha_number.to_int n
+    |> Ezfio.set_electrons_elec_alpha_num 
+  ;;
+
+
   let read_elec_beta_num() = 
     Ezfio.get_electrons_elec_beta_num ()
     |> Elec_beta_number.of_int
   ;;
 
+  let write_elec_beta_num n = 
+    Elec_beta_number.to_int n
+    |> Ezfio.set_electrons_elec_beta_num 
+  ;;
+
   let read_elec_num () = 
     let na = Ezfio.get_electrons_elec_alpha_num ()
     and nb = Ezfio.get_electrons_elec_beta_num  ()
@@ -42,19 +55,57 @@ end = struct
   let read () = 
     { elec_alpha_num      = read_elec_alpha_num ();
       elec_beta_num       = read_elec_beta_num ();
-      elec_num            = read_elec_num ();
     }
   ;;
 
-  let to_string b =
+  let write { elec_alpha_num ; elec_beta_num } =
+    write_elec_alpha_num elec_alpha_num;
+    write_elec_beta_num  elec_beta_num;
+  ;;
+
+
+  let to_rst b =
     Printf.sprintf "
-elec_alpha_num     = %s
+Spin multiplicity is %s.
+
+Number of alpha and beta electrons ::
+
+  elec_alpha_num = %s
+  elec_beta_num  = %s
+
+"
+        (Multiplicity.of_alpha_beta b.elec_alpha_num b.elec_beta_num
+         |> Multiplicity.to_string)
+        (Elec_alpha_number.to_string b.elec_alpha_num)
+        (Elec_beta_number.to_string b.elec_beta_num)
+    |> Rst_string.of_string
+  ;;
+
+  let to_string b =
+    Printf.sprintf "elec_alpha_num     = %s
 elec_beta_num      = %s
 elec_num           = %s
 "
         (Elec_alpha_number.to_string b.elec_alpha_num)
         (Elec_beta_number.to_string b.elec_beta_num)
-        (Elec_number.to_string b.elec_num)
+        (Elec_number.to_string (read_elec_num ()))
+  ;;
+
+  let of_rst s =
+    let s = Rst_string.to_string s
+    |> String.split ~on:'\n'
+    |> List.filter ~f:(fun line ->
+        String.contains line '=')
+    |> List.map ~f:(fun line ->
+        "("^(
+        String.tr line ~target:'=' ~replacement:' '
+        )^")" )
+    |> String.concat
+    in
+    Sexp.of_string ("("^s^")")
+    |> t_of_sexp
+  ;;
+
 end
 
 

ocaml/Input_full_ci.ml

@@ -4,19 +4,22 @@ open Core.Std;;
 
 module Full_ci : sig
   type t = 
-    { n_det_max_fci      : Det_number.t;
+    { n_det_max_fci      : Det_number_max.t;
       pt2_max            : PT2_energy.t;
       do_pt2_end         : bool;
-    }
+    } with sexp
   ;;
-  val read : unit -> t
+  val read  : unit -> t
+  val write : t-> unit
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
+  val of_rst : Rst_string.t -> t
 end = struct
   type t = 
-    { n_det_max_fci      : Det_number.t;
+    { n_det_max_fci      : Det_number_max.t;
       pt2_max            : PT2_energy.t;
       do_pt2_end         : bool;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "full_ci";;
@@ -28,7 +31,12 @@ end = struct
        |> Ezfio.set_full_ci_n_det_max_fci
     ;
     Ezfio.get_full_ci_n_det_max_fci ()
-    |> Det_number.of_int
+    |> Det_number_max.of_int
+  ;;
+
+  let write_n_det_max_fci ndet = 
+    Det_number_max.to_int ndet
+    |> Ezfio.set_full_ci_n_det_max_fci
   ;;
 
   let read_pt2_max () = 
@@ -41,6 +49,11 @@ end = struct
     |> PT2_energy.of_float
   ;;
 
+  let write_pt2_max pt2_max = 
+    PT2_energy.to_float pt2_max
+    |> Ezfio.set_full_ci_pt2_max
+  ;;
+
   let read_do_pt2_end () = 
     if not (Ezfio.has_full_ci_do_pt2_end ()) then
        get_default "do_pt2_end"
@@ -50,6 +63,10 @@ end = struct
     Ezfio.get_full_ci_do_pt2_end ()
   ;;
 
+  let write_do_pt2_end =
+    Ezfio.set_full_ci_do_pt2_end
+  ;;
+
 
   let read () = 
     { n_det_max_fci    = read_n_det_max_fci ();
@@ -58,15 +75,64 @@ end = struct
     }
   ;;
 
+
+  let write { n_det_max_fci ;
+              pt2_max       ;
+              do_pt2_end    ;
+            } = 
+     write_n_det_max_fci n_det_max_fci;
+     write_pt2_max       pt2_max;
+     write_do_pt2_end    do_pt2_end;
+  ;;
+
   let to_string b =
     Printf.sprintf "
-n_det_max_fci      = %s
+n_det_max_fci = %s
-pt2_max            = %s
+pt2_max = %s
-do_pt2_end         = %s
+do_pt2_end = %s
 "
-        (Det_number.to_string b.n_det_max_fci)
+        (Det_number_max.to_string b.n_det_max_fci)
         (PT2_energy.to_string b.pt2_max)
         (Bool.to_string b.do_pt2_end)
+    ;;
+
+  let to_rst b =
+    Printf.sprintf "
+Stop when the `n_det` > `n_det_max_fci` ::
+
+  n_det_max_fci = %s
+
+Stop when -E(PT2) < `pt2_max` ::
+
+  pt2_max = %s
+
+Compute E(PT2) at the end ::
+
+  do_pt2_end = %s
+
+"
+        (Det_number_max.to_string b.n_det_max_fci)
+        (PT2_energy.to_string b.pt2_max)
+        (Bool.to_string b.do_pt2_end)
+    |> Rst_string.of_string
+    ;;
+
+  let of_rst s =
+    let s = Rst_string.to_string s
+    |> String.split ~on:'\n'
+    |> List.filter ~f:(fun line ->
+        String.contains line '=')
+    |> List.map ~f:(fun line ->
+        "("^(
+        String.tr line ~target:'=' ~replacement:' '
+        )^")" )
+    |> String.concat
+    in
+    Sexp.of_string ("("^s^")")
+    |> t_of_sexp
+  ;;
+
+
 end
 
 

ocaml/Input_hartree_fock.ml

@@ -6,15 +6,18 @@ module Hartree_fock : sig
   type t = 
     { n_it_scf_max       : Strictly_positive_int.t;
       thresh_scf         : Threshold.t;
-    }
+    } with sexp
   ;;
-  val read : unit -> t
+  val read  : unit -> t
+  val write : t -> unit
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
+  val of_rst : Rst_string.t -> t 
 end = struct
   type t = 
     { n_it_scf_max       : Strictly_positive_int.t;
       thresh_scf         : Threshold.t;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "hartree_fock";;
@@ -29,14 +32,25 @@ end = struct
     |> Strictly_positive_int.of_int
   ;;
 
-  let read_thresh_scf() = 
+  let write_n_it_scf_max n_it_scf_max = 
+    Strictly_positive_int.to_int n_it_scf_max
+    |> Ezfio.set_hartree_fock_n_it_scf_max
+  ;;
+
+  let read_thresh_scf () = 
     if not (Ezfio.has_hartree_fock_thresh_scf()) then
        get_default "thresh_scf"
        |> Float.of_string
        |> Ezfio.set_hartree_fock_thresh_scf
     ;
     Ezfio.get_hartree_fock_thresh_scf ()
-    |> Threshold.of_float ;;
+    |> Threshold.of_float
+  ;;
+
+  let write_thresh_scf thresh_scf =
+    Threshold.to_float thresh_scf
+    |> Ezfio.set_hartree_fock_thresh_scf
+  ;;
 
 
   let read () = 
@@ -45,13 +59,56 @@ end = struct
     }
   ;;
 
+
+  let write { n_it_scf_max ;
+              thresh_scf ;
+            } =
+    write_n_it_scf_max n_it_scf_max;
+    write_thresh_scf   thresh_scf
+  ;;
+
+
   let to_string b =
     Printf.sprintf "
-n_it_scf_max       = %s
+n_it_scf_max = %s
-thresh_scf         = %s
+thresh_scf = %s
+"
+      (Strictly_positive_int.to_string b.n_it_scf_max)
+      (Threshold.to_string b.thresh_scf)
+  ;;
+
+  let to_rst b =
+    Printf.sprintf "
+Max number of SCF iterations ::
+
+  n_it_scf_max = %s
+
+SCF convergence criterion (on energy) ::
+
+  thresh_scf = %s
+
 "
         (Strictly_positive_int.to_string b.n_it_scf_max)
         (Threshold.to_string b.thresh_scf)
+    |> Rst_string.of_string
+  ;;
+
+  let of_rst s =
+    let s = Rst_string.to_string s
+    |> String.split ~on:'\n'
+    |> List.filter ~f:(fun line ->
+        String.contains line '=')
+    |> List.map ~f:(fun line ->
+        "("^(
+        String.tr line ~target:'=' ~replacement:' '
+        )^")" )
+    |> String.concat
+    in
+    Sexp.of_string ("("^s^")")
+    |> t_of_sexp
+  ;;
+
+
 end
 
 

ocaml/Input_mo_basis.ml

@@ -5,30 +5,31 @@ open Core.Std;;
 module Mo_basis : sig
   type t = 
     { mo_tot_num      : MO_number.t ;
-      mo_label        : Non_empty_string.t;
+      mo_label        : MO_label.t;
-      mo_occ          : Positive_float.t array;
+      mo_occ          : MO_occ.t array;
-      mo_coef         : MO_coef.t array;
+      mo_coef         : (MO_coef.t array) array;
-    }
+    } with sexp
   ;;
   val read : unit -> t
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
 end = struct
   type t = 
     { mo_tot_num      : MO_number.t ;
-      mo_label        : Non_empty_string.t;
+      mo_label        : MO_label.t;
-      mo_occ          : Positive_float.t array;
+      mo_occ          : MO_occ.t array;
-      mo_coef         : MO_coef.t array;
+      mo_coef         : (MO_coef.t array) array;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "mo_basis";;
 
   let read_mo_label () = 
     if not (Ezfio.has_mo_basis_mo_label ()) then
-       Ezfio.set_mo_basis_mo_label "Unknown"
+       Ezfio.set_mo_basis_mo_label "None"
     ;
     Ezfio.get_mo_basis_mo_label ()
-    |> Non_empty_string.of_string
+    |> MO_label.of_string
   ;;
 
   let read_mo_tot_num () =
@@ -50,15 +51,20 @@ end = struct
           ~dim:[| mo_tot_num |] ~data:data
         |> Ezfio.set_mo_basis_mo_occ
       end;
-    (Ezfio.get_mo_basis_mo_occ () ).Ezfio.data
+    Ezfio.flattened_ezfio (Ezfio.get_mo_basis_mo_occ () )
-    |> Ezfio.flattened_ezfio_data
+    |> Array.map ~f:MO_occ.of_float
-    |> Array.map ~f:Positive_float.of_float
   ;;
 
   let read_mo_coef () =
-    (Ezfio.get_mo_basis_mo_coef () ).Ezfio.data
+    let a = Ezfio.get_mo_basis_mo_coef () 
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.map ~f:MO_coef.of_float
+    in
+    let mo_tot_num = read_mo_tot_num () |> MO_number.to_int in
+    let ao_num = (Array.length a)/mo_tot_num in
+    Array.init mo_tot_num ~f:(fun j ->
+      Array.sub ~pos:(j*ao_num) ~len:(ao_num) a
+    )
   ;;
 
   let read () =
@@ -69,6 +75,91 @@ end = struct
     }
   ;;
 
+  let mo_coef_to_string mo_coef =
+    let ao_num = Array.length mo_coef.(0) 
+    and mo_tot_num = Array.length mo_coef in
+    let rec print_five imin imax =
+      match (imax-imin+1) with 
+      | 1 ->
+          let header = [ Printf.sprintf "  #%15d" (imin+1) ; ] in
+          let new_lines = 
+            List.init ao_num ~f:(fun i ->
+              Printf.sprintf "  %3d %15.10f " (i+1)
+              (MO_coef.to_float mo_coef.(imin  ).(i)) )
+          in header @ new_lines
+      | 2 ->
+          let header = [ Printf.sprintf "  #%15d %15d" (imin+1) (imin+2) ; ] in
+          let new_lines = 
+            List.init ao_num ~f:(fun i ->
+              Printf.sprintf "  %3d %15.10f %15.10f" (i+1)
+              (MO_coef.to_float mo_coef.(imin  ).(i))
+              (MO_coef.to_float mo_coef.(imin+1).(i)) )
+          in header @ new_lines
+      | 3 ->
+          let header = [ Printf.sprintf "  #%15d %15d %15d"
+              (imin+1) (imin+2) (imin+3); ] in
+          let new_lines = 
+            List.init ao_num ~f:(fun i ->
+              Printf.sprintf "  %3d %15.10f %15.10f %15.10f" (i+1)
+              (MO_coef.to_float mo_coef.(imin  ).(i))
+              (MO_coef.to_float mo_coef.(imin+1).(i))
+              (MO_coef.to_float mo_coef.(imin+2).(i)) )
+          in header @ new_lines
+      | 4 ->
+          let header = [ Printf.sprintf "  #%15d %15d %15d %15d"
+              (imin+1) (imin+2) (imin+3) (imin+4) ; ] in
+          let new_lines = 
+            List.init ao_num ~f:(fun i ->
+              Printf.sprintf "  %3d %15.10f %15.10f %15.10f %15.10f" (i+1)
+              (MO_coef.to_float mo_coef.(imin  ).(i))
+              (MO_coef.to_float mo_coef.(imin+1).(i))
+              (MO_coef.to_float mo_coef.(imin+2).(i))
+              (MO_coef.to_float mo_coef.(imin+3).(i)) )
+          in header @ new_lines
+      | 5 ->
+          let header = [ Printf.sprintf "  #%15d %15d %15d %15d %15d"
+              (imin+1) (imin+2) (imin+3) (imin+4) (imin+5) ; ] in
+          let new_lines = 
+            List.init ao_num ~f:(fun i ->
+              Printf.sprintf "  %3d %15.10f %15.10f %15.10f %15.10f %15.10f" (i+1)
+              (MO_coef.to_float mo_coef.(imin  ).(i))
+              (MO_coef.to_float mo_coef.(imin+1).(i))
+              (MO_coef.to_float mo_coef.(imin+2).(i))
+              (MO_coef.to_float mo_coef.(imin+3).(i))
+              (MO_coef.to_float mo_coef.(imin+4).(i)) )
+          in header @ new_lines 
+      | _ -> assert false
+    in
+    let rec create_list accu i = 
+      if (i+4 < mo_tot_num) then
+        create_list ( (print_five i (i+3) |> String.concat ~sep:"\n")::accu ) (i+4)
+      else
+        (print_five i (mo_tot_num-1) |> String.concat ~sep:"\n")::accu |> List.rev
+    in
+    create_list [] 0 |> String.concat ~sep:"\n\n"
+  ;;
+
+  let to_rst b =
+    Printf.sprintf "
+Label of the molecular orbitals ::
+
+  mo_label = %s
+
+Total number of MOs ::
+
+  mo_tot_num = %s
+
+MO coefficients ::
+
+%s
+"
+    (MO_label.to_string b.mo_label)
+    (MO_number.to_string b.mo_tot_num)
+    (mo_coef_to_string b.mo_coef)
+    |> Rst_string.of_string
+
+  ;;
+
   let to_string b =
     Printf.sprintf "
 mo_label        = %s
@@ -76,12 +167,15 @@ mo_tot_num      = \"%s\"
 mo_occ          = %s
 mo_coef         = %s
 "
-    (Non_empty_string.to_string b.mo_label)
+    (MO_label.to_string b.mo_label)
     (MO_number.to_string b.mo_tot_num)
     (b.mo_occ |> Array.to_list |> List.map
-      ~f:(Positive_float.to_string) |> String.concat ~sep:", " )
+      ~f:(MO_occ.to_string) |> String.concat ~sep:", " )
-    (b.mo_coef |> Array.to_list |> List.map
+    (b.mo_coef |> Array.map
-      ~f:(MO_coef.to_string) |> String.concat ~sep:", " )
+      ~f:(fun x-> Array.map ~f:MO_coef.to_string x |> String.concat_array
+      ~sep:"," ) |>
+      String.concat_array ~sep:"\n" )
+  ;;
 
 end
 

ocaml/Input_nuclei.ml

@@ -8,43 +8,80 @@ module Nuclei : sig
       nucl_label      : Element.t array;
       nucl_charge     : Charge.t array;
       nucl_coord      : Point3d.t array;
-    }
+    } with sexp
   ;;
-  val read : unit -> t
+  val read  : unit -> t
+  val write : t -> unit 
   val to_string : t -> string
+  val to_rst : t -> Rst_string.t
+  val of_rst : Rst_string.t -> t
 end = struct
   type t = 
     { nucl_num        : Nucl_number.t ;
       nucl_label      : Element.t array;
       nucl_charge     : Charge.t array;
       nucl_coord      : Point3d.t array;
-    }
+    } with sexp
   ;;
 
   let get_default = Qpackage.get_ezfio_default "nuclei";;
 
   let read_nucl_num () = 
-    Ezfio.get_nuclei_nucl_num ()
+    let nmax = Nucl_number.get_max () in
-    |> Nucl_number.of_int
+    Nucl_number.of_int ~max:nmax nmax
   ;;
 
+  let write_nucl_num n = 
+    Nucl_number.to_int n
+    |> Ezfio.set_nuclei_nucl_num
+  ;;
+    
+
   let read_nucl_label () =
-    (Ezfio.get_nuclei_nucl_label ()).Ezfio.data
+    Ezfio.get_nuclei_nucl_label ()
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.map ~f:Element.of_string 
   ;;
 
+  let write_nucl_label ~nucl_num labels =
+    let nucl_num = 
+       Nucl_number.to_int nucl_num
+    in
+    let labels =
+       Array.to_list labels 
+       |> List.map ~f:Element.to_string
+    in
+    Ezfio.ezfio_array_of_list ~rank:1 
+       ~dim:[| nucl_num |] ~data:labels
+    |> Ezfio.set_nuclei_nucl_label
+  ;;
+
+
   let read_nucl_charge () =
-    (Ezfio.get_nuclei_nucl_charge () ).Ezfio.data
+    Ezfio.get_nuclei_nucl_charge ()
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     |> Array.map ~f:Charge.of_float
   ;;
 
+  let write_nucl_charge ~nucl_num charges =
+    let nucl_num = 
+       Nucl_number.to_int nucl_num
+    in
+    let charges =
+       Array.to_list charges
+       |> List.map ~f:Charge.to_float
+    in
+    Ezfio.ezfio_array_of_list ~rank:1 
+       ~dim:[| nucl_num |] ~data:charges
+    |> Ezfio.set_nuclei_nucl_charge
+  ;;
+
+
   let read_nucl_coord () =
     let nucl_num = Nucl_number.to_int (read_nucl_num ()) in
     let raw_data = 
-    (Ezfio.get_nuclei_nucl_coord() ).Ezfio.data
+    Ezfio.get_nuclei_nucl_coord()
-    |> Ezfio.flattened_ezfio_data
+    |> Ezfio.flattened_ezfio
     in
     let zero = Point3d.of_string Units.Bohr "0. 0. 0." in
     let result = Array.create nucl_num zero in
@@ -57,6 +94,22 @@ end = struct
     result
   ;;
 
+  let write_nucl_coord ~nucl_num coord =
+    let nucl_num = 
+       Nucl_number.to_int nucl_num
+    in
+    let coord = Array.to_list coord in
+    let coord = 
+       (List.map ~f:(fun x-> x.Point3d.x) coord) @
+       (List.map ~f:(fun x-> x.Point3d.y) coord) @
+       (List.map ~f:(fun x-> x.Point3d.z) coord) 
+    in
+    Ezfio.ezfio_array_of_list ~rank:2 
+       ~dim:[| nucl_num ; 3 |] ~data:coord
+    |> Ezfio.set_nuclei_nucl_coord
+  ;;
+
+
   let read () =
     { nucl_num        = read_nucl_num ();
       nucl_label      = read_nucl_label () ;
@@ -65,6 +118,18 @@ end = struct
     }
   ;;
 
+  let write { nucl_num    ;
+              nucl_label  ;
+              nucl_charge ;
+              nucl_coord  ;
+            } =
+    write_nucl_num nucl_num ;
+    write_nucl_label  ~nucl_num:nucl_num nucl_label;
+    write_nucl_charge ~nucl_num:nucl_num nucl_charge;
+    write_nucl_coord  ~nucl_num:nucl_num nucl_coord;
+  ;;         
+
+
   let to_string b =
     Printf.sprintf "
 nucl_num         = %s
@@ -79,6 +144,75 @@ nucl_coord       = %s
       ~f:(Charge.to_string) |> String.concat ~sep:", " )
     (b.nucl_coord  |> Array.to_list |> List.map
       ~f:(Point3d.to_string Units.Bohr) |> String.concat ~sep:"\n" )
+  ;;
+
+
+   let to_rst b = 
+     let nucl_num = Nucl_number.to_int b.nucl_num in
+     let text = 
+       ( Printf.sprintf "  %d\n  "
+         nucl_num 
+       ) :: (
+       List.init nucl_num ~f:(fun i->
+         Printf.sprintf "  %-3s  %d   %s"
+          (b.nucl_label.(i)  |> Element.to_string)
+          (b.nucl_charge.(i) |> Charge.to_int )
+          (b.nucl_coord.(i)  |> Point3d.to_string Units.Angstrom) )
+      ) |> String.concat ~sep:"\n"
+     in
+     Printf.sprintf "
+Nuclear coordinates in xyz format (Angstroms) ::
+
+%s
+
+" text
+     |> Rst_string.of_string
+  ;;
+     
+  let of_rst s =
+    let l = Rst_string.to_string s 
+    |> String.split ~on:'\n'
+    in
+    (* Find lines containing the xyz data *)
+    let rec extract_begin = function
+    | [] -> raise Not_found
+    | line::tail ->
+      let line = String.strip line in
+      if (String.length line > 3) &&
+        (String.sub line ~pos:((String.length line)-2)
+            ~len:2 = "::") then
+           tail
+      else
+         extract_begin tail
+    in
+    (* Create a list of Atom.t *)
+    let nmax = Nucl_number.get_max () in
+    let atom_list = 
+      match (extract_begin l) with 
+      | _ :: nucl_num :: title :: lines ->
+        begin
+          let nucl_num = nucl_num
+          |> String.strip
+          |> Int.of_string
+          |> Nucl_number.of_int ~max:nmax
+          and lines = Array.of_list lines
+          in
+          List.init (Nucl_number.to_int nucl_num) ~f:(fun i ->
+            Atom.of_string Units.Angstrom lines.(i))
+        end
+      | _ -> failwith "Error in xyz format"
+    in
+    (* Create the Nuclei.t data structure *)
+    { nucl_num = List.length atom_list
+        |> Nucl_number.of_int ~max:nmax;
+      nucl_label = List.map atom_list ~f:(fun x ->
+        x.Atom.element) |> Array.of_list ;
+      nucl_charge = List.map atom_list ~f:(fun x ->
+        x.Atom.charge ) |> Array.of_list ;
+      nucl_coord = List.map atom_list ~f:(fun x ->
+        x.Atom.coord ) |> Array.of_list ;
+    }
+  ;;
 
 end
 

ocaml/Long_basis.ml

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

ocaml/Long_basis.mli

@@ -5,12 +5,16 @@ 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
+type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list with sexp
 
 (** Transform a basis to a long basis *)
 val of_basis :
   (Gto.t * Nucl_number.t) list -> (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list
 
+(** Transform a long basis to a basis *)
+val to_basis :
+  (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list 
+
 (** Convert the basis into its string representation *)
 val to_string :
   (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> string

ocaml/MO_class.ml

@@ -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
 
 
 let to_string  x = 

ocaml/MO_class.mli

@@ -0,0 +1,21 @@
+type t =
+  | Core     of Qptypes.MO_number.t list
+  | Inactive of Qptypes.MO_number.t list
+  | Active   of Qptypes.MO_number.t list
+  | Virtual  of Qptypes.MO_number.t list
+  | Deleted  of Qptypes.MO_number.t list
+with sexp
+
+
+(** Create different excitation classes *)
+val create_core     : string -> t
+val create_inactive : string -> t
+val create_active   : string -> t
+val create_virtual  : string -> t
+val create_deleted  : string -> t
+
+(** Convert to a Bitlist.t *)
+val to_bitlist      : Qptypes.N_int_number.t -> t -> Bitlist.t
+
+(** Convert to string for printing *)
+val to_string : t -> string

ocaml/MO_label.ml

@@ -0,0 +1,32 @@
+open Core.Std;;
+
+type t =
+| Guess
+| Canonical
+| Natural
+| Localized
+| Orthonormalized
+| None
+with sexp
+;;
+
+let to_string = function
+  | Guess     -> "Guess"
+  | Canonical -> "Canonical"
+  | Orthonormalized -> "Orthonormalized"
+  | Natural   -> "Natural"
+  | Localized -> "Localized"
+  | None      -> "None"
+;;
+
+let of_string  s = 
+  match String.lowercase s with 
+  | "guess"     -> Guess
+  | "canonical" -> Canonical
+  | "natural"   -> Natural
+  | "localized" -> Localized
+  | "orthonormalized" -> Orthonormalized
+  | "none"      -> None
+  | _ -> failwith "MO_label should be one of:
+    Guess | Orthonormalized | Canonical | Natural | Localized | None."
+;;

ocaml/MO_label.mli

@@ -0,0 +1,15 @@
+type t =
+  | Guess
+  | Canonical
+  | Natural
+  | Localized
+  | Orthonormalized
+  | None
+with sexp
+
+(** String representation *)
+val to_string : t -> string
+
+(** Build from string representation *)
+val of_string : string -> t
+

ocaml/Makefile

@@ -1,5 +1,3 @@
-#TODO : Opam auto-installer in makefile
-
 # Check if QPACKAGE_ROOT is defined
 
 ifndef QPACKAGE_ROOT
@@ -12,8 +10,8 @@ endif
 
 LIBS=
 PKGS=
-OCAMLCFLAGS=-g 
+OCAMLCFLAGS="-g"
-OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags -g 
+OCAMLBUILD=ocamlbuild -j 0 -syntax camlp4o -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS)
 MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml
 MLIFILES=$(wildcard *.mli) 
 ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml))
@@ -21,11 +19,21 @@ ALL_EXE=$(patsubst %.ml,%.native,$(wildcard qp_*.ml))
 
 .PHONY: executables default
 
+
 default: $(ALL_TESTS) $(ALL_EXE)
 
 executables: 
 	$(MAKE) -C $(QPACKAGE_ROOT)/data executables
 
+external_libs:
+	opam install cryptokit core
+
+qpackage.odocl: $(MLIFILES)
+	ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl
+
+doc: qpackage.odocl
+	$(OCAMLBUILD) qpackage.docdir/index.html -use-ocamlfind  $(PKGS) 
+
 %.inferred.mli: $(MLFILES)
 	$(OCAMLBUILD) $*.inferred.mli -use-ocamlfind  $(PKGS) 
 	mv _build/$*.inferred.mli .

ocaml/Molecule.ml

@@ -7,7 +7,7 @@ type t = {
   nuclei     : Atom.t list ;
   elec_alpha : Elec_alpha_number.t ;
   elec_beta  : Elec_beta_number.t ;
-}
+} with sexp
 
 let get_charge { nuclei  ; elec_alpha ; elec_beta } =
   let result = (Elec_alpha_number.to_int elec_alpha) +
@@ -25,7 +25,8 @@ let get_multiplicity m =
 ;;
 
 let get_nucl_num m =
-  Nucl_number.of_int (List.length m.nuclei)
+  let nmax = (List.length m.nuclei) in
+  Nucl_number.of_int nmax ~max:nmax
 ;;
 
 let name m = 
@@ -78,7 +79,7 @@ let to_string m =
 ;;
 
 let of_xyz_string
-    ?(charge=0) ?(multiplicity=(Multiplicity.of_int 1))
+    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
     ?(units=Units.Angstrom)
     s =
   let l = String.split s ~on:'\n'
@@ -90,7 +91,7 @@ let of_xyz_string
         elec_alpha=(Elec_alpha_number.of_int 1) ;
         elec_beta=(Elec_beta_number.of_int 0) } 
       |> Charge.to_int 
-      ) + 1 - charge 
+      ) + 1 - (Charge.to_int charge)
       |> Elec_number.of_int 
   in
   let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity in
@@ -112,10 +113,16 @@ let of_xyz_string
 
 
 let of_xyz_file
-    ?(charge=0) ?(multiplicity=(Multiplicity.of_int 1))
+    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
+    ?(units=Units.Angstrom)
     filename =
   let (_,buffer) = In_channel.read_all filename 
   |> String.lsplit2_exn ~on:'\n' in
   let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in
-  of_xyz_string ~charge:charge ~multiplicity:multiplicity buffer
+  of_xyz_string ~charge:charge ~multiplicity:multiplicity 
+    ~units:units buffer
+;;
+
+include To_md5;;
+let to_md5 = to_md5 sexp_of_t
 ;;

ocaml/Molecule.mli

@@ -0,0 +1,38 @@
+exception MultiplicityError of string
+
+type t = {
+  nuclei     : Atom.t list;
+  elec_alpha : Qptypes.Elec_alpha_number.t;
+  elec_beta  : Qptypes.Elec_beta_number.t;
+} with sexp
+
+(** Returns the charge of the molecule *)
+val get_charge       : t -> Charge.t
+
+(** Returns the multiplicity of the molecule *)
+val get_multiplicity : t -> Multiplicity.t
+
+(** Returns the number of nuclei *)
+val get_nucl_num : t -> Qptypes.Nucl_number.t
+
+(** The name of the molecule *)
+val name : t -> string
+
+(** Conversion for printing *)
+val to_string : t -> string
+
+
+(** Creates a molecule from an xyz file *)
+val of_xyz_file :
+  ?charge:Charge.t -> 
+  ?multiplicity:Multiplicity.t ->
+  ?units:Units.units -> string -> t
+
+(** Creates a molecule from an xyz file in a string *)
+val of_xyz_string :
+  ?charge:Charge.t ->
+  ?multiplicity:Multiplicity.t ->
+  ?units:Units.units -> string -> t
+
+(** Computes the MD5 hash *)
+val to_md5 : t -> Qptypes.MD5.t

ocaml/Multiplicity.ml

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

ocaml/Multiplicity.mli

@@ -0,0 +1,19 @@
+type t = Qptypes.Strictly_positive_int.t with sexp
+
+(** Conversion from int *)
+val of_int : int -> t
+val to_int : t -> int
+
+(** Computation from the number of alpha and beta electrons *)
+val of_alpha_beta :
+  Qptypes.Elec_alpha_number.t ->
+  Qptypes.Elec_beta_number.t -> t
+
+(** Generation of the number of alpha and beta electrons *)
+val to_alpha_beta :
+  Qptypes.Elec_number.t ->  t ->
+  Qptypes.Elec_alpha_number.t * Qptypes.Elec_beta_number.t
+
+(** Conversion to string for printing *)
+val to_string : t-> string
+

ocaml/Point3d.ml

@@ -1,10 +1,11 @@
 open Core.Std;;
+open Qptypes;;
 
 type t = {
   x : float ;
   y : float ;
   z : float ;
-}
+} with sexp
 
 (** Read x y z coordinates in string s with units u *)
 let of_string u s =
@@ -28,9 +29,10 @@ let distance2 p1 p2 =
   let { x=x1 ; y=y1 ; z=z1 } = p1
   and { x=x2 ; y=y2 ; z=z2 } = p2 in
   (x2-.x1)*.(x2-.x1) +. (y2-.y1)*.(y2-.y1) +. (z2-.z1)*.(z2-.z1) 
+  |> Positive_float.of_float
 ;;
 
-let distance p1 p2 = sqrt (distance2 p1 p2)
+let distance p1 p2 = sqrt (Positive_float.to_float (distance2 p1 p2))
 ;;
 
 let to_string u p =
@@ -39,6 +41,6 @@ let to_string u p =
     | Units.Angstrom -> Units.bohr_to_angstrom
   in
   let { x=x ; y=y ; z=z } = p in
-  Printf.sprintf "%f  %f  %f" (x*.f) (y*.f) (z*.f)
+  Printf.sprintf "%16.8f %16.8f %16.8f" (x*.f) (y*.f) (z*.f)
 ;;
   

ocaml/Point3d.mli

@@ -0,0 +1,17 @@
+type t =
+{ x : float;
+  y : float;
+  z : float;
+} with sexp
+
+(** Create from an xyz string *)
+val of_string : Units.units -> string -> t
+
+(** Convert to a string for printing *)
+val to_string : Units.units -> t -> string
+
+(** Computes the squared distance between 2 points *)
+val distance2 : t -> t -> Qptypes.Positive_float.t
+
+(** Computes the distance between 2 points *)
+val distance : t -> t -> float

ocaml/Primitive.ml

@@ -4,7 +4,7 @@ open Core.Std;;
 type t =
 { sym  : Symmetry.t ;
   expo : AO_expo.t ;
-}
+} with sexp
 
 let to_string p = 
   let { sym = s ; expo = e } = p in
@@ -13,3 +13,6 @@ let to_string p =
    (AO_expo.to_float e)
 ;;
   
+let of_sym_expo s e = 
+  { sym=s ; expo=e}
+;;

ocaml/Primitive.mli

@@ -0,0 +1,11 @@
+type t =
+{ sym : Symmetry.t;
+  expo : Qptypes.AO_expo.t;
+} with sexp
+
+(** Conversion to string for printing *)
+val to_string : t -> string
+
+(** Creation *)
+val of_sym_expo : Symmetry.t -> Qptypes.AO_expo.t -> t
+

ocaml/Qputils.ml

@@ -1,3 +1,6 @@
+open Core.Std;;
+
+(*
 let rec transpose = function
 | []          -> []
 | []::tail    -> transpose tail
@@ -7,5 +10,20 @@ let rec transpose = function
     in 
     new_head @ new_tail
 ;;
+*)
 
 
+let input_to_sexp s =
+    let result = 
+      String.split_lines s 
+      |> List.filter ~f:(fun x->
+         (String.strip x) <> "")
+      |> List.map ~f:(fun x->
+         "("^(String.tr '=' ' ' x)^")")
+      |> String.concat 
+    in
+    print_endline ("("^result^")");
+    "("^result^")"
+    |> Sexp.of_string
+;;
+

ocaml/README.rst

@@ -4,3 +4,9 @@ Ocaml scripts
 
 This directory contains all the scripts that control the input/output
 with the user.
+
+All executables start with `qp_` and all tests start with `test_`. Modules
+file names start with a capital letter.
+
+Info on how to extend the `qp_edit` tool is given in `README_qp_edit.rst`.
+

ocaml/README_qp_edit.rst

@@ -0,0 +1,220 @@
+Adding a new block
+==================
+
+In this section, we assume we will add the `New_keyword` keyword.
+
+Create the `Input_new_keyword.ml` file
+--------------------------------------
+
+Copy for example the `Input_full_ci.ml` file as a starting point.
+
+The template is the following, where `r_x`, `r_y`, ..., `last_r` are the records
+of the block.
+
+.. code-block:: ocaml
+
+  module New_keyword : sig
+    type t =
+      {  r_x    : Type_of_x.t
+         r_y    : Y_type.t
+        ...
+         last_r : bool
+      } with sexp
+    ;;
+    val read  : unit -> t
+    val write : t -> unit
+    val to_string : t -> string
+    val to_rst : t -> Rst_string.t
+    val of_rst : Rst_string.t -> t
+  end = struct
+    type t =
+      {  r_x    : Type_of_x.t
+         r_y    : Y_type.t
+        ...
+         last_r : bool
+      } with sexp
+    ;;
+
+    let get_default = Qpackage.get_ezfio_default "new_keyword";;
+
+
+    ...
+
+  end
+
+The following functions need to be defined::
+
+    val read   : unit -> t
+    val write  : t -> unit
+    val to_rst : t -> Rst_string.t
+    val of_rst : Rst_string.t -> t
+
+
+The type `t` has to be defined in a same way in the `sig` and the `struct`.
+
+For each record of the type `t`, use types defined in the `Qptypes.ml` file as
+much as possible.
+
+The `get_default` function will fetch the default values in the `ezfio_defaults` file
+in the `new_keyword` block.
+
+For each record `r_x` of the type `t`, create a `read_r_x ()` function
+and a `write_r_x r_x` function that performs the I/O in the EZFIO.
+To set a default value in the `read_r_x` function, use the following template
+(assuming that the `Type_of_x` is built from a `double precision` value in
+the EZFIO file).
+
+.. code-block:: ocaml
+
+  let read_r_x () =
+    if not (Ezfio.has_new_keyword_r_x ()) then
+       get_default "r_x"
+       |> Float.of_string
+       |> Ezfio.set_new_keyword_r_x
+    ;
+    Ezfio.get_new_keyword_r_x ()
+    |> Type_of_x.of_float
+  ;;
+
+  let write_r_x r_x =
+    Type_of_x.to_float r_x
+    |> Ezfio.set_new_keyword_r_x
+  ;;
+
+
+Then, create a `read` and a `write` function as
+
+.. code-block:: ocaml
+
+  let read () =
+    { r_x      = read_r_x () ;
+      r_y      = read_r_y () ;
+      ...
+      last_r   = read_last_r () ;
+    }
+  ;;
+
+  let write { r_x ;
+              r_y
+              ...
+              last_r ;
+            } = 
+    write_r_x r_x;
+    write_r_y r_y;
+    ...
+    write_last_r last_r;
+  ;;
+
+Finally, create the functions to write an RST string as
+
+.. code-block:: ocaml
+
+  let to_rst b =
+    Printf.sprintf "
+  You can put here some documentation as long as there is no equal sign.
+  The record entries should be indented on the right with  a blank line
+  before and a blank line after, as they would be in a rst file.
+
+  Here is the text for r_x
+
+    r_x = %s
+
+  And here is the text for r_y
+
+    r_y = %s
+
+  ...
+  Finally, the text for last_r
+
+    last_r = %s
+  "
+      (Type_of_x.to_string  b.r_x)
+      (Y_type.to_string     b.r_y)
+      ...
+      (Bool.to_string       b.last_r)
+  ;;
+
+
+and you can use this function to read it back:
+
+.. code-block:: ocaml
+
+  let of_rst s =
+    let s = Rst_string.to_string s
+    |> String.split ~on:'\n'
+    |> List.filter ~f:(fun line ->
+        String.contains line '=')
+    |> List.map ~f:(fun line ->
+        "("^(
+        String.tr line ~target:'=' ~replacement:' '
+        )^")" )
+    |> String.concat
+    in
+    Sexp.of_string ("("^s^")")
+    |> t_of_sexp
+  ;;
+
+
+
+  
+
+Add module to `Input.ml` file
+-----------------------------
+
+Append module to the `Input.ml` file. Use the name of the `Input_new_keyword.ml` without the
+`.ml` suffix.
+
+.. code-block:: ocaml
+
+  include Input_new_keyword;;
+
+
+In the `qp_edit.ml` file
+------------------------
+
+vim search strings are given in brackets.
+
+1. (`/type keyword`) : Add a new entry to the keyword type corresponding to the block to add:
+
+.. code-block:: ocaml
+
+  type keyword =
+  ...
+  | New_keyword
+  ;;
+
+
+
+2. (`/keyword_to_string`) : Add a new entry to the `keyword_to_string` function for the title of the block
+
+.. code-block:: ocaml
+
+  let keyword_to_string = function
+  ...
+  | New_keyword -> "My new keyword"
+  ;;
+
+
+3. (`/let get s`) : Add a new call to the to_rst function of the `Input.New_keyword` module
+
+.. code-block:: ocaml
+
+  let get s =
+    let header = (make_header s)
+      and rst = match s with
+      ...
+      | New_keyword ->
+        Input.New_keyword.(to_rst (read ()))
+      ...
+      
+
+4. (`/let set s`) : Add a new call to the of_rst function of the `Input.New_keyword` module
+
+.. code-block:: ocaml
+
+    match s with
+    ...
+    | New_keyword ->
+        Input.New_keyword.(write (of_rst str))
+    ;;  
+

ocaml/Range.ml

@@ -12,7 +12,7 @@ open Core.Std;;
 *)
 
 
-type t = int list ;;
+type t = int list with sexp
 
 let expand_range r =
   match String.lsplit2 ~on:'-' r with

ocaml/Range.mli

@@ -0,0 +1,10 @@
+type t = int list with sexp
+
+(** A range is a sorted list of ints in an interval.
+    It is created using a string :
+    "[a-b]" : range between a and b (included)
+    "[a]"   : the list with only one integer a
+    "a"     : equivalent to "[a]"
+ *)
+val of_string : string -> t
+val to_string : t -> string

ocaml/Symmetry.ml

@@ -1,7 +1,7 @@
 open Qptypes;;
 open Core.Std;;
 
-type t = S|P|D|F|G|H|I|J|K|L
+type t = S|P|D|F|G|H|I|J|K|L with sexp
 
 let to_string = function
   | S -> "S"
@@ -53,22 +53,31 @@ let to_l = function
   | J -> Positive_int.of_int 7
   | K -> Positive_int.of_int 8
   | L -> Positive_int.of_int 9
+;;
+
+let of_l i = 
+  let i = Positive_int.to_int i in
+  match i with
+  | 0 -> S
+  | 1 -> P
+  | 2 -> D
+  | 3 -> F
+  | 4 -> G
+  | 5 -> H
+  | 6 -> I
+  | 7 -> J
+  | 8 -> K
+  | 9 -> L
+  | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L"))
+;;
 
 type st = t
 ;;
 
-module Xyz : sig
+module Xyz = struct
   type t = { x: Positive_int.t ;
              y: Positive_int.t ;
-             z: Positive_int.t }
+             z: Positive_int.t } with sexp
-  val  of_string : string -> t
-  val  to_string : t -> string 
-  val  get_l     : t -> Positive_int.t
-  val  of_symmetry : st -> t list
-end = struct
-  type t = { x: Positive_int.t ;
-             y: Positive_int.t ;
-             z: Positive_int.t }
   type state_type = Null | X | Y | Z
 
   (** Builds an XYZ triplet from a string.
@@ -127,7 +136,9 @@ end = struct
     | 1 -> "z"
     | i -> Printf.sprintf "z%d" i
     in
-    x^y^z
+    let result = (x^y^z) in
+    if (result = "") then "s"
+    else result
   ;;
 
  (** Returns the l quantum number from a XYZ powers triplet *)
@@ -167,7 +178,10 @@ end = struct
    in
    create_x [] { x=(to_l sym) ; y=Positive_int.of_int 0 ;
      z=Positive_int.of_int 0 }
-   ;;
+ ;;
 
+ (** Returns the symmetry corresponding to the XYZ triplet *)
+ let to_symmetry sym = of_l (get_l sym)
+ ;;
 end
 

ocaml/Symmetry.mli

@@ -0,0 +1,36 @@
+type t = S | P | D | F | G | H | I | J | K | L with sexp
+
+(** Creatio from strings *)
+val to_string : t -> string
+val of_string : string -> t
+val of_char : char -> t
+
+(** Connexion with l quantum number *)
+val to_l : t -> Qptypes.Positive_int.t
+val of_l : Qptypes.Positive_int.t -> t
+
+type st = t
+module Xyz :
+  sig
+    type t = {
+      x : Qptypes.Positive_int.t;
+      y : Qptypes.Positive_int.t;
+      z : Qptypes.Positive_int.t;
+    } with sexp
+
+    (** The string format contains the powers of x,y and z in a
+        format like "x2z3" *)
+
+    val of_string : string -> t
+    val to_string : t -> string
+
+    (** Returns the quantum number l *)
+    val get_l : t -> Qptypes.Positive_int.t
+
+    (** Returns a list of XYZ powers for a given symmetry *)
+    val of_symmetry : st -> t list
+
+    (** Returns the symmetry corresponding to the XYZ powers *)
+    val to_symmetry : t -> st
+
+  end

ocaml/To_md5.ml

@@ -0,0 +1,11 @@
+open Core.Std;;
+open Qptypes;;
+
+let to_md5 sexp_of_t t =
+  sexp_of_t t
+  |> Sexp.to_string
+  |> Cryptokit.hash_string (Cryptokit.Hash.md5 ())
+  |> Cryptokit.transform_string (Cryptokit.Hexa.encode ())
+  |> MD5.of_string
+;;
+

ocaml/Units.mli

@@ -0,0 +1,7 @@
+type units = Bohr | Angstrom
+
+(** Conversion functions *)
+val angstrom_to_bohr : float
+val bohr_to_angstrom : float
+
+

ocaml/_tags

@@ -1,3 +1,2 @@
-true: package(core)
+true: package(core,sexplib.syntax,cryptokit)
 true: thread
-

ocaml/qp_create_ezfio_from_xyz.ml

@@ -26,8 +26,8 @@ let run ?o b c m xyz_file =
 
   (* Read molecule *)
   let molecule =
-    Molecule.of_xyz_file xyz_file ~charge:c
+    (Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c)
-      ~multiplicity:(Multiplicity.of_int m)
+      ~multiplicity:(Multiplicity.of_int m) )
   in
 
   (* Build EZFIO File name *)
@@ -74,9 +74,10 @@ let run ?o b c m xyz_file =
 
   (* Write Basis set *)
   let basis =
+    let nmax = Nucl_number.get_max () in
     let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
     | [] -> accu
-    | e::tail -> let new_accu = (e,(Nucl_number.of_int n))::accu in
+    | e::tail -> let new_accu = (e,(Nucl_number.of_int ~max:nmax n))::accu in
       do_work new_accu (n+1) tail
     in
     do_work [] 1  nuclei

ocaml/qp_edit.ml

@@ -0,0 +1,209 @@
+open Qputils;;
+open Qptypes;;
+open Core.Std;;
+
+let file_header filename = Printf.sprintf
+"
+==================================================================
+                       Quantum Package
+==================================================================
+
+Editing file `%s`
+
+" filename
+  
+type keyword = 
+| Ao_basis
+| Bielec_integrals
+| Cisd_sc2
+| Determinants
+| Electrons
+| Full_ci
+| Hartree_fock
+| Mo_basis
+| Nuclei
+;;
+
+let keyword_to_string = function
+| Ao_basis         -> "AO basis"
+| Bielec_integrals -> "Two electron integrals"
+| Cisd_sc2         -> "CISD (SC)^2"
+| Determinants     -> "Determinants"
+| Electrons        -> "Electrons"
+| Full_ci          -> "Selected Full-CI"
+| Hartree_fock     -> "Hartree-Fock"
+| Mo_basis         -> "MO basis"
+| Nuclei           -> "Molecule"
+;;
+
+let make_header kw =
+  let s = keyword_to_string kw in
+  let l = String.length s in
+  "\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n"
+;;
+
+let get s = 
+  let header = (make_header s) 
+  and rst = match s with
+  | Full_ci ->
+    Input.Full_ci.(to_rst (read ()))
+  | Hartree_fock ->
+    Input.Hartree_fock.(to_rst (read ()))
+  | Mo_basis ->
+    Input.Mo_basis.(to_rst (read ()))
+  | Electrons ->
+    Input.Electrons.(to_rst (read ()))
+  | Determinants ->
+    Input.Determinants.(to_rst (read ()))
+  | Cisd_sc2 ->
+    Input.Cisd_sc2.(to_rst (read ()))
+  | Nuclei ->
+    Input.Nuclei.(to_rst (read ()))
+  | Ao_basis ->
+    Input.Ao_basis.(to_rst (read ()))
+  | Bielec_integrals -> 
+    Input.Bielec_integrals.(to_rst (read ()))
+ 
+  in header^(Rst_string.to_string rst)
+;;
+
+let set str s = 
+  let header = (make_header s) in
+  let index_begin = String.substr_index_exn ~pos:0 ~pattern:header str in
+  let index_begin = index_begin + (String.length header) in
+  let index_end   = 
+    match ( String.substr_index ~pos:(index_begin+(String.length header)+1)
+      ~pattern:"==" str) with
+      | Some i -> i
+      | None -> String.length str
+  in
+  let l = index_end - index_begin in
+  let str = String.sub ~pos:index_begin ~len:l str
+  |> Rst_string.of_string
+  in
+  match s with
+  (*
+  | Mo_basis ->
+    *)
+  | Hartree_fock ->
+      Input.Hartree_fock.(write (of_rst str ))
+  | Full_ci ->
+      Input.Full_ci.(write (of_rst str))
+  | Electrons ->
+      Input.Electrons.(write (of_rst str))
+  | Determinants ->
+      Input.Determinants.(write (of_rst str))
+  | Cisd_sc2 ->
+      Input.Cisd_sc2.(write (of_rst str))
+  | Nuclei ->
+      Input.Nuclei.(write (of_rst str))
+  | Bielec_integrals -> 
+      Input.Bielec_integrals.(write (of_rst str))
+    (*
+  | Ao_basis ->
+    *)
+ 
+;;
+
+
+let create_temp_file ezfio_filename fields =
+  let temp_filename  = Filename.temp_file "qp_edit_" ".rst" in
+  Out_channel.with_file temp_filename ~f:(fun out_channel ->
+    (file_header ezfio_filename) :: (List.map ~f:get fields) 
+    |> String.concat ~sep:"\n" 
+    |> Out_channel.output_string out_channel 
+  );
+  temp_filename
+;;
+
+let run ezfio_filename =
+
+  (* Open EZFIO *)
+  if (not (Sys.file_exists_exn ezfio_filename)) then
+    failwith (ezfio_filename^" does not exists");
+
+  Ezfio.set_file ezfio_filename;
+
+  (*
+  let output = (file_header ezfio_filename) :: (
+    List.map ~f:get [
+      Ao_basis ; 
+      Mo_basis ; 
+    ])
+   in
+  String.concat output
+  |> print_string
+  *)
+  
+  let tasks = [
+      Nuclei ;
+      Electrons ;
+      Bielec_integrals ;
+      Hartree_fock ;
+      Cisd_sc2 ;
+      Full_ci ;
+      Determinants ;
+  ]
+  in
+
+  (* Create the temp file *)
+  let temp_filename = create_temp_file ezfio_filename tasks in
+
+  (* Open the temp file with external editor *)
+  let editor = 
+    match Sys.getenv "EDITOR" with
+    | Some editor -> editor
+    | None -> "vi"
+  in
+  let command = Printf.sprintf "%s %s" editor temp_filename in
+  Sys.command_exn command;
+
+  (* Re-read the temp file *)
+  let temp_string  =
+    In_channel.with_file temp_filename ~f:(fun in_channel ->
+      In_channel.input_all in_channel) 
+  in
+  List.iter ~f:(fun x -> set temp_string x) tasks;
+
+  (* Remove temp_file *)
+  Sys.remove temp_filename;
+;;
+
+
+
+let spec =
+  let open Command.Spec in
+  empty 
+(*
+  +> flag "i"  (optional string)
+     ~doc:"Prints input data"
+  +> flag "o" (optional string)
+     ~doc:"Prints output data"
+*)
+  +> anon ("ezfio_file" %: string)
+;;
+
+let command = 
+    Command.basic 
+    ~summary: "Quantum Package command"
+    ~readme:(fun () ->
+      "
+Edit input data
+      ")
+    spec
+(*    (fun i o ezfio_file () -> *)
+    (*fun ezfio_file () -> 
+       try 
+           run ezfio_file
+       with
+       | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n")
+    *)
+    (fun ezfio_file () -> run ezfio_file)
+;;
+
+let () =
+    Command.run command
+;;
+
+
+

ocaml/qp_print.ml

@@ -52,13 +52,14 @@ let run_i ~action ezfio_filename =
 
   let compute_charge () =
     let input = Input.Electrons.read () in
-    let nucl_charge = Ezfio.((get_nuclei_nucl_charge ()).data)
+    let nucl_charge = Ezfio.get_nuclei_nucl_charge ()
-      |> Ezfio.flattened_ezfio_data |> Array.map ~f:(Float.to_int) 
+      |> Ezfio.flattened_ezfio |> Array.map ~f:(Float.to_int) 
     and n_alpha = input.Input.Electrons.elec_alpha_num
       |> Elec_alpha_number.to_int
     and n_beta  = input.Input.Electrons.elec_beta_num 
       |> Elec_beta_number.to_int
     in Array.fold ~init:(-n_alpha-n_beta) ~f:(fun x y -> x+y) nucl_charge
+    |> Charge.of_int
   in
 
   let compute_multiplicity () =
@@ -70,17 +71,17 @@ let run_i ~action ezfio_filename =
 
   let create_molecule () =
     let nucl_num = Ezfio.get_nuclei_nucl_num () 
-    and nucl_charge = Ezfio.((get_nuclei_nucl_charge ()).data)
+    and nucl_charge = Ezfio.get_nuclei_nucl_charge ()
-       |> Ezfio.flattened_ezfio_data
+       |> Ezfio.flattened_ezfio
-    and nucl_coord = Ezfio.((get_nuclei_nucl_coord ()).data )
+    and nucl_coord = Ezfio.get_nuclei_nucl_coord ()
-       |> Ezfio.flattened_ezfio_data
+       |> Ezfio.flattened_ezfio
     in
     let nucl_label =
       if (Ezfio.has_nuclei_nucl_label ()) then
-       Ezfio.((get_nuclei_nucl_label ()).data) |> Ezfio.flattened_ezfio_data
+       Ezfio.get_nuclei_nucl_label () |> Ezfio.flattened_ezfio
       else
         Array.map ~f:(fun x-> x
-          |> Float.to_int
+          |> Charge.of_float
           |> Element.of_charge 
           |> Element.to_string ) nucl_charge
     in

ocaml/qp_set_ddci.ml

@@ -0,0 +1,303 @@
+open Qputils;;
+open Qptypes;;
+open Core.Std;;
+
+(*
+ * Command-line arguments
+ * ----------------------
+ *)
+
+let build_mask from upto n_int =
+  let from  = MO_number.to_int from
+  and upto  = MO_number.to_int upto
+  and n_int = N_int_number.to_int n_int
+  in
+  let rec build_mask bit = function
+    | 0 -> []
+    | i -> 
+        if ( i = upto ) then
+          Bit.One::(build_mask Bit.One (i-1))
+        else if ( i = from ) then
+          Bit.One::(build_mask Bit.Zero (i-1))
+        else
+          bit::(build_mask bit (i-1))
+  in
+  let starting_bit = 
+    if ( (upto >= n_int*64) || (upto < 0) ) then Bit.One
+    else Bit.Zero
+  in
+  build_mask starting_bit (n_int*64)
+  |> List.rev
+;;
+
+
+let failure s = raise (Failure s)
+;;
+
+type t = 
+  | Core
+  | Inactive
+  | Active
+  | Virtual
+  | Deleted
+  | None
+;;
+
+let t_to_string = function
+  | Core -> "core"
+  | Inactive -> "inactive"
+  | Active -> "active"
+  | Virtual -> "virtual"
+  | Deleted -> "deleted"
+  | None -> assert false
+;;
+
+let run ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_filename =
+
+  Ezfio.set_file ezfio_filename ;
+  if not (Ezfio.has_mo_basis_mo_tot_num ()) then
+    failure "mo_basis/mo_tot_num not found" ;
+
+  let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in
+  let n_int =
+     try  N_int_number.of_int (Ezfio.get_determinants_n_int ())
+     with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num 
+  in
+
+
+  let mo_class = Array.init mo_tot_num ~f:(fun i -> None) in
+
+  (* Check input data *)
+  let apply_class l = 
+    let rec apply_class t = function
+    | [] -> ()
+    | k::tail -> let i = MO_number.to_int k in
+        begin
+          match mo_class.(i-1) with
+          | None -> mo_class.(i-1) <- t ;
+            apply_class t tail;
+          | x -> failure 
+             (Printf.sprintf "Orbital %d is defined both in the %s and %s spaces"
+             i (t_to_string x) (t_to_string t))
+        end
+    in
+    match l with
+    | MO_class.Core     x -> apply_class Core      x
+    | MO_class.Inactive x -> apply_class Inactive  x
+    | MO_class.Active   x -> apply_class Active    x
+    | MO_class.Virtual  x -> apply_class Virtual   x
+    | MO_class.Deleted  x -> apply_class Deleted   x
+  in
+
+  let core_input = core in
+  let core  = MO_class.create_core     core in
+  let inact = MO_class.create_inactive inact in
+  let act   = MO_class.create_active   act in
+  let virt  = MO_class.create_virtual  virt in
+  let del   = MO_class.create_deleted  del in
+
+  apply_class core  ;
+  apply_class inact ;
+  apply_class act   ;
+  apply_class virt  ;
+  apply_class del   ;
+
+  for i=1 to (Array.length mo_class)
+  do
+    if (mo_class.(i-1) = None) then
+      failure (Printf.sprintf "Orbital %d is not specified (mo_tot_num = %d)" i mo_tot_num)
+  done;
+  
+  
+  (* Debug output *)
+  MO_class.to_string core  |> print_endline ;
+  MO_class.to_string inact |> print_endline ;
+  MO_class.to_string act   |> print_endline ;
+  MO_class.to_string virt  |> print_endline ;
+  MO_class.to_string del   |> print_endline ;
+
+  (* Create masks *)
+  let ia = Excitation.create_single inact act 
+  and aa = Excitation.create_single act act 
+  and av = Excitation.create_single act virt
+  and iv = Excitation.create_single inact virt
+  in
+  let single_excitations = [| ia ; aa ; av ; iv |]
+    |> Array.map ~f:Excitation.(fun x ->
+       match x with
+       | Single (x,y) -> 
+         ( MO_class.to_bitlist n_int (Hole.to_mo_class x),
+           MO_class.to_bitlist n_int (Particle.to_mo_class y) ) 
+       | Double _ -> assert false
+       )
+       
+  and double_excitations = [|
+    Excitation.double_of_singles ia ia ;
+    Excitation.double_of_singles ia aa ;
+    Excitation.double_of_singles ia iv ;
+    Excitation.double_of_singles ia av ;
+
+    Excitation.double_of_singles aa aa ;
+    Excitation.double_of_singles aa iv ;
+    Excitation.double_of_singles aa av ;
+
+    Excitation.double_of_singles iv aa ;
+    Excitation.double_of_singles iv av ;
+
+(*    Excitation.double_of_singles iv iv ; *)
+ |]
+
+    |> Array.map ~f:Excitation.(fun x ->
+       match x with
+       | Single _ -> assert false
+       | Double (x,y,z,t) -> 
+         ( MO_class.to_bitlist n_int (Hole.to_mo_class x),
+           MO_class.to_bitlist n_int (Particle.to_mo_class y) , 
+           MO_class.to_bitlist n_int (Hole.to_mo_class z),
+           MO_class.to_bitlist n_int (Particle.to_mo_class t) )
+       )
+  in
+  
+  let extract_hole (h,_) = h 
+  and extract_particle (_,p) = p 
+  and extract_hole1 (h,_,_,_) = h 
+  and extract_particle1 (_,p,_,_) = p 
+  and extract_hole2 (_,_,h,_) = h 
+  and extract_particle2 (_,_,_,p) = p 
+  in
+  let result_ref = 
+    let core = MO_class.create_inactive core_input in
+    let cv = Excitation.create_single core virt in
+    let cv = match cv with
+    | Excitation.Single (x,y) ->  
+      ( MO_class.to_bitlist n_int (Excitation.Hole.to_mo_class x),
+        MO_class.to_bitlist n_int (Excitation.Particle.to_mo_class y) ) 
+    | Excitation.Double _ -> assert false
+    in
+    let iv = match iv with
+    | Excitation.Single (x,y) ->  
+      ( MO_class.to_bitlist n_int (Excitation.Hole.to_mo_class x),
+        MO_class.to_bitlist n_int (Excitation.Particle.to_mo_class y) ) 
+    | Excitation.Double _ -> assert false
+    in
+    [ Bitlist.or_operator (extract_hole iv) (extract_hole cv);
+      extract_particle iv ]
+  in
+
+  let n_single = Array.length single_excitations in
+  let n_mask = Array.length double_excitations in
+  let zero = List.init (N_int_number.to_int n_int) ~f:(fun i -> 0L)
+   |> Bitlist.of_int64_list in
+  let result_gen = (List.init n_single ~f:(fun i-> [
+    extract_hole      single_excitations.(i) ; 
+    extract_particle  single_excitations.(i) ; 
+    extract_hole1     double_excitations.(i) ; 
+    extract_particle1 double_excitations.(i) ; 
+    extract_hole2     double_excitations.(i) ; 
+    extract_particle2 double_excitations.(i) ; ])
+   )@(List.init (n_mask-n_single) ~f:(fun i-> [
+    zero ; zero ;
+    extract_hole1     double_excitations.(n_single+i) ; 
+    extract_particle1 double_excitations.(n_single+i) ; 
+    extract_hole2     double_excitations.(n_single+i) ; 
+    extract_particle2 double_excitations.(n_single+i) ; ])
+   )
+   |> List.concat
+  in
+
+  (* Print bitmasks *)
+  print_endline "Reference  Bitmasks:";
+  List.iter ~f:(fun x-> print_endline (Bitlist.to_string x)) result_ref;
+
+  print_endline "Generators Bitmasks:";
+  List.iter ~f:(fun x-> print_endline (Bitlist.to_string x)) result_gen;
+
+  (* Transform to int64 *)
+  let result_gen =  List.map ~f:(fun x ->
+     let y = Bitlist.to_int64_list x in y@y )
+     result_gen 
+  |> List.concat
+  in
+  let result_ref =  List.map ~f:(fun x ->
+     let y = Bitlist.to_int64_list x in y@y )
+     result_ref 
+  |> List.concat
+  in
+
+  (* Write generators masks *)
+  Ezfio.set_bitmasks_n_int (N_int_number.to_int n_int);
+  Ezfio.set_bitmasks_bit_kind 8;
+  Ezfio.set_bitmasks_n_mask_gen n_mask;
+  Ezfio.ezfio_array_of_list ~rank:4 ~dim:([| (N_int_number.to_int n_int) ; 2; 6; n_mask|]) ~data:result_gen
+  |> Ezfio.set_bitmasks_generators ; 
+
+  (* Write reference masks *)
+  Ezfio.set_bitmasks_n_mask_ref 1;
+  Ezfio.ezfio_array_of_list ~rank:4 ~dim:([| (N_int_number.to_int n_int) ; 2; 2; 1|]) ~data:result_ref
+  |> Ezfio.set_bitmasks_reference ; 
+
+
+;;
+
+let ezfio_file =
+  let failure filename = 
+        eprintf "'%s' is not an EZFIO file.\n%!" filename;
+        exit 1
+  in
+  Command.Spec.Arg_type.create
+  (fun filename ->
+    match Sys.is_directory filename with
+    | `Yes -> 
+        begin
+          match Sys.is_file (filename ^ "/.version") with
+          | `Yes -> filename
+          | _ -> failure filename
+        end
+    | _ -> failure filename
+  )
+;;
+
+let default range =
+  let failure filename = 
+        eprintf "'%s' is not a regular file.\n%!" filename;
+        exit 1
+  in
+  Command.Spec.Arg_type.create
+  (fun filename ->
+    match Sys.is_directory filename with
+    | `Yes -> 
+        begin
+          match Sys.is_file (filename^"/.version") with
+          | `Yes -> filename
+          | _ -> failure filename
+        end
+    | _ -> failure filename
+  )
+;;
+
+let spec =
+  let open Command.Spec in
+  empty 
+  +> flag "core"   (optional string) ~doc:"range Range of core orbitals"
+  +> flag "inact"  (optional string) ~doc:"range Range of inactive orbitals"
+  +> flag "act"    (optional string) ~doc:"range Range of active orbitals"
+  +> flag "virt"   (optional string) ~doc:"range Range of virtual orbitals"
+  +> flag "del"    (optional string) ~doc:"range Range of deleted orbitals"
+  +> anon ("ezfio_filename" %: ezfio_file)
+;;
+
+let command = 
+    Command.basic 
+    ~summary: "Quantum Package command"
+    ~readme:(fun () ->
+     "Set the orbital classes in an EZFIO directory
+      The range of MOs has the form : \"[36-53,72-107,126-131]\"
+        ")
+    spec
+    (fun core inact act virt del ezfio_filename () -> run ?core ?inact ?act ?virt ?del ezfio_filename )
+;;
+
+let () =
+    Command.run command
+
+

ocaml/qptypes_generator.ml

@@ -36,40 +36,11 @@ let input_data = "
 * Non_empty_string : string
   assert (x <> \"\") ;
 
-* MO_number : int  
-  assert (x > 0) ; 
-  if (x > 1000) then
-    warning \"More than 1000 MOs\";
-  if (Ezfio.has_mo_basis_mo_tot_num ()) then
-    assert (x <= (Ezfio.get_mo_basis_mo_tot_num ()));
 
-* AO_number : int  
+* Det_number_max : int 
-  assert (x > 0) ; 
-  if (x > 1000) then
-    warning \"More than 1000 AOs\";
-  if (Ezfio.has_ao_basis_ao_num ()) then
-    assert (x <= (Ezfio.get_ao_basis_ao_num ()));
-
-* Nucl_number : int  
-  assert (x > 0) ; 
-  if (x > 1000) then
-    warning \"More than 1000 atoms\";
-  if (Ezfio.has_nuclei_nucl_num ()) then
-    assert (x <= (Ezfio.get_nuclei_nucl_num ()));
-
-* N_int_number : int 
-  assert (x > 0) ; 
-  if (x > 100) then
-    warning \"N_int > 100\";
-  if (Ezfio.has_determinants_n_int ()) then
-    assert (x = (Ezfio.get_determinants_n_int ()));
-
-* Det_number : int 
   assert (x > 0) ; 
   if (x > 100000000) then
     warning \"More than 100 million determinants\";
-  if (Ezfio.has_determinants_det_num ()) then
-    assert (x <= (Ezfio.get_determinants_det_num ()));
 
 * States_number : int 
   assert (x > 0) ; 
@@ -97,6 +68,9 @@ let input_data = "
 
 * MO_coef : float
 
+* MO_occ : float
+  assert (x >= 0.); 
+
 * AO_coef : float
 
 * AO_expo : float  
@@ -121,42 +95,56 @@ let input_data = "
 * Elec_number : int
   assert (x > 0) ;
 
+* MD5 : string
+  assert ((String.length x) = 32);
+
+* Rst_string : string
+
+"
+;;
+
+let input_ezfio = "
+* MO_number : int
+  mo_basis_mo_tot_num
+  1 : 10000
+  More than 10000 MOs
+
+* AO_number : int
+  ao_basis_ao_num
+  1 : 10000
+  More than 10000 AOs
+
+* Nucl_number : int
+  nuclei_nucl_num
+  1 : 10000
+  More than 10000 nuclei
+
+* N_int_number : int
+  determinants_n_int
+  1 : 30
+  N_int > 30
+
+* Det_number : int
+  determinants_n_det
+  1 : 100000000
+  More than 100 million determinants
+
 "
 ;;
 
 let untouched = "
-module Determinant : sig
-  type t
-  val to_int64_array : t -> int64 array
-  val of_int64_array : int64 array -> t
-  val to_string : t -> string
-end = struct
-  type t = int64 array
-  let to_int64_array x = x
-  let of_int64_array x = 
-    if (Ezfio.has_determinants_n_int ()) then
-      begin
-        let n_int = Ezfio.get_determinants_n_int () in
-        assert ((Array.length x) = n_int*2) 
-      end
-      ; x
-  let to_string x = Array.to_list x
-    |> List.map ~f:Int64.to_string
-    |> String.concat ~sep:\", \"
-end
-
 "
 
 let template = format_of_string "
 module %s : sig
-  type t
+  type t with sexp
   val to_%s : t -> %s
-  val of_%s : %s -> t
+  val of_%s : %s %s -> t
   val to_string : t -> string
 end = struct
-  type t = %s
+  type t = %s with sexp
   let to_%s x = x
-  let of_%s x = ( %s x )
+  let of_%s %s x = ( %s x )
   let to_string x = %s.to_string x
 end
 
@@ -169,13 +157,18 @@ let parse_input input=
     | [] -> result
     | ( "" , ""   )::tail -> parse result tail
     | ( t  , text )::tail -> 
-        let  name , typ  = String.lsplit2_exn ~on:':' t
+        let name,typ,params,params_val = 
+          match String.split ~on:':' t with
+          | [name;typ] -> (name,typ,"","")
+          | name::typ::params::params_val -> (name,typ,params,
+            (String.concat params_val ~sep:":") )
+          | _ -> assert false
         in
         let typ  = String.strip typ
         and name = String.strip name in
         let typ_cap = String.capitalize typ in
-        let newstring = Printf.sprintf template name typ typ typ typ typ typ typ 
+        let newstring = Printf.sprintf template name typ typ typ params_val typ typ 
-          ( String.strip text ) typ_cap
+          typ typ params ( String.strip text ) typ_cap
         in
         List.rev (parse (newstring::result) tail )
   in
@@ -186,9 +179,76 @@ let parse_input input=
   |> print_string
 ;;
 
-let () = 
+
-  parse_input input_data ;
+let ezfio_template = format_of_string "
-  print_endline untouched
+module %s : sig
+  type t with 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
+  let to_string x = %s.to_string x
+  let get_max () =
+    if (Ezfio.has_%s ()) then
+      Ezfio.get_%s ()
+    else
+      %s
+  let get_min () =
+      %s
+  let to_%s x = x
+  let of_%s ?(min=get_min ()) ?(max=get_max ()) x = 
+    begin
+      assert (x >= min) ;
+      if (x > %s) then
+        warning \"%s\";
+      begin
+        match max with
+        | %s -> ()
+        | i  -> assert ( x <= i )
+      end ;
+      x
+    end
+end
+"
 ;;
 
+let parse_input_ezfio input=
+  let parse s = 
+    match (
+      String.split s ~on:'\n'
+      |> List.filter ~f:(fun x -> (String.strip x) <> "")
+    ) with
+    | [] -> ""
+    | a :: b :: c :: d :: [] ->
+      begin
+        let (name,typ) = String.lsplit2_exn ~on:':' a
+        and ezfio_func = b
+        and (min, max) = String.lsplit2_exn ~on:':' c
+        and msg = d
+        in 
+        let name :: typ :: ezfio_func :: min :: max :: msg :: [] = 
+        match (name :: typ :: ezfio_func :: min :: max :: msg :: []) with
+        | l -> List.map ~f:String.strip  l
+        | _ -> assert false
+        in
+        Printf.sprintf ezfio_template 
+          name typ typ typ typ typ typ typ typ (String.capitalize typ)
+          ezfio_func ezfio_func max min typ typ max msg min
+      end
+    | _ -> failwith "Error in input_ezfio"
+  in
+     String.split ~on:'*' input
+  |> List.map ~f:parse
+  |> String.concat 
+  |> print_string
+;;
+
+let () = 
+  parse_input input_data ;
+  parse_input_ezfio input_ezfio;
+  print_endline untouched;
+
+
 

ocaml/test_basis.ml

@@ -21,9 +21,21 @@ let test_module () =
       (Basis.read_element basis_channel (Nucl_number.of_int 2) Element.F) 
     in
 
-    Long_basis.of_basis basis 
+    print_string "Long basis\n==========\n";
-    |> Long_basis.to_string 
+    let long_basis = 
-    |> print_endline
+      Long_basis.of_basis basis 
+    in
+    print_endline (Long_basis.to_string long_basis);
+
+    let short_basis = 
+      Long_basis.to_basis long_basis
+    in
+    if (short_basis <> basis) then
+      print_endline "(short_basis <> basis)"
+    ;
+    print_string "Short basis\n===========\n";
+    print_endline (Basis.to_string basis);
+    print_endline ("MD5: "^(Basis.to_md5 basis |> MD5.to_string));
 ;;
 
 test_module ();

ocaml/test_determinants.ml

@@ -0,0 +1,15 @@
+open Qptypes;;
+
+let test_module () = 
+  let mo_tot_num = MO_number.of_int 10 in
+  let det =
+   [| 15L ; 7L |]
+   |> Determinant.of_int64_array 
+      ~n_int:(N_int_number.of_int 1)
+      ~alpha:(Elec_alpha_number.of_int 4)
+      ~beta:(Elec_beta_number.of_int 3)
+  in
+  Printf.printf "%s\n" (Determinant.to_string (~mo_tot_num:mo_tot_num) det)
+;;
+
+test_module ();;

ocaml/test_elements.ml

@@ -1,6 +1,6 @@
 let test_module () = 
   let atom = Element.of_string "Cobalt" in
-  Printf.printf "%s %d\n" (Element.to_string atom) (Element.to_charge atom)
+  Printf.printf "%s %d\n" (Element.to_string atom) (Charge.to_int (Element.to_charge atom))
 ;;
 
 test_module ();;

ocaml/test_gto.ml

@@ -13,11 +13,20 @@ let test_gto_1 () =
   let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in
   ignore (input_line in_channel);
   let gto = Gto.read_one in_channel in
-  print_string (Gto.to_string gto);
+  print_endline (Gto.to_string gto);
-  let gto = Gto.read_one in_channel in
+  In_channel.seek in_channel 0L;
-  print_string (Gto.to_string gto);
+  ignore (input_line in_channel);
-  let gto = Gto.read_one in_channel in
+  let gto2 = Gto.read_one in_channel in
-  print_string (Gto.to_string gto);
+  print_endline (Gto.to_string gto2);
+  let gto3 = Gto.read_one in_channel in
+  print_endline (Gto.to_string gto3);
+  if (gto2 = gto) then
+    print_endline "gto2 = gto";
+  if (gto3 = gto) then
+    print_endline "gto3 = gto";
+  if (gto3 = gto3) then
+    print_endline "gto3 = gto3";
+  
 ;;
 
 let test_gto_2 () =
@@ -34,7 +43,7 @@ let test_gto () =
 ;;
 
 let test_module () =
-  test_gto()
+  test_gto_1()
 ;;
 
 test_module ();;

ocaml/test_input.ml

@@ -1,15 +1,26 @@
+open Qptypes;;
+
 let test_ao () =
   Ezfio.set_file "F2.ezfio" ;
   let b = Input.Ao_basis.read ()
   in
   print_endline (Input.Ao_basis.to_string b);
+  print_endline (Input.Ao_basis.to_rst b |> Rst_string.to_string);
 ;;
 
 let test_bielec_intergals () =
   Ezfio.set_file "F2.ezfio" ;
   let b = Input.Bielec_integrals.read ()
   in
-  print_endline (Input.Bielec_integrals.to_string b);
+  let output = Input.Bielec_integrals.to_string b
+  in
+  print_endline output;
+  let rst = Input.Bielec_integrals.to_rst b in
+  let b2 = Input.Bielec_integrals.of_rst rst in
+  if (b = b2) then
+    print_endline "OK"
+  else
+    print_endline "rst failed";
 ;;
 
 let test_bitmasks () =
@@ -30,7 +41,14 @@ let test_dets () =
   Ezfio.set_file "F2.ezfio" ;
   let b = Input.Determinants.read ()
   in
-  print_endline (Input.Determinants.to_string b);
+  print_endline (Input.Determinants.to_rst b |> Rst_string.to_string ) ;
+  print_endline (Input.Determinants.sexp_of_t b |> Sexplib.Sexp.to_string ) ;
+  let r = Input.Determinants.to_rst b in
+  let b2 = Input.Determinants.of_rst r in
+  if (b2 = b) then
+    print_endline "OK"
+  else
+    print_endline "Failed"
 ;;
 
 let test_cisd_sc2 () =
@@ -38,6 +56,13 @@ let test_cisd_sc2 () =
   let b = Input.Cisd_sc2.read ()
   in
   print_endline (Input.Cisd_sc2.to_string b);
+  let rst = Input.Cisd_sc2.to_rst b in
+  let b2 = Input.Cisd_sc2.of_rst rst in
+  if (b = b2) then
+    print_endline "OK"
+  else
+    print_endline "rst failed";
+
 ;;
 
 let test_electrons () =
@@ -45,6 +70,12 @@ let test_electrons () =
   let b = Input.Electrons.read ()
   in
   print_endline (Input.Electrons.to_string b);
+  let rst = Input.Electrons.to_rst b in
+  let new_b = Input.Electrons.of_rst rst in
+  if (b = new_b) then
+    print_endline "OK"
+  else
+    print_endline "Failed in rst"
 ;;
 
 let test_fci () =
@@ -52,6 +83,13 @@ let test_fci () =
   let b = Input.Full_ci.read ()
   in
   print_endline (Input.Full_ci.to_string b);
+  let rst = Input.Full_ci.to_rst b in
+  let new_b = Input.Full_ci.of_rst rst in
+  print_endline (Input.Full_ci.to_string b);
+  if (b = new_b) then
+    print_endline "OK"
+  else
+    print_endline "Failed in rst"
 ;;
 
 let test_hf () =
@@ -59,6 +97,13 @@ let test_hf () =
   let b = Input.Hartree_fock.read ()
   in
   print_endline (Input.Hartree_fock.to_string b);
+  let rst = Input.Hartree_fock.to_rst b in
+  let new_b = Input.Hartree_fock.of_rst rst in
+  print_endline (Input.Hartree_fock.to_string b);
+  if (b = new_b) then
+    print_endline "OK"
+  else
+    print_endline "Failed in rst"
 ;;
 
 let test_mo () =
@@ -70,19 +115,27 @@ let test_mo () =
 
 let test_nucl () =
   Ezfio.set_file "F2.ezfio" ;
-  let b = Input.Nuclei.read ()
+  let b = Input.Nuclei.read () in
-  in
+  let rst = Input.Nuclei.to_rst b in
+  let new_b = Input.Nuclei.of_rst rst in
   print_endline (Input.Nuclei.to_string b);
+  if (b = new_b) then
+    print_endline "OK"
+  else
+    print_endline "Failed in rst"
 ;;
 
 (*
-test_hf ();;
 test_ao ();;
-test_bielec_intergals ();;
 test_bitmasks ();
 test_cis ();
-test_dets ();
 test_cisd_sc2 ();
+test_dets ();
+test_hf ();;
 test_mo ();;
-*)
 test_nucl ();
+test_bielec_intergals ();;
+test_electrons();
+*)
+test_dets ();
+

ocaml/test_mo_label.ml

@@ -0,0 +1,5 @@
+let () =
+  let m = MO_label.of_string "canonical" in
+  let s = MO_label.to_string m in
+  print_string s
+;;

ocaml/test_molecule.ml

@@ -21,7 +21,7 @@ H           1.0       0.54386314      0.00000000     -0.92559535
   print_string "---\n";
   let m = Molecule.of_xyz_string xyz 
   in print_endline (Molecule.name m) ;
-  let m = Molecule.of_xyz_string xyz ~charge:1 ~multiplicity:(Multiplicity.of_int 2)
+  let m = Molecule.of_xyz_string xyz ~charge:(Charge.of_int 1) ~multiplicity:(Multiplicity.of_int 2)
   in print_endline (Molecule.name m) ;
 
   let xyz =
@@ -31,7 +31,7 @@ O        1.65102147      0.00000000     -2.35602344
 H        0.54386314      0.00000000     -0.92559535
 "
   in
-  let m = Molecule.of_xyz_string xyz ~charge:(-2) 
+  let m = Molecule.of_xyz_string xyz ~charge:(Charge.of_int (-2))
   in print_endline (Molecule.name m) ;
   print_endline (Molecule.to_string m);
   print_string "---------\n";

scripts/generate_h_apply.py

@@ -222,6 +222,9 @@ class H_apply(object):
       self.data["size_max"] = str(1024*128) 
       self.data["copy_buffer"] = """
       call copy_H_apply_buffer_to_wf
+      if (s2_eig) then
+        call make_s2_eigenfunction
+      endif
       SOFT_TOUCH psi_det psi_coef N_det
       selection_criterion_min = min(selection_criterion_min, maxval(select_max))*0.1d0
       selection_criterion = selection_criterion_min

scripts/install_ocaml.sh

@@ -4,16 +4,22 @@
 # Thu Oct 23 21:58:40 CEST 2014
 
 QPACKAGE_ROOT=${PWD}
+PACKAGES="core cryptokit"
 
+if [[ -f quantum_package.rc ]]
+then
+  source quantum_package.rc
+fi
 make -C ocaml Qptypes.ml &> /dev/null
 if [[ $? -ne 0 ]]
 then
 
+  rm -rf -- ${HOME}/ocamlbrew
   scripts/fetch_from_web.py "https://raw.github.com/hcarty/ocamlbrew/master/ocamlbrew-install" ocamlbrew-install.sh 
   cat < ocamlbrew-install.sh | env OCAMLBREW_FLAGS="-r" bash | tee ocamlbrew_install.log
   grep "source " ocamlbrew_install.log | grep "etc/ocamlbrew.bashrc"  >> quantum_package.rc
   source quantum_package.rc
-  echo Y | opam install core 
+  echo Y | opam install ${PACKAGES}
 fi
 make -C ocaml Qptypes.ml
 

scripts/upgrade_ezfio.sh

@@ -0,0 +1,26 @@
+#!/bin/bash
+# 
+# Upgrades the EZFIO library from the web.
+# Tue Nov  4 00:53:13 CET 2014
+
+if [[ -z ${QPACKAGE_ROOT} ]]
+then
+  print "The QPACKAGE_ROOT environment variable is not set."
+  print "Please reload the quantum_package.rc file."
+fi
+
+cd -- ${QPACKAGE_ROOT}
+mv -- ${QPACKAGE_ROOT}/EZFIO ${QPACKAGE_ROOT}/EZFIO.old
+
+make EZFIO
+
+if [[ $? -eq 0 ]]
+then
+  rm -rf -- ${QPACKAGE_ROOT}/EZFIO.old
+  echo "Successfully updated EZFIO"
+else
+  rm -rf -- ${QPACKAGE_ROOT}/EZFIO
+  mv --  ${QPACKAGE_ROOT}/EZFIO.old ${QPACKAGE_ROOT}/EZFIO
+  echo "Failed to update EZFIO"
+fi
+

scripts/upgrade_irpf90.sh

@@ -0,0 +1,25 @@
+#!/bin/bash
+# 
+# Upgrades IRPF90 from the web.
+# Tue Nov  4 00:53:13 CET 2014
+
+if [[ -z ${QPACKAGE_ROOT} ]]
+then
+  print "The QPACKAGE_ROOT environment variable is not set."
+  print "Please reload the quantum_package.rc file."
+fi
+
+cd -- ${QPACKAGE_ROOT}
+mv -- ${QPACKAGE_ROOT}/irpf90 ${QPACKAGE_ROOT}/irpf90.old
+
+make irpf90 
+
+if [[ $? -eq 0 ]]
+then
+  rm -rf -- ${QPACKAGE_ROOT}/irpf90.old
+  echo "Successfully updated IRPF90"
+else
+  rm -rf -- ${QPACKAGE_ROOT}/irpf90
+  mv --  ${QPACKAGE_ROOT}/irpf90.old ${QPACKAGE_ROOT}/irpf90
+  echo "Failed to update IRPF90"
+fi

setup_environment.sh

@@ -53,7 +53,7 @@ echo $RED "
 To complete the installation, add the following line to
 your ~/.bashrc:
 
-source quantum_package.rc
+source ${QPACKAGE_ROOT}/quantum_package.rc
 
 =======================================================
 " $BLACK

src/BiInts/map_integrals.irp.f

@@ -91,6 +91,17 @@ subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
       endif
     enddo
   enddo
+  do ii=1,8
+    if (i(ii) /= 0) then
+      call bielec_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
+      if (i1 /= i2) then
+        print *,  i1, i2
+        print *,  i(ii), j(jj), k(jj), l(jj)
+        stop 'bielec_integrals_index_reverse failed'
+      endif
+    endif
+  enddo
+
 
 end
 

src/Bitmask/bitmasks.ezfio_config

@@ -4,5 +4,5 @@ bitmasks
   N_mask_gen    integer
   generators    integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_gen)
   N_mask_ref    integer
-  reference     integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_ref)
+  reference     integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,2,bitmasks_N_mask_ref)
 

src/Bitmask/bitmasks.irp.f

@@ -175,10 +175,6 @@ BEGIN_PROVIDER [ integer(bit_kind), reference_bitmask, (N_int,2,2,N_reference_bi
  else
    reference_bitmask(:,:,s_hole ,1) = HF_bitmask
    reference_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
-   reference_bitmask(:,:,d_hole1,1) = HF_bitmask
-   reference_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
-   reference_bitmask(:,:,d_hole2,1) = HF_bitmask
-   reference_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
  endif
 
 END_PROVIDER

src/Dets/H_apply.irp.f

@@ -1,4 +1,5 @@
 use bitmasks
+use omp_lib
 
 type H_apply_buffer_type
 integer                        :: N_det
@@ -11,7 +12,8 @@ end type H_apply_buffer_type
 type(H_apply_buffer_type), pointer :: H_apply_buffer(:)
 
 
-BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
+ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
+&BEGIN_PROVIDER [ integer(omp_lock_kind), H_apply_buffer_lock, (64,0:nproc-1) ]
   use omp_lib
   implicit none
   BEGIN_DOC
@@ -24,7 +26,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
     allocate(H_apply_buffer(0:nproc-1))
     iproc = 0
     !$OMP PARALLEL PRIVATE(iproc) DEFAULT(NONE)  &
-    !$OMP SHARED(H_apply_buffer,N_int,sze,N_states)
+    !$OMP SHARED(H_apply_buffer,N_int,sze,N_states,H_apply_buffer_lock)
     !$   iproc = omp_get_thread_num()
     H_apply_buffer(iproc)%N_det = 0
     H_apply_buffer(iproc)%sze = sze
@@ -36,6 +38,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
     H_apply_buffer(iproc)%det  = 0_bit_kind
     H_apply_buffer(iproc)%coef = 0.d0
     H_apply_buffer(iproc)%e2   = 0.d0
+    call omp_init_lock(H_apply_buffer_lock(1,iproc))
     !$OMP END PARALLEL
   endif
   
@@ -56,6 +59,7 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   ASSERT (iproc >= 0)
   ASSERT (iproc < nproc)
   
+  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   allocate ( buffer_det(N_int,2,new_size),                           &
       buffer_coef(new_size,N_states),                                &
       buffer_e2(new_size,N_states) )
@@ -89,6 +93,7 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   
   H_apply_buffer(iproc)%sze = new_size
   H_apply_buffer(iproc)%N_det = min(new_size,H_apply_buffer(iproc)%N_det)
+  call omp_unset_lock(H_apply_buffer_lock(1,iproc))
   
 end
 
@@ -174,6 +179,7 @@ subroutine copy_H_apply_buffer_to_wf
   H_apply_buffer(j)%N_det = 0
   !$OMP END PARALLEL
   call normalize(psi_coef,N_det)
+  SOFT_TOUCH N_det psi_det psi_coef
   
 end
 
@@ -194,6 +200,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
   if (new_size > H_apply_buffer(iproc)%sze) then
     call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc)
   endif
+  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   do i=1,H_apply_buffer(iproc)%N_det
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
@@ -216,6 +223,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
   enddo
+  call omp_unset_lock(H_apply_buffer_lock(1,iproc))
 end
 
 

src/Dets/utils.irp.f

@@ -1,7 +1,7 @@
 BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
  implicit none
  BEGIN_DOC
- ! H matrix on the basis of the slater deter;inants defined by psi_det
+ ! H matrix on the basis of the slater determinants defined by psi_det
  END_DOC
  integer :: i,j
  double precision :: hij

src/MOGuess/H_CORE_guess.irp.f

@@ -3,7 +3,7 @@ program H_CORE_guess
   character*(64)                 :: label
   mo_coef = ao_ortho_lowdin_coef
   TOUCH mo_coef
-  label = "H_CORE_GUESS"
+  label = "Guess"
   call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),label)
   call save_mos
   

src/Makefile.common

@@ -1,3 +1,9 @@
+# Required EZFIO version
+EZFIO_VERSION=1.1
+
+# Required IRPF90 version
+IRPF90_VERSION=1.4
+
 # Check if QPACKAGE_ROOT is defined
 
 ifndef QPACKAGE_ROOT
@@ -28,10 +34,12 @@ include $(QPACKAGE_ROOT)/src/Makefile.config
 
 # Check the  version of IRPF90
 
-IRP_VERSION_OK=$(shell $(IRPF90) -v | python -c "import sys ; print float(sys.stdin.readline().rsplit('.',1)[0]) >= 1.3")
+IRP_VERSION_OK=$(shell $(IRPF90) -v | python -c "import sys ; print float(sys.stdin.readline().rsplit('.',1)[0]) >= $(IRPF90_VERSION)")
 ifeq ($(IRP_VERSION_OK),False)
 $(info  -------------------- Error --------------------)
-$(info  IRPF90 version >= 1.3 is required)
+$(info  IRPF90 version >= $(IRPF90_VERSION) is required)
+$(info  To upgrade IRPF90, run :                       )
+$(info    $(QPACKAGE_ROOT)/scripts/upgrade_irpf90.sh   )
 $(info  -----------------------------------------------)
 $(error )
 endif
@@ -85,12 +93,26 @@ $(info  -----------------------------------------------)
 endif
 
 
-# Define the Makefile common variables and rules
+# Define the Makefile common variables
 
 EZFIO_DIR=$(QPACKAGE_ROOT)/EZFIO
 EZFIO=$(EZFIO_DIR)/lib/libezfio_irp.a
 INCLUDE_DIRS=$(NEEDED_MODULES) include
 
+# Check EZFIO version
+
+EZFIO_VERSION_OK=$(shell cat $(EZFIO_DIR)/version | cut -d '=' -f 2 | python -c "import sys ; print float(sys.stdin.readline().rsplit('.',1)[0]) >= $(EZFIO_VERSION)")
+ifeq ($(EZFIO_VERSION_OK),False)
+$(info  -------------------- Error --------------------)
+$(info  EZFIO version >= $(EZFIO_VERSION) is required  )
+$(info  To upgrade EZFIO, run :                        )
+$(info    $(QPACKAGE_ROOT)/scripts/upgrade_ezfio.sh    )
+$(info  -----------------------------------------------)
+$(error )
+endif
+
+# Define the EZFIO rules
+
 $(EZFIO): $(wildcard $(QPACKAGE_ROOT)/src/*.ezfio_config) $(wildcard $(QPACKAGE_ROOT)/src/*/*.ezfio_config)
 	@echo Building EZFIO library
 	@cp $(wildcard $(QPACKAGE_ROOT)/src/*.ezfio_config) $(wildcard $(QPACKAGE_ROOT)/src/*/*.ezfio_config) $(EZFIO_DIR)/config 

src/Perturbation/selection.irp.f

@@ -26,6 +26,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
   if (new_size > h_apply_buffer(iproc)%sze) then
     call resize_h_apply_buffer(max(h_apply_buffer(iproc)%sze*2,new_size),iproc)
   endif
+  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   do i=1,H_apply_buffer(iproc)%N_det
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
@@ -63,6 +64,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
   enddo
+  call omp_unset_lock(H_apply_buffer_lock(1,iproc))
   !$OMP CRITICAL
   selection_criterion = max(selection_criterion,smax)
   selection_criterion_min = min(selection_criterion_min,smin)
@@ -146,6 +148,7 @@ subroutine make_s2_eigenfunction
   integer, parameter             :: bufsze = 1000
   logical, external              :: is_in_wavefunction
 
+  print *,  irp_here
 !  !TODO DEBUG
 !  do i=1,N_det
 !   do j=i+1,N_det
@@ -172,9 +175,11 @@ subroutine make_s2_eigenfunction
 
   do i=1,N_occ_pattern
     call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int)
+    s += 1
     if (s > smax) then
       deallocate(d)
       allocate ( d(N_int,2,s) )
+      smax = s
     endif
     call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int)
     do j=1,s