diff --git a/Basis/Basis.ml b/Basis/Basis.ml
index 7c1b869..d19fe2d 100644
--- a/Basis/Basis.ml
+++ b/Basis/Basis.ml
@@ -4,22 +4,23 @@ type t = Contracted_shell.t array
 let of_nuclei_and_general_basis n b =
   let result = 
     Array.map (fun (e, center) ->
-      List.assoc e b
-      |> Array.map (fun (totAngMom, shell) ->
-        let expo = Array.map (fun General_basis.{exponent ; coefficient} ->
-          exponent) shell
-        and coef = Array.map (fun General_basis.{exponent ; coefficient} ->
-          coefficient) shell
-        in
-        Contracted_shell.create ~expo ~coef ~totAngMom ~center ~indice:0)
-    ) n
+        List.assoc e b
+        |> Array.map (fun (totAngMom, shell) ->
+            let expo = Array.map (fun General_basis.{exponent ; coefficient} ->
+                exponent) shell
+            and coef = Array.map (fun General_basis.{exponent ; coefficient} ->
+                coefficient) shell
+            in
+            Contracted_shell.create ~expo ~coef ~totAngMom ~center ~index:0)
+      ) n
     |> Array.to_list
     |> Array.concat
   in
   Array.iteri (fun i x ->
-    if (i > 0) then
-      result.(i) <- {x with Contracted_shell.indice= (
-    result.(i-1).Contracted_shell.indice + (Array.length result.(i-1).Contracted_shell.powers)) }
+      if (i > 0) then
+        result.(i) <- Contracted_shell.with_index x (
+            (Contracted_shell.index result.(i-1)) +
+            (Array.length (Contracted_shell.powers result.(i-1)))) 
     ) result ;
   result
 
@@ -38,12 +39,12 @@ let to_string b =
 -----------------------------------------------------------------------
 "
   ^
-    ( Array.map (fun i -> 
-          Contracted_shell.to_string i) b
-      |> Array.to_list
-      |> String.concat line
-    )
-    ^ line
+  ( Array.map (fun i -> 
+        Contracted_shell.to_string i) b
+    |> Array.to_list
+    |> String.concat line
+  )
+  ^ line
 
 let file : string option ref = ref None
 
diff --git a/Basis/Basis.mli b/Basis/Basis.mli
index 9626e93..30be9c8 100644
--- a/Basis/Basis.mli
+++ b/Basis/Basis.mli
@@ -1,7 +1,30 @@
 type t = Contracted_shell.t array
+
+
+(** Returns an array of the basis set per atom *)
 val of_nuclei_and_general_basis : Nuclei.t -> General_basis.t list -> t
+
+
+(** Pretty prints the basis set in a string *)
 val to_string : t -> string
+
+
+(** Mutates the state of the file variable to Some f. Required for command-line
+    interface.
+*)
 val set_file : string -> unit
+
+
+(** Basis set file read and parsed. Requires the set_file function to have
+    been called before.
+*)
 val general_basis :
   (Element.t * General_basis.general_contracted_shell array) list lazy_t
+
+
+(** Global variable which sets the basis set of the current run.
+    Lazy evaluated from the nuclear coordinates (Nuclei.nuclei) and
+    the basis set file (general_basis). The set_file function has
+    to be called before the basis is used.
+*)
 val basis : Contracted_shell.t array lazy_t
diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml
index 78442e1..2b294ff 100644
--- a/Basis/Contracted_shell.ml
+++ b/Basis/Contracted_shell.ml
@@ -9,7 +9,7 @@ type t = {
   size      : int;
   norm_coef : float array;
   norm_coef_scale : float array;
-  indice    : int;
+  index    : int;
   powers   : Zkey.t array;
 }
 
@@ -20,7 +20,8 @@ let center a = a.center
 let totAngMom a = a.totAngMom
 let norm_coef a i = a.norm_coef.(i)
 let norm_coef_scale a  = a.norm_coef_scale
-let indice a  = a.indice
+let index a  = a.index
+let with_index a i = { a with index = i }
 let powers a = a.powers
 
 
@@ -30,8 +31,8 @@ let to_string s =
   in
   let open Printf in
   (match s.totAngMom with
-   | Angular_momentum.S -> sprintf "%3d    " (s.indice+1)
-   | _                  -> sprintf "%3d-%-3d" (s.indice+1) (s.indice+(Array.length s.powers))
+   | Angular_momentum.S -> sprintf "%3d    " (s.index+1)
+   | _                  -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.powers))
   ) ^
   ( sprintf "%1s %8.3f %8.3f %8.3f " (Angular_momentum.to_string s.totAngMom)
       (Coordinate.x coord) (Coordinate.y coord) (Coordinate.z coord) ) ^
@@ -71,7 +72,7 @@ let compute_norm_coef expo totAngMom =
   Array.map (fun x -> let f a = x *. (factor a) in f) expo
 
 
-let create ~indice ~expo ~coef ~center ~totAngMom =
+let create ~index ~expo ~coef ~center ~totAngMom =
   assert (Array.length expo = Array.length coef);
   assert (Array.length expo > 0);
   let norm_coef_func =
@@ -88,7 +89,7 @@ let create ~indice ~expo ~coef ~center ~totAngMom =
       (norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0)
     ) powers
   in
-  { indice ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
+  { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
     norm_coef_scale ; powers }
 
 
diff --git a/Basis/Contracted_shell.mli b/Basis/Contracted_shell.mli
new file mode 100644
index 0000000..e3f6b2b
--- /dev/null
+++ b/Basis/Contracted_shell.mli
@@ -0,0 +1,67 @@
+type t
+(*
+type t = {
+  expo : float array;
+  coef : float array;
+  center : Coordinate.t;
+  totAngMom : Angular_momentum.t;
+  size : int;
+  norm_coef : float array;
+  norm_coef_scale : float array;
+  indice : int;
+  powers : Zkey.t array;
+}
+*)
+
+(** Returns the number of contracted Gaussians *)
+val size : t -> int
+
+
+(** Returns the i-th exponent *)
+val expo : t -> int -> float
+
+
+(** Returns the i-th contraction coefficient *)
+val coef : t -> int -> float
+
+
+(** Point on which all the Gaussians are centered *)
+val center : t -> Coordinate.t
+
+
+(** Total angular momentum *)
+val totAngMom : t -> Angular_momentum.t
+
+
+(** Normalization coefficient of the class corresponding to the i-th contraction *)
+val norm_coef : t -> int -> float
+
+
+(** Inside a class, the norm is the norm of the function with (totAngMom,0,0) *. 
+    this scaling factor *)
+val norm_coef_scale : t -> float array
+
+
+(** The index in the array of contracted shells *)
+val index  : t -> int
+
+
+(** Returns a copy of the contracted shell with a modified index *)
+val with_index  : t -> int -> t
+
+
+(** The array of Zkeys corresponding to the powers of (x,y,z) in the class *)
+val powers : t -> Zkey.t array
+
+
+(** Pretty-printing of the contracted shell in a string *)
+val to_string : t -> string
+
+(** Creates a contracted shell *)
+val create :
+  index:int ->
+  expo:float array ->
+  coef:float array ->
+  center:Coordinate.t -> totAngMom:Angular_momentum.t -> t
+
+
diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index 2897046..1f4772c 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -94,7 +94,7 @@ let to_file ~filename basis =
 
   let icount = ref 0 in
   for i=0 to (Array.length basis) - 1 do
-    print_int basis.(i).Contracted_shell.indice ; print_newline ();
+    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
     for j=0 to i do
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
         if (schwartz_p_max >= cutoff) then
@@ -107,7 +107,7 @@ let to_file ~filename basis =
 
     let n = ref 0 in
     for i=0 to (Array.length basis) - 1 do
-      n := !n + (Array.length (basis.(i).Contracted_shell.powers))
+      n := !n + (Array.length (Contracted_shell.powers (basis.(i))))
     done;
     let n = !n in
     Genarray.create Float64 c_layout [| n ; n ; n ; n|]
@@ -120,7 +120,7 @@ let to_file ~filename basis =
   let inn = ref 0 and out = ref 0 in
 
   for i=0 to (Array.length basis) - 1 do
-    print_int basis.(i).Contracted_shell.indice ; print_newline ();
+    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
     for j=0 to i do
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
       try
@@ -159,16 +159,16 @@ let to_file ~filename basis =
 
             (* Write the data in the output file *)
             Array.iteri (fun i_c powers_i ->
-                let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
+                let i_c = (Contracted_shell.index basis.(i)) + i_c + 1 in
                 let xi = to_int_tuple powers_i in
                 Array.iteri (fun j_c powers_j ->
-                    let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
+                    let j_c = (Contracted_shell.index basis.(j)) + j_c + 1 in
                     let xj = to_int_tuple powers_j in
                     Array.iteri (fun k_c powers_k ->
-                        let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in 
+                        let k_c = (Contracted_shell.index basis.(k)) + k_c + 1 in 
                         let xk = to_int_tuple powers_k in
                         Array.iteri (fun l_c powers_l ->
-                            let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in
+                            let l_c = (Contracted_shell.index basis.(l)) + l_c + 1 in
                             let xl = to_int_tuple powers_l in
                             let key = 
                               if swap then
@@ -192,10 +192,10 @@ let to_file ~filename basis =
                               )
                             else
                               out := !out + 1;
-                          ) basis.(l).Contracted_shell.powers 
-                      ) basis.(k).Contracted_shell.powers
-                  ) basis.(j).Contracted_shell.powers
-              ) basis.(i).Contracted_shell.powers;
+                          ) (Contracted_shell.powers basis.(l)) 
+                      ) (Contracted_shell.powers basis.(k))
+                  ) (Contracted_shell.powers basis.(j))
+              ) (Contracted_shell.powers basis.(i));
               with NullIntegral ->  ()
           done;
         done;
diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
index 18dff5c..31acd68 100644
--- a/Basis/KinInt.ml
+++ b/Basis/KinInt.ml
@@ -108,7 +108,7 @@ let to_file ~filename basis =
 
   let oc = open_out filename in
   for i=0 to (Array.length basis) - 1 do
-    print_int basis.(i).Contracted_shell.indice ; print_newline ();
+    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
     for j=0 to i do
       (* Compute all the integrals of the class *)
       let cls =
@@ -117,10 +117,10 @@ let to_file ~filename basis =
 
       (* Write the data in the output file *)
       Array.iteri (fun i_c powers_i ->
-          let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
+          let i_c = Contracted_shell.index basis.(i) + i_c + 1 in
           let xi = to_int_tuple powers_i in
           Array.iteri (fun j_c powers_j ->
-              let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
+              let j_c = Contracted_shell.index basis.(j) + j_c + 1 in
               let xj = to_int_tuple powers_j in
               let key = 
                 Zkey.of_int_tuple (Zkey.Six (xi,xj))
@@ -133,8 +133,8 @@ let to_file ~filename basis =
               if (abs_float value > cutoff) then
                 Printf.fprintf oc "%4d %4d %20.12e\n"
                   i_c j_c value
-            ) basis.(j).Contracted_shell.powers
-        ) basis.(i).Contracted_shell.powers;
+            ) (Contracted_shell.powers basis.(j))
+        ) (Contracted_shell.powers basis.(i));
     done;
   done;
   close_out oc
diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index 12c1a59..0e4bda4 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -52,7 +52,7 @@ let to_file ~filename basis geometry =
   let eni_array = 
     let n = ref 0 in
     for i=0 to (Array.length basis) - 1 do
-      n := !n + (Array.length (basis.(i).Contracted_shell.powers))
+      n := !n + (Array.length (Contracted_shell.powers basis.(i)))
     done;
     let n = !n in
     Array2.create Float64 c_layout n  n
@@ -65,7 +65,7 @@ let to_file ~filename basis geometry =
   let inn = ref 0 and out = ref 0 in
 
   for i=0 to (Array.length basis) - 1 do
-    print_int basis.(i).Contracted_shell.indice ; print_newline ();
+    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
     for j=0 to i do
       let
         shell_p = shell_pairs.(i).(j)
@@ -78,10 +78,10 @@ let to_file ~filename basis geometry =
 
       (* Write the data in the output file *)
       Array.iteri (fun i_c powers_i ->
-          let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
+          let i_c = (Contracted_shell.index basis.(i)) + i_c + 1 in
           let xi = to_int_tuple powers_i in
           Array.iteri (fun j_c powers_j ->
-              let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
+              let j_c = (Contracted_shell.index basis.(j)) + j_c + 1 in
               let xj = to_int_tuple powers_j in
               let key = 
                 Zkey.of_int_tuple (Zkey.Six (xi,xj))
@@ -95,8 +95,8 @@ let to_file ~filename basis geometry =
                 )
               else
                 out := !out + 1;
-            ) basis.(j).Contracted_shell.powers
-        ) basis.(i).Contracted_shell.powers;
+            ) (Contracted_shell.powers basis.(j))
+        ) (Contracted_shell.powers basis.(i));
     done;
   done;
 
diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml
index bac802c..ea4a0b5 100644
--- a/Basis/OneElectronRR.ml
+++ b/Basis/OneElectronRR.ml
@@ -168,7 +168,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
           shell_p.(ab).Shell_pair.center
         in
         let center_pa = 
-          Coordinate.(center_p |- shell_a.Contracted_shell.center) 
+          Coordinate.(center_p |- Contracted_shell.center shell_a) 
         in
 
         for c=0 to Array.length geometry - 1 do
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index 4c31890..f75d013 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -80,7 +80,7 @@ let to_file ~filename basis =
 
   let oc = open_out filename in
   for i=0 to (Array.length basis) - 1 do
-    print_int basis.(i).Contracted_shell.indice ; print_newline ();
+    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
     for j=0 to i do
       (* Compute all the integrals of the class *)
       let cls =
@@ -89,10 +89,10 @@ let to_file ~filename basis =
 
       (* Write the data in the output file *)
       Array.iteri (fun i_c powers_i ->
-          let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
+          let i_c = Contracted_shell.index basis.(i) + i_c + 1 in
           let xi = to_int_tuple powers_i in
           Array.iteri (fun j_c powers_j ->
-              let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
+              let j_c = Contracted_shell.index basis.(j) + j_c + 1 in
               let xj = to_int_tuple powers_j in
               let key = 
                 Zkey.of_int_tuple (Zkey.Six (xi,xj))
@@ -105,8 +105,8 @@ let to_file ~filename basis =
               if (abs_float value > cutoff) then
                 Printf.fprintf oc "%4d %4d %20.12e\n"
                   i_c j_c value
-            ) basis.(j).Contracted_shell.powers
-        ) basis.(i).Contracted_shell.powers;
+            ) (Contracted_shell.powers basis.(j))
+        ) (Contracted_shell.powers basis.(i));
     done;
   done;
   close_out oc
diff --git a/README.rst b/README.rst
index 09f23ea..233261b 100644
--- a/README.rst
+++ b/README.rst
@@ -1,7 +1,6 @@
 Requirements
 ------------
 
-* gsl: GNU Scientific Library
 * gmp : GNU Multiple Precision arithmetic library
 * zarith : Arbitrary-precision integers
 * BLAS/LAPACK : Linear algebra