diff --git a/.ocamlinit b/.ocamlinit
new file mode 100644
index 0000000..55aa028
--- /dev/null
+++ b/.ocamlinit
@@ -0,0 +1,28 @@
+let () =
+  try Topdirs.dir_directory (Sys.getenv "OCAML_TOPLEVEL_PATH")
+  with Not_found -> ()
+;;
+
+#use "topfind";;
+#require "lacaml";;
+#directory "_build";;
+#directory "_build/Utils";;
+#directory "_build/Basis";;
+#directory "_build/HartreeFock";;
+
+#load "Powers.cmo";;
+#load "Zkey.cmo";;
+#load "AngularMomentum.cmo";;
+#load "Constants.cmo";;
+#load "Util.cmo";;
+#load "Coordinate.cmo";;
+#load "Zmap.cmo";;
+#load "Zkey.cmo";;
+#load "ContractedShell.cmo";;
+#load "ShellPair.cmo";;
+#load "TwoElectronRR.cmo";;
+#load "ERI.cmo";;
+
+#use "topfind";;
+
+
diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index b707cfa..fa3385c 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -10,6 +10,8 @@ module Bs  = Basis
 module Cs  = ContractedShell
 module Csp = ContractedShellPair
 
+let cutoff = integrals_cutoff
+
 (** (00|00)^m : Fundamental electron repulsion integral
     $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $
 
diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
index e82b3a0..3e6b02e 100644
--- a/Basis/KinInt.ml
+++ b/Basis/KinInt.ml
@@ -11,6 +11,7 @@ module Sp  = ShellPair
 
 type t = Mat.t
 
+let cutoff = integrals_cutoff
 
 (** Computes all the kinetic integrals of the contracted shell pair *)
 let contracted_class shell_a shell_b : float Zmap.t = 
diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index feacf29..8428032 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -33,6 +33,8 @@ let contracted_class_shell_pair shell_p geometry: float Zmap.t =
   OneElectronRR.contracted_class_shell_pair ~zero_m shell_p geometry
 
 
+let cutoff = integrals_cutoff
+
 let cutoff2 = cutoff *. cutoff
 exception NullIntegral
 
diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml
index 4b76d03..e6b1345 100644
--- a/Basis/OneElectronRR.ml
+++ b/Basis/OneElectronRR.ml
@@ -10,11 +10,6 @@ module Csp = ContractedShellPair
 module Po  = Powers
 module Sp  = ShellPair
 
-(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *)
-let chop f g =
-  if (abs_float f) < cutoff then 0.
-  else f *. (g ())
-
 
 
 (** Horizontal and Vertical Recurrence Relations (HVRR) *)
@@ -84,7 +79,7 @@ let hvrr_one_e  angMom_a  angMom_b
           let bm = Po.decr xyz angMom_b in
           let h1 = hrr ap bm in
           let f2 = Co.get xyz center_ab in
-          if abs_float f2 < cutoff then h1 else
+          if abs_float f2 < integrals_cutoff then h1 else
               let h2 = hrr angMom_a bm in
               h1 +. f2 *. h2
 
@@ -130,7 +125,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
         let coef_prod = shell_p.Csp.coef.(ab) in
 
         (** Screening on the product of coefficients *)
-        if abs_float coef_prod < 1.e-3 *. cutoff then
+        if abs_float coef_prod < 1.e-3 *. integrals_cutoff then
           raise NullPair;
 
 
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index e2c6d71..c5985e2 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -11,6 +11,8 @@ module Cs  = ContractedShell
 module Csp = ContractedShellPair
 module Sp  = ShellPair
 
+let cutoff = integrals_cutoff
+
 (** Computes all the overlap integrals of the contracted shell pair *)
 let contracted_class shell_a shell_b : float Zmap.t = 
 
@@ -45,7 +47,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
       shell_p.Csp.coef.(ab) 
     in
     (** Screening on thr product of coefficients *)
-    if (abs_float coef_prod) > 1.e-4*.cutoff then
+    if (abs_float coef_prod) > 1.e-3*.cutoff then
       begin
         let expo_inv = 
           shell_p.Csp.expo_inv.(ab)
diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml
index 9684f04..99c48bd 100644
--- a/Basis/TwoElectronRR.ml
+++ b/Basis/TwoElectronRR.ml
@@ -1,5 +1,4 @@
 open Util
-open Constants
 
 module Am  = AngularMomentum
 module Co  = Coordinate
@@ -8,8 +7,7 @@ module Csp = ContractedShellPair
 module Sp  = ShellPair
 module Po  = Powers
 
-let debug=false
-
+let cutoff = Constants.integrals_cutoff 
 let cutoff2 = cutoff *. cutoff
 
 exception NullQuartet
diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 903403e..3e40d37 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -11,7 +11,7 @@ module Po  = Powers
 exception NullQuartet
 exception Found
 
-let cutoff = Constants.cutoff 
+let cutoff = Constants.integrals_cutoff 
 let cutoff2 = cutoff *. cutoff
 
 let at_least_one_valid arr =
diff --git a/Makefile b/Makefile
index b7b083a..7449bba 100644
--- a/Makefile
+++ b/Makefile
@@ -11,7 +11,9 @@ MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) Utils/math_functions.c
 MLYFILES=$(wildcard */*.mly) $(wildcard *.mly) 
 MLFILES= $(wildcard */*.ml)  $(wildcard *.ml) 
 MLIFILES=$(wildcard */*.mli) $(wildcard *.mli) 
-ALL_EXE=$(patsubst %.ml,%.native,$(wildcard run_*.ml)) 
+ALL_NATIVE=$(patsubst %.ml,%.native,$(wildcard run_*.ml)) 
+ALL_BYTE=$(patsubst %.ml,%.byte,$(wildcard run_*.ml)) 
+ALL_EXE=$(ALL_BYTE) $(ALL_NATIVE)
 
 .PHONY: default 
 
@@ -32,7 +34,6 @@ doc: qpackage.odocl
 %.byte: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES)
 	rm -f -- $*
 	$(OCAMLBUILD) $*.byte  -use-ocamlfind  $(PKGS) 
-	ln -s $*.byte $*
 
 %.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) 
 	rm -f -- $*
@@ -42,12 +43,10 @@ doc: qpackage.odocl
 %.p.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) 
 	rm -f -- $*
 	$(OCAMLBUILD) $*.p.native -use-ocamlfind $(PKGS)
-	ln -s $*.p.native $*
 
 %.p.byte: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) 
 	rm -f -- $*
 	$(OCAMLBUILD) -ocamlc ocamlcp $*.byte -use-ocamlfind $(PKGS)
-	ln -s $*.byte $*
 
 clean: 
 	rm -rf _build $(ALL_EXE) $(ALL_TESTS) *.native *.byte
diff --git a/Nuclei/Element.ml b/Nuclei/Element.ml
index 159671b..50275e1 100644
--- a/Nuclei/Element.ml
+++ b/Nuclei/Element.ml
@@ -141,7 +141,7 @@ let covalent_radius x =
   |  Te  ->  1.38  |  I   ->  1.39  |  Xe  ->  1.40  |  Pt  ->  1.30
   in
   Constants.a0 *. (result x)
-  |> Radius.of_float
+  |> NonNegativeFloat.of_float
 
 
 let vdw_radius x = 
@@ -162,7 +162,7 @@ let vdw_radius x =
   |  Te  ->  2.06  |  I   ->  1.98  |  Xe  ->  2.16  |  Pt  ->  1.75
   in
   Constants.a0 *. (result x)
-  |> Radius.of_float
+  |> NonNegativeFloat.of_float
 
 
 let mass c =
diff --git a/Nuclei/Element.mli b/Nuclei/Element.mli
index 3c78084..2f8cba7 100644
--- a/Nuclei/Element.mli
+++ b/Nuclei/Element.mli
@@ -20,6 +20,6 @@ val to_int : t -> int
 val of_int : int -> t
 val to_charge : t -> Charge.t
 val of_charge : Charge.t -> t
-val covalent_radius : t -> Radius.t
-val vdw_radius : t -> Radius.t
+val covalent_radius : t -> NonNegativeFloat.t
+val vdw_radius : t -> NonNegativeFloat.t
 val mass : t -> Mass.t
diff --git a/Nuclei/Mass.ml b/Nuclei/Mass.ml
index 035c570..98c4393 100644
--- a/Nuclei/Mass.ml
+++ b/Nuclei/Mass.ml
@@ -1 +1 @@
-include PositiveFloat
+include NonNegativeFloat
diff --git a/Utils/Angstrom.ml b/Utils/Angstrom.ml
index b553427..0d3883b 100644
--- a/Utils/Angstrom.ml
+++ b/Utils/Angstrom.ml
@@ -1,12 +1 @@
-type t = {
-  x : float ;
-  y : float ;
-  z : float ;
-}
-
-
-let make { Point.x ; y ; z } = { x ; y ; z } 
-(* = %identity *)
-  
-
-
+include Bohr
diff --git a/Utils/Angstrom.mli b/Utils/Angstrom.mli
index 77aa9e2..fa6fd98 100644
--- a/Utils/Angstrom.mli
+++ b/Utils/Angstrom.mli
@@ -1,3 +1,5 @@
+(** 3D point in cartesian coordinates, where units are Angstroms. *)
+
 type t = private {
   x : float ;
   y : float ;
@@ -7,4 +9,3 @@ type t = private {
 
 val make : Point.t -> t
 
-
diff --git a/Utils/AngularMomentum.mli b/Utils/AngularMomentum.mli
new file mode 100644
index 0000000..dd0bea9
--- /dev/null
+++ b/Utils/AngularMomentum.mli
@@ -0,0 +1,87 @@
+(** Azimuthal quantum number, represented as s,p,d,... *)
+
+type t = S | P | D | F | G | H | I | J | K | L | M | N | O
+
+exception AngularMomentumError of string
+(** Raised when the {!AngularMomentum.t} element can't be created.
+  *)
+
+val of_char : char -> t
+(** Returns an {!AngularMomentum.t} when a shell is given as a character (case
+  insensitive).
+
+  Example:
+
+  [AngularMomentum.of_char 'p' -> AngularMomentum.P]
+ *)
+
+val to_string : t -> string
+(** 
+  [AngularMomentum.(to_string D) -> "D"]
+ *)
+
+
+val to_char : t -> char
+(** 
+  [AngularMomentum.(to_char D) -> 'D']
+ *)
+
+
+val to_int : t -> int
+(** 
+  Returns the l{_max} value of the shell.
+
+  Example:
+
+  [AngularMomentum.of_int 3 -> AngularMomentum.F]
+ *)
+
+val of_int : int -> t
+(** 
+  Opposite of {!of_int}.
+
+  Example:
+
+  [AngularMomentum.of_int 3 -> AngularMomentum.F]
+ *)
+
+type kind =
+    Singlet of t
+  | Doublet of (t * t)
+  | Triplet of (t * t * t)
+  | Quartet of (t * t * t * t)
+
+
+val n_functions : t -> int
+(** Number of cartesian functions in shell.
+
+  Example:
+
+  [AngularMomentum.n_functions D -> 6]
+  *)
+
+
+val zkey_array : kind -> Zkey.t array
+(** Array of {!Zkey.t}, where each element is a a key associated with the
+  the powers of x,y,z.
+
+  Example:
+
+  {[
+  AngularMomentum.( zkey_array Doublet (P,S) ) ->
+    [| {Zkey.left = 0; right = 1125899906842624} ;
+       {Zkey.left = 0; right = 1099511627776}    ;
+       {Zkey.left = 0; right = 1073741824}       |] 
+  =
+ 
+  let s,x,y,z =
+      Powers.( of_int_tuple (0,0,0),
+               of_int_tuple (1,0,0),
+               of_int_tuple (0,1,0), 
+               of_int_tuple (0,0,1) ) 
+  in
+  Array.map (fun (a,b) -> {!Zkey.of_powers_six} a b)
+      [| (x,s) ; (y,s) ; (z,s) |]
+  ]}
+
+  *)
diff --git a/Utils/Bohr.ml b/Utils/Bohr.ml
index b553427..c51cda6 100644
--- a/Utils/Bohr.ml
+++ b/Utils/Bohr.ml
@@ -5,8 +5,7 @@ type t = {
 }
 
 
-let make { Point.x ; y ; z } = { x ; y ; z } 
-(* = %identity *)
+external make : Point.t -> t = "%identity"
   
 
 
diff --git a/Utils/Bohr.mli b/Utils/Bohr.mli
index 77aa9e2..dd46b40 100644
--- a/Utils/Bohr.mli
+++ b/Utils/Bohr.mli
@@ -1,3 +1,5 @@
+(** 3D point in cartesian coordinates, where units are atomic units (Bohr). *)
+
 type t = private {
   x : float ;
   y : float ;
diff --git a/Utils/Constants.ml b/Utils/Constants.ml
index 9d8f9f9..d2b1b25 100644
--- a/Utils/Constants.ml
+++ b/Utils/Constants.ml
@@ -1,4 +1,5 @@
-let cutoff = 1.e-15
+let integrals_cutoff  = 1.e-15
+let epsilon = 1.e-20
 
 (** Constants *)
 let pi = acos (-1.)
@@ -7,5 +8,5 @@ let sq_pi_over_two = sq_pi *. 0.5
 let pi_inv = 1. /. pi
 let two_over_sq_pi = 2. /. (sqrt pi)
 
-let a0 = 0.529_177_210_67
+let a0 = 0.529_177_210_671_2
 let a0_inv = 1. /. a0
diff --git a/Utils/Constants.mli b/Utils/Constants.mli
index 48541b6..73a1b08 100644
--- a/Utils/Constants.mli
+++ b/Utils/Constants.mli
@@ -1,8 +1,26 @@
-val cutoff : float
+val epsilon : float
+(** Value below which a float is considered null. Default is 10{^-20}. *)
+
+val integrals_cutoff : float
+(** Cutoff value for integrals. Default is 10{^-15}. *)
+
 val pi : float
+(** pi = 3.141_592_653_589_793_12 *)
+
 val sq_pi : float
+(** [sqrt pi] *)
+
 val sq_pi_over_two : float
+(** [(sqrt pi) /. 2.] *)
+
 val pi_inv : float
+(** [ 1. /. pi ] *)
+
 val two_over_sq_pi : float
+(** [ 2. /. (sqrt pi) ] *)
+
 val a0 : float
+(** Bohr radius : a{_0} = 0.529_177_210_671_2 Angstrom *)
+
 val a0_inv : float
+(** [ 1. /. a0 ] *)
diff --git a/Utils/Coordinate_type.ml b/Utils/Coordinate_type.ml
index b2af9a8..25f8504 100644
--- a/Utils/Coordinate_type.ml
+++ b/Utils/Coordinate_type.ml
@@ -1,12 +1,6 @@
-type point = Point.t
-
-type bohr = Bohr.t
-type angstrom = Angstrom.t
-
-let a_to_b a = Constants.a0_inv *. a
-let b_to_a b = Constants.a0     *. b
 
 let bohr_to_angstrom { Bohr.x ; y ; z } = 
+  let b_to_a b = Constants.a0 *. b in
   Angstrom.make 
   {
     Point.
@@ -15,7 +9,10 @@ let bohr_to_angstrom { Bohr.x ; y ; z } =
     z = b_to_a z ;
   }
 
+
+
 let angstrom_to_bohr { Angstrom.x ; y ; z } =  
+  let a_to_b a = Constants.a0_inv *. a in
   Bohr.make 
   {
     Point.
@@ -24,3 +21,4 @@ let angstrom_to_bohr { Angstrom.x ; y ; z } =
     z = a_to_b z ;
   }
 
+
diff --git a/Utils/Coordinate_type.mli b/Utils/Coordinate_type.mli
index e3f1c00..5d2cbb2 100644
--- a/Utils/Coordinate_type.mli
+++ b/Utils/Coordinate_type.mli
@@ -1,3 +1,6 @@
+(** Coordinate_type can be Bohr or Angstrom. This module provides the conversion
+  functions from {!Bohr.t} to {!Angstrom.t} and from {!Angstrom.t} to {!Bohr.t} *)
+
 type bohr = Bohr.t
 type angstrom = Angstrom.t
 
diff --git a/Utils/Electrons.mli b/Utils/Electrons.mli
new file mode 100644
index 0000000..c5f750c
--- /dev/null
+++ b/Utils/Electrons.mli
@@ -0,0 +1,17 @@
+(** Information related to electrons. *)
+
+type t = {
+  n_alpha : int ;  (** Number of alpha electrons *)
+  n_beta  : int ;  (** Number of beta  electrons *)
+  multiplicity : int; (** Spin multiplicity: 2S+1 *)
+}
+
+
+val make : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t
+(** Creates the data relative to electrons for a molecular system
+  described by {!Nuclei.t} for a given total charge and spin multiplicity.
+  @param multiplicity default is 1
+  @param charge default is 0
+  @raise Invalid_argument if the spin multiplicity is not compatible with
+    the molecule and the total charge.
+  *)
diff --git a/Utils/PositiveFloat.ml b/Utils/NonNegativeFloat.ml
similarity index 81%
rename from Utils/PositiveFloat.ml
rename to Utils/NonNegativeFloat.ml
index fbf46c0..fe2bc8b 100644
--- a/Utils/PositiveFloat.ml
+++ b/Utils/NonNegativeFloat.ml
@@ -5,6 +5,7 @@ let of_float x =
   x
 
 external to_float : t -> float = "%identity"
+external unsafe_of_float : float -> t = "%identity"
 
 let to_string x =
   let f = to_float x in string_of_float f
diff --git a/Utils/PositiveFloat.mli b/Utils/NonNegativeFloat.mli
similarity index 100%
rename from Utils/PositiveFloat.mli
rename to Utils/NonNegativeFloat.mli
diff --git a/Utils/Radius.ml b/Utils/Radius.ml
deleted file mode 100644
index 035c570..0000000
--- a/Utils/Radius.ml
+++ /dev/null
@@ -1 +0,0 @@
-include PositiveFloat
diff --git a/Utils/Util.ml b/Utils/Util.ml
index 5feb872..68819b2 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -1,13 +1,22 @@
-(** Functions from libm *)
-external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed]  [@@noalloc]
-external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed]  [@@noalloc]
-external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed]  [@@noalloc]
+(** All utilities which should be included in all source files are defined here *)
+
+(** {1 Functions from libm} *)
 
 open Constants
 open Lacaml.D
 
 let factmax = 150
 
+external erf_float : float -> float = "erf_float_bytecode" "erf_float"
+  [@@unboxed] [@@noalloc]
+
+external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float"
+  [@@unboxed] [@@noalloc]
+
+external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float"
+  [@@unboxed] [@@noalloc]
+
+
 
 
 
@@ -79,8 +88,6 @@ let fact_memo =
 
 
 
-(** Factorial function. 
-    @raise Invalid_argument for negative or arguments >100. *)
 let fact = function
   | i when (i < 0) ->
     raise (Invalid_argument "Argument of factorial should be non-negative")
@@ -89,24 +96,23 @@ let fact = function
   | i -> fact_memo.(i)
 
 
-
-(** Integer powers of floats *)
 let rec pow a = function
   | 0 -> 1.
   | 1 -> a
   | 2 -> a *. a
   | 3 -> a *. a *. a
   | -1 -> 1. /. a
-  | n when (n<0) -> pow (1./.a) (-n)
-  | n ->
+  | n when  n > 0  ->
     let b = pow a (n / 2) in
     b *. b *. (if n mod 2 = 0 then 1. else a)
-;;
+  | n when  n < 0  -> pow (1./.a) (-n)
+  | _ -> assert false
+
+
 
 
-(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *)
 let chop f g =
-  if (abs_float f) < cutoff then 0.
+  if (abs_float f) < Constants.epsilon then 0.
   else f *. (g ())
 
 
@@ -162,7 +168,7 @@ let diagonalize_symm h =
   in
   v, result
 
-let xt_o_x o x =
+let xt_o_x ~o ~x =
   gemm o x
   |> gemm ~transa:`T x
 
diff --git a/Utils/Util.mli b/Utils/Util.mli
new file mode 100644
index 0000000..4560f6c
--- /dev/null
+++ b/Utils/Util.mli
@@ -0,0 +1,67 @@
+(** All utilities which should be included in all source files are defined here *)
+
+(** {1 Functions from libm} *)
+
+external erf_float : float -> float = "erf_float_bytecode" "erf_float"
+  [@@unboxed] [@@noalloc]
+  (** Error function [erf] from [libm] *)
+
+external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float"
+  [@@unboxed] [@@noalloc]
+  (** Complementary error function [erfc] from [libm] *)
+
+external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float"
+  [@@unboxed] [@@noalloc]
+  (** Gamma function [gamma] from [libm] *)
+
+
+
+(** {1 General functions} *)
+
+val fact : int -> float
+(** Factorial function.
+  @raise Invalid_argument for negative arguments or arguments >100.
+  *)
+
+val pow : float -> int -> float
+(** Fast implementation of the power function for small integer powers *)
+
+val chop : float -> (unit -> float) -> float
+(** In [chop a f], evaluate [f] only if the absolute value of [a] is larger
+  than {!Constants.epsilon}, and return [a *. f ()].
+  *)
+
+
+
+(** {1 Functions related to the Boys function} *)
+
+val incomplete_gamma : alpha:float -> float -> float
+(** {{:https://en.wikipedia.org/wiki/Incomplete_gamma_function}
+    Lower incomplete gamma function}
+   @raise Failure when the calculation doesn't converge.
+  *)
+
+val boys_function : maxm:int -> float -> float array
+(** {{:https://link.springer.com/article/10.1007/s10910-005-9023-3}
+    Generalized Boys function}.
+  @param maxm Maximum total angular momentum.
+  *)
+
+
+(** {1 Extension of the Array module} *)
+
+val array_sum : float array -> float
+(** Returns the sum of all the elements of the array *)
+
+val array_product : float array -> float
+(** Returns the product of all the elements of the array *)
+
+
+(** {1 Linear algebra } *)
+
+val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec
+(** Diagonalize a symmetric matrix. Returns the eigenvectors and the eigenvalues. *)
+
+val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
+(** Computes X{^T}.O.X *)
+