diff --git a/Basis/Basis.ml b/Basis/Basis.ml
index 6ff464b..5259dc7 100644
--- a/Basis/Basis.ml
+++ b/Basis/Basis.ml
@@ -4,8 +4,8 @@ type t =
     contracted_shells : Contracted_shell.t array;
   }
 
-let size b = b.size
-let contracted_shells b = b.contracted_shells
+module Cs = Contracted_shell
+module Gb = General_basis
 
 (** Returns an array of the basis set per atom *)
 let of_nuclei_and_general_basis n b =
@@ -13,27 +13,27 @@ let of_nuclei_and_general_basis n b =
     Array.map (fun (e, center) ->
         List.assoc e b
         |> Array.map (fun (totAngMom, shell) ->
-            let expo = Array.map (fun General_basis.{exponent ; coefficient} ->
+            let expo = Array.map (fun Gb.{exponent ; coefficient} ->
                 exponent) shell
-            and coef = Array.map (fun General_basis.{exponent ; coefficient} ->
+            and coef = Array.map (fun Gb.{exponent ; coefficient} ->
                 coefficient) shell
             in
-            Contracted_shell.make ~expo ~coef ~totAngMom ~center ~index:0)
+            Cs.make ~expo ~coef ~totAngMom ~center ~index:0)
       ) n
     |> Array.to_list
     |> Array.concat
   in
   Array.iteri (fun i x ->
       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.(i) <- Cs.with_index x (
+            result.(i-1).index  +
+            (Array.length result.(i-1).powers )) 
     ) result ;
 
   let size =
     let n = ref 0 in
     for i=0 to (Array.length result) - 1 do
-      n := !n + (Array.length (Contracted_shell.powers (result.(i))))
+      n := !n + (Array.length (result.(i).powers ))
     done; !n
   in
   { contracted_shells = result ; size }
@@ -56,7 +56,7 @@ let to_string b =
 "
   ^
   ( Array.map (fun i -> 
-        Contracted_shell.to_string i) b
+        Cs.to_string i) b
     |> Array.to_list
     |> String.concat line
   )
diff --git a/Basis/Basis.mli b/Basis/Basis.mli
index 198ed57..5aa65cb 100644
--- a/Basis/Basis.mli
+++ b/Basis/Basis.mli
@@ -1,10 +1,11 @@
-type t  
+type t = private
+  {
+    (** Number of contracted Gaussians *)
+    size : int;
 
-(** Number of contracted Gaussians *)
-val size : t -> int
-
-(** Array of contracted shells *)
-val contracted_shells : t -> Contracted_shell.t array
+    (** Array of contracted shells *)
+    contracted_shells : Contracted_shell.t array;
+  }
 
 
 (** Returns an array of the basis set per atom *)
diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml
index 98476aa..9fb51f7 100644
--- a/Basis/ContractedShellPair.ml
+++ b/Basis/ContractedShellPair.ml
@@ -18,6 +18,11 @@ type t =
 }
 
 
+module Am = Angular_momentum
+module Co = Coordinate
+module Cs = Contracted_shell
+module Sp = ShellPair
+
 (** Creates an contracted shell pair : an array of pairs of primitive shells.
     A contracted shell with N functions combined with a contracted
     shell with M functions generates a NxM array of shell pairs.
@@ -29,16 +34,15 @@ let create ?cutoff p_a p_b =
     | Some cutoff -> cutoff, -. (log cutoff)
   in
 
-  let center_ab = Coordinate.(
-    Contracted_shell.center p_a |- Contracted_shell.center p_b )
+  let center_ab = Co.( p_a.Cs.center |- p_b.Cs.center )
   in
   let norm_sq =
-    Coordinate.dot center_ab center_ab
+    Co.dot center_ab center_ab
   in
   let norm_coef_scale_a =
-    Contracted_shell.norm_coef_scale p_a 
+    p_a.norm_coef_scale 
   and norm_coef_scale_b =
-    Contracted_shell.norm_coef_scale p_b 
+    p_b.norm_coef_scale 
   in
   let norm_coef_scale =
     Array.map (fun v1 ->
@@ -48,56 +52,53 @@ let create ?cutoff p_a p_b =
     |> Array.concat
   in
   let shell_pairs = 
-    Array.init (Contracted_shell.size p_a) (fun i ->
-        let p_a_expo_center = Coordinate.(
-          Contracted_shell.expo p_a i |. Contracted_shell.center p_a )
-        in
-        let norm_coef_a =
-          Contracted_shell.norm_coef p_a i
-        in
+    Array.init p_a.size (fun i ->
+        let p_a_expo_center = Co.(p_a.Cs.expo.(i) |. p_a.Cs.center) in
+        let norm_coef_a = p_a.norm_coef.(i)  in
 
-        Array.init (Contracted_shell.size p_b) (fun j ->
+        Array.init p_b.size (fun j ->
             try
-              let norm_coef_b =
-                Contracted_shell.norm_coef p_b j
+              let norm_coef_b = p_b.norm_coef.(j) in
+              let norm_coef   = norm_coef_a *. norm_coef_b
               in
-              let norm_coef = 
-                norm_coef_a *. norm_coef_b
-              in
-              if (norm_coef < cutoff) then
+              if norm_coef < cutoff  then
                 raise Null_contribution;
-              let p_b_expo_center = Coordinate.(
-                Contracted_shell.expo p_b j |. Contracted_shell.center p_b )
-              in
-              let expo = Contracted_shell.(expo p_a i +. expo p_b j) in
+              let p_b_expo_center = Co.(p_b.expo.(j) |. p_b.center) in
+              let expo = p_a.expo.(i) +. p_b.expo.(j) in
               let expo_inv = 1. /. expo in
-              let center = Coordinate.(
-                expo_inv |. (p_a_expo_center |+ p_b_expo_center ) )
+              let center = Co.(expo_inv |. (p_a_expo_center |+ p_b_expo_center ) )
               in
               let argexpo =
-                Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv
+                p_a.Cs.expo.(i) *. p_b.Cs.expo.(j) *. norm_sq *. expo_inv
               in
               if (argexpo > log_cutoff) then
                 raise Null_contribution;
               let g =
-                (pi *. expo_inv)**(1.5) *. exp(-. argexpo)
+                (pi *. expo_inv)**(1.5) *. exp (-. argexpo)
               in
               let coef =
-                norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g
+                norm_coef *. p_a.coef.(i) *. p_b.coef.(j) *. g
               in
-              if (abs_float coef  < cutoff) then
+              if abs_float coef  < cutoff then
                 raise Null_contribution;
               let center_a =
-                Coordinate.(center |- Contracted_shell.center p_a)
+                Co.(center |- p_a.center)
               in
               let monocentric =
-                Contracted_shell.center p_a = Contracted_shell.center p_b
+                p_a.center  = p_b.center 
               in
               let totAngMomInt = 
-                (Angular_momentum.to_int (Contracted_shell.totAngMom p_a))
-                + (Angular_momentum.to_int (Contracted_shell.totAngMom p_b))
+                Am.to_int p_a.totAngMom +
+                Am.to_int p_b.totAngMom 
               in
-              Some ShellPair.{ i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; monocentric ; totAngMomInt}
+              Some {
+                Sp.i ; j ;
+                shell_a=p_a ; shell_b=p_b ;
+                norm_coef ; coef ;
+                expo ; expo_inv ;
+                center ; center_a ; center_ab ;
+                norm_sq ; monocentric ; totAngMomInt
+              }
             with
             | Null_contribution -> None
           )
@@ -109,15 +110,15 @@ let create ?cutoff p_a p_b =
     |> List.map (function Some x -> x | None -> assert false)
     |> Array.of_list
   in
-  let coef     = Array.map (fun x -> (fun y -> y.ShellPair.coef) x) shell_pairs
-  and expo_inv = Array.map (fun x -> (fun y -> y.ShellPair.expo_inv) x) shell_pairs
+  let coef     = Array.map (fun x -> (fun y -> y.Sp.coef) x) shell_pairs
+  and expo_inv = Array.map (fun x -> (fun y -> y.Sp.expo_inv) x) shell_pairs
   in
   {
     shell_a = p_a ; shell_b = p_b ; coef ; expo_inv ;
     shell_pairs ; center_ab=shell_pairs.(0).center_ab;
     norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq;
-    totAngMomInt = shell_pairs.(0).ShellPair.totAngMomInt;
-    monocentric =  shell_pairs.(0).ShellPair.monocentric;
+    totAngMomInt = shell_pairs.(0).Sp.totAngMomInt;
+    monocentric =  shell_pairs.(0).Sp.monocentric;
   }
 
 
@@ -159,7 +160,7 @@ let shell_pairs basis =
 
 let equivalent x y =
   (Array.length x = Array.length y) &&
-  (Array.init (Array.length x) (fun k -> ShellPair.equivalent x.(k) y.(k))
+  (Array.init (Array.length x) (fun k -> Sp.equivalent x.(k) y.(k))
   |> Array.fold_left (fun accu x -> x && accu) true)
 
 
diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml
index 5004d19..d47af6a 100644
--- a/Basis/Contracted_shell.ml
+++ b/Basis/Contracted_shell.ml
@@ -1,32 +1,26 @@
 open Util
 open Constants
-open Contracted_shell_type 
 open Coordinate
 
-type t = Contracted_shell_type.t
+type shell_contracted = {
+  expo      : float array;         (* Gaussian exponents *)
+  coef      : float array;         (* Contraction coefficients *)
+  center    : Coordinate.t;        (* Center of all the Gaussians *)
+  totAngMom : Angular_momentum.t;  (* Total angular momentum *)
+  size      : int;                 (* Number of contracted Gaussians *)
+  norm_coef : float array;         (* Normalization coefficient of the class
+                                      corresponding to the i-th contraction *)
+  norm_coef_scale : float array;   (* Inside a class, the norm is the norm
+                                      of the function with (totAngMom,0,0) *. 
+                                      this scaling factor *)
+  index    : int;                  (* Index in the array of contracted shells *)
+  powers   : Zkey.t array;         (* Array of Zkeys corresponding to the
+                                      powers of (x,y,z) in the class *)
+}
 
-let size a = a.size
-let expo a i = a.expo.(i)
-let coef a i = a.coef.(i)
-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 index a  = a.index
-let with_index = Contracted_shell_type.with_index
-let powers a = a.powers
+type t = shell_contracted 
 
-let to_string s =
-  let coord = s.center in
-  let open Printf in
-  (match s.totAngMom with
-   | 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)
-      (get X coord) (get Y coord) (get Z coord) ) ^
-  (Array.map2 (fun e c -> sprintf "%16.8e  %16.8e" e c) s.expo s.coef
-   |> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) 
+module Am = Angular_momentum
 
 (** Normalization coefficient of contracted function i, which depends on the
     exponent and the angular momentum. Two conventions can be chosen : a single
@@ -37,7 +31,7 @@ let to_string s =
 *)
 let compute_norm_coef expo totAngMom =
   let atot =
-    Angular_momentum.to_int totAngMom
+    Am.to_int totAngMom
   in
   let factor int_array = 
        let dfa = Array.map (fun j ->
@@ -61,5 +55,75 @@ let compute_norm_coef expo totAngMom =
   Array.map (fun x -> let f a = x *. (factor a) in f) expo
 
 
-let make = Contracted_shell_type.make
+let make ~index ~expo ~coef ~center ~totAngMom =
+  assert (Array.length expo = Array.length coef);
+  assert (Array.length expo > 0);
+  let norm_coef_func =
+    compute_norm_coef expo totAngMom
+  in
+  let powers =
+    Am.zkey_array (Am.Singlet totAngMom) 
+  in
+  let norm_coef =
+    Array.map (fun f -> f [| Am.to_int totAngMom ; 0 ; 0 |])  norm_coef_func
+  in
+  let norm_coef_scale = 
+    Array.map (fun a ->
+      (norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0)
+    ) powers
+  in
+  { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
+    norm_coef_scale ; powers }
+
+
+let with_index a i =
+  { a with index = i }
+
+
+let to_string s =
+  let coord = s.center in
+  let open Printf in
+  (match s.totAngMom with
+   | Am.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 " (Am.to_string s.totAngMom)
+      (get X coord) (get Y coord) (get Z coord) ) ^
+  (Array.map2 (fun e c -> sprintf "%16.8e  %16.8e" e c) s.expo s.coef
+   |> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) 
+
+
+(** Normalization coefficient of contracted function i, which depends on the
+    exponent and the angular momentum. Two conventions can be chosen : a single
+    normalisation factor for all functions of the class, or a coefficient which
+    depends on the powers of x,y and z.
+    Returns, for each contracted function, an array of functions taking as
+    argument the [|x;y;z|] powers.
+*)
+let compute_norm_coef expo totAngMom =
+  let atot =
+    Am.to_int totAngMom
+  in
+  let factor int_array = 
+       let dfa = Array.map (fun j ->
+           (float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
+         ) int_array
+       in
+       sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
+  in
+  let expo = 
+    if atot mod 2 = 0 then
+      Array.map (fun alpha ->
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
+      ) expo
+    else
+      Array.map (fun alpha ->
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
+      ) expo
+  in
+  Array.map (fun x -> let f a = x *. factor a in f) expo
+
+
 
diff --git a/Basis/Contracted_shell.mli b/Basis/Contracted_shell.mli
index 55bf153..6f1a73b 100644
--- a/Basis/Contracted_shell.mli
+++ b/Basis/Contracted_shell.mli
@@ -1,45 +1,16 @@
-type t = Contracted_shell_type.t
+type shell_contracted = private {
+  expo      : float array;
+  coef      : float array;
+  center    : Coordinate.t;
+  totAngMom : Angular_momentum.t;
+  size      : int;
+  norm_coef : float array;
+  norm_coef_scale : float array;
+  index    : 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
+type t = shell_contracted
 
 
 (** Pretty-printing of the contracted shell in a string *)
@@ -52,4 +23,5 @@ val make :
   coef:float array ->
   center:Coordinate.t -> totAngMom:Angular_momentum.t -> t
 
-
+(** Returns a copy of the contracted shell with a modified index *)
+val with_index  : t -> int -> t
diff --git a/Basis/Contracted_shell_type.ml b/Basis/Contracted_shell_type.ml
deleted file mode 100644
index ab7d494..0000000
--- a/Basis/Contracted_shell_type.ml
+++ /dev/null
@@ -1,73 +0,0 @@
-open Util
-open Constants
-
-type shell_contracted = {
-  expo      : float array;
-  coef      : float array;
-  center    : Coordinate.t;
-  totAngMom : Angular_momentum.t;
-  size      : int;
-  norm_coef : float array;
-  norm_coef_scale : float array;
-  index    : int;
-  powers   : Zkey.t array;
-}
-
-type t = shell_contracted 
-
-(** Normalization coefficient of contracted function i, which depends on the
-    exponent and the angular momentum. Two conventions can be chosen : a single
-    normalisation factor for all functions of the class, or a coefficient which
-    depends on the powers of x,y and z.
-    Returns, for each contracted function, an array of functions taking as
-    argument the [|x;y;z|] powers.
-*)
-let compute_norm_coef expo totAngMom =
-  let atot =
-    Angular_momentum.to_int totAngMom
-  in
-  let factor int_array = 
-       let dfa = Array.map (fun j ->
-           ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
-         ) int_array
-       in
-       sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
-  in
-  let expo = 
-    if atot mod 2 = 0 then
-      Array.map (fun alpha ->
-        let alpha_2 = alpha +. alpha in
-        (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
-      ) expo
-    else
-      Array.map (fun alpha ->
-        let alpha_2 = alpha +. alpha in
-        (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
-      ) expo
-  in
-  Array.map (fun x -> let f a = x *. (factor a) in f) expo
-
-
-let make ~index ~expo ~coef ~center ~totAngMom =
-  assert (Array.length expo = Array.length coef);
-  assert (Array.length expo > 0);
-  let norm_coef_func =
-    compute_norm_coef expo totAngMom
-  in
-  let powers =
-    Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom) 
-  in
-  let norm_coef =
-    Array.map (fun f -> f [| Angular_momentum.to_int totAngMom ; 0 ; 0 |])  norm_coef_func
-  in
-  let norm_coef_scale = 
-    Array.map (fun a ->
-      (norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0)
-    ) powers
-  in
-  { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
-    norm_coef_scale ; powers }
-
-
-let with_index a i = { a with index = i }
-
diff --git a/Basis/Contracted_shell_type.mli b/Basis/Contracted_shell_type.mli
deleted file mode 100644
index 1ec55fb..0000000
--- a/Basis/Contracted_shell_type.mli
+++ /dev/null
@@ -1,24 +0,0 @@
-type shell_contracted = private {
-  expo      : float array;
-  coef      : float array;
-  center    : Coordinate.t;
-  totAngMom : Angular_momentum.t;
-  size      : int;
-  norm_coef : float array;
-  norm_coef_scale : float array;
-  index    : int;
-  powers   : Zkey.t array;
-}
-
-type t = shell_contracted
-
-
-val make :
-  index:int ->
-  expo:float array ->
-  coef:float array ->
-  center:Coordinate.t -> totAngMom:Angular_momentum.t -> t
-
-
-(** Returns a copy of the contracted shell with a modified index *)
-val with_index  : t -> int -> t
diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index 7d67fb0..044aca5 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -6,6 +6,10 @@ open Bigarray
 
 type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t
 
+module Bs  = Basis
+module Cs  = Contracted_shell
+module Csp = ContractedShellPair
+
 (** (00|00)^m : Fundamental electron repulsion integral
     $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $
 
@@ -77,16 +81,14 @@ let of_basis basis =
     | _ -> assert false
   in
 
-  let n =
-    Basis.size basis
-  and shell =
-    Basis.contracted_shells basis
+  let n = basis.Bs.size
+  and shell = basis.Bs.contracted_shells 
   in
 
 
   (* Pre-compute all shell pairs *)
   let shell_pairs =
-    ContractedShellPair.shell_pairs shell
+    Csp.shell_pairs shell
   in
 
   (* Pre-compute diagonal integrals for Schwartz *)
@@ -103,7 +105,7 @@ let of_basis basis =
 
   let icount = ref 0 in
   for i=0 to (Array.length shell) - 1 do
-    print_int (Contracted_shell.index shell.(i)) ; print_newline ();
+    print_int shell.(i).index ; print_newline ();
     for j=0 to i do
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
         if (schwartz_p_max >= cutoff) then
@@ -112,6 +114,9 @@ let of_basis basis =
   done;
   Printf.printf "%d shell pairs computed in %f seconds\n" !icount (Unix.gettimeofday () -. t0);
 
+  (* Group shell pairs by common pairs of atoms *)
+  
+
   (* 4D data initialization *)
   let eri_array = 
     Genarray.create Float32 fortran_layout [| n ; n ; n ; n|]
@@ -124,7 +129,7 @@ let of_basis basis =
   let inn = ref 0 and out = ref 0 in
 
   for i=0 to (Array.length shell) - 1 do
-    print_int (Contracted_shell.index shell.(i)) ; print_newline ();
+    print_int shell.(i).index ; print_newline ();
     for j=0 to i do
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
       try
@@ -135,7 +140,7 @@ let of_basis basis =
         in
 
         let sp =
-          shell_p.ContractedShellPair.shell_pairs
+          shell_p.Csp.shell_pairs
         in
 
         for k=0 to i do
@@ -148,7 +153,7 @@ let of_basis basis =
               shell_q = shell_pairs.(k).(l)
             in
             let sq =
-              shell_q.ContractedShellPair.shell_pairs
+              shell_q.Csp.shell_pairs
             in
 
             let swap = 
@@ -172,16 +177,16 @@ let of_basis basis =
 
             (* Write the data in the output file *)
             Array.iteri (fun i_c powers_i ->
-                let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in
+                let i_c = shell.(i).index + i_c + 1 in
                 let xi = to_powers powers_i in
                 Array.iteri (fun j_c powers_j ->
-                    let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in
+                    let j_c = shell.(j).index + j_c + 1 in
                     let xj = to_powers powers_j in
                     Array.iteri (fun k_c powers_k ->
-                        let k_c = (Contracted_shell.index shell.(k)) + k_c + 1 in 
+                        let k_c = shell.(k).index + k_c + 1 in 
                         let xk = to_powers powers_k in
                         Array.iteri (fun l_c powers_l ->
-                            let l_c = (Contracted_shell.index shell.(l)) + l_c + 1 in
+                            let l_c = shell.(l).index + l_c + 1 in
                             let xl = to_powers powers_l in
                             let key = 
                               if swap then
@@ -205,10 +210,10 @@ let of_basis basis =
                               )
                             else
                               out := !out + 1;
-                          ) (Contracted_shell.powers shell.(l)) 
-                      ) (Contracted_shell.powers shell.(k))
-                  ) (Contracted_shell.powers shell.(j))
-              ) (Contracted_shell.powers shell.(i));
+                          ) shell.(l).powers 
+                      ) shell.(k).powers
+                  ) shell.(j).powers
+              ) shell.(i).powers
             with NullIntegral ->  ()
           done;
         done;
@@ -264,7 +269,7 @@ let to_file ~filename basis =
   let oc = open_out filename in
   let zkey = Array.map (fun b -> 
       let result =
-        Angular_momentum.(zkey_array (Kind_1 b.Contracted_shell.totAngMom))
+        Angular_momentum.(zkey_array (Kind_1 b.totAngMom))
       in
       { n=Array.length result ;  cls=result }
     ) basis
@@ -362,7 +367,7 @@ let to_file ~filename basis =
 let xto_file ~filename basis =
   let zkey = Array.map (fun b -> 
       let result =
-        Angular_momentum.(zkey_array (Kind_1 b.Contracted_shell.totAngMom))
+        Angular_momentum.(zkey_array (Kind_1 b.totAngMom))
       in
       { n=Array.length result ;  cls=result }
     ) basis
@@ -373,10 +378,10 @@ let xto_file ~filename basis =
 
   let (i,j,k,l) = (1,1,1,18) in
   let (i,j,k,l) = (i-1,j-1,k-1,l-1) in
-  basis.(i) |> Contracted_shell.to_string |> print_endline;
-  basis.(j) |> Contracted_shell.to_string |> print_endline;
-  basis.(k) |> Contracted_shell.to_string |> print_endline;
-  basis.(l) |> Contracted_shell.to_string |> print_endline;
+  basis.(i) |> Cs.to_string |> print_endline;
+  basis.(j) |> Cs.to_string |> print_endline;
+  basis.(k) |> Cs.to_string |> print_endline;
+  basis.(l) |> Cs.to_string |> print_endline;
   let bi, bj, bk, bl =
     basis.(i), basis.(j), basis.(k), basis.(l)
   in
diff --git a/Basis/General_basis.ml b/Basis/General_basis.ml
index edb20d3..752d9b1 100644
--- a/Basis/General_basis.ml
+++ b/Basis/General_basis.ml
@@ -9,6 +9,8 @@ type general_contracted_shell = Angular_momentum.t * (primitive array)
 type t = Element.t * (general_contracted_shell array)
 
 
+module Am = Angular_momentum
+
 let string_of_primitive ?id prim =
   match id with
   | None -> (string_of_float prim.exponent)^" "^(string_of_float prim.coefficient)
@@ -20,7 +22,7 @@ let string_of_contracted_shell (angular_momentum, prim_array) =
     Array.length prim_array
   in
   Printf.sprintf "%s %d\n%s"
-    (Angular_momentum.to_string angular_momentum) n
+    (Am.to_string angular_momentum) n
     (Array.init n (fun i -> string_of_primitive ~id:(i+1) prim_array.(i))
      |> Array.to_list
      |> String.concat "\n")
diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
index 0200ddb..78c9325 100644
--- a/Basis/KinInt.ml
+++ b/Basis/KinInt.ml
@@ -1,8 +1,13 @@
 open Util
 open Constants
 open Lacaml.D
-open Powers
-open Coordinate
+
+module Am  = Angular_momentum
+module Bs  = Basis
+module Co  = Coordinate
+module Cs  = Contracted_shell
+module Csp = ContractedShellPair
+module Sp  = ShellPair
 
 type t = Mat.t
 
@@ -11,13 +16,12 @@ type t = Mat.t
 let contracted_class shell_a shell_b : float Zmap.t = 
 
   let shell_p =
-    ContractedShellPair.create shell_a shell_b
+    Csp.create shell_a shell_b
   in
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
-    Angular_momentum.zkey_array (Angular_momentum.Doublet
-                                   Contracted_shell.(totAngMom shell_a, totAngMom shell_b))
+    Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom))
   in
 
   let contracted_class = 
@@ -26,35 +30,35 @@ let contracted_class shell_a shell_b : float Zmap.t =
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
-  let sp = shell_p.ContractedShellPair.shell_pairs in
+  let sp = shell_p.Csp.shell_pairs in
   let center_ab = 
-    shell_p.ContractedShellPair.center_ab
+    shell_p.Csp.center_ab
   in
   let norm_coef_scale =
-    shell_p.ContractedShellPair.norm_coef_scale
+    shell_p.Csp.norm_coef_scale
   in
 
   for ab=0 to (Array.length sp - 1)
   do
     let coef_prod =
-      shell_p.ContractedShellPair.coef.(ab) 
+      shell_p.Csp.coef.(ab) 
     in
     (** Screening on thr product of coefficients *)
     if (abs_float coef_prod) > 1.e-4*.cutoff then
       begin
         let center_pa = 
-          sp.(ab).ShellPair.center_a
+          sp.(ab).Sp.center_a
         in
         let expo_inv = 
-          shell_p.ContractedShellPair.expo_inv.(ab)
+          shell_p.Csp.expo_inv.(ab)
         in
         let i, j = 
-          sp.(ab).ShellPair.i, sp.(ab).ShellPair.j 
+          sp.(ab).Sp.i, sp.(ab).Sp.j 
         in
         let expo_a = 
-          Contracted_shell.expo sp.(ab).ShellPair.shell_a i
+          sp.(ab).Sp.shell_a.expo.(i)
         and expo_b = 
-          Contracted_shell.expo sp.(ab).ShellPair.shell_b j
+          sp.(ab).Sp.shell_b.expo.(j)
         in
 
         Array.iteri (fun i key ->
@@ -65,14 +69,14 @@ let contracted_class shell_a shell_b : float Zmap.t =
             in
             let ov a b k = 
               let xyz = match k with
-              | 0 -> X
-              | 1 -> Y
-              | _ -> Z
+              | 0 -> Co.X
+              | 1 -> Co.Y
+              | _ -> Co.Z
               in
               Overlap_primitives.hvrr (a, b)
                 expo_inv
-                (Coordinate.get xyz center_ab,
-                 Coordinate.get xyz center_pa) 
+                (Co.get xyz center_ab,
+                 Co.get xyz center_pa) 
             in
             let f k = 
               ov angMomA.(k) angMomB.(k) k
@@ -116,10 +120,8 @@ let of_basis basis =
     | _ -> assert false
   in
 
-  let n =
-    Basis.size basis
-  and shell =
-    Basis.contracted_shells basis
+  let n     = basis.Bs.size 
+  and shell = basis.Bs.contracted_shells
   in
 
   let result = Mat.create n n in
@@ -131,10 +133,10 @@ let of_basis basis =
       in
 
       Array.iteri (fun j_c powers_j ->
-          let j_c = Contracted_shell.index shell.(j) + j_c + 1 in
+          let j_c = shell.(j).index + j_c + 1 in
           let xj = to_powers powers_j in
           Array.iteri (fun i_c powers_i ->
-              let i_c = Contracted_shell.index shell.(i) + i_c + 1 in
+              let i_c = shell.(i).index + i_c + 1 in
               let xi = to_powers powers_i in
               let key =
                 Zkey.of_powers (Zkey.Six (xi,xj))
@@ -145,8 +147,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) (Contracted_shell.powers shell.(i));
-        ) (Contracted_shell.powers shell.(j))
+            ) shell.(i).powers
+        ) shell.(j).powers 
     done;
   done;
   Mat.detri result;
diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index f5e0e69..feacf29 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -6,6 +6,8 @@ open Lacaml.D
 
 type t = Mat.t
 
+module Bs = Basis
+
 (** (0|0)^m : Fundamental electron-nucleus repulsion integral
     $ \int \phi_p(r1) 1/r_{C} dr_1 $
 
@@ -43,10 +45,8 @@ let of_basis_nuclei basis nuclei =
     | _ -> assert false
   in
 
-  let n =
-    Basis.size basis
-  and shell =
-    Basis.contracted_shells basis
+  let n     = basis.Bs.size 
+  and shell = basis.Bs.contracted_shells 
   in
 
   let eni_array = Mat.create n n in
@@ -71,10 +71,10 @@ let of_basis_nuclei basis nuclei =
 
       (* Write the data in the output file *)
       Array.iteri (fun i_c powers_i ->
-          let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in
+          let i_c = shell.(i).index + i_c + 1 in
           let xi = to_powers powers_i in
           Array.iteri (fun j_c powers_j ->
-              let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in
+              let j_c = shell.(j).index + j_c + 1 in
               let xj = to_powers powers_j in
               let key = 
                 Zkey.of_powers (Zkey.Six (xi,xj))
@@ -84,8 +84,8 @@ let of_basis_nuclei basis nuclei =
               in
               eni_array.{j_c,i_c} <- value;
               eni_array.{i_c,j_c} <- value;
-            ) (Contracted_shell.powers shell.(j))
-        ) (Contracted_shell.powers shell.(i));
+            ) shell.(j).powers
+        ) shell.(i).powers 
     done;
   done;
   Mat.detri eni_array;
diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml
index bd4c114..065240b 100644
--- a/Basis/OneElectronRR.ml
+++ b/Basis/OneElectronRR.ml
@@ -1,10 +1,15 @@
 open Util
 open Constants
-open Powers
-open Coordinate
 
 exception NullPair
 
+module Am  = Angular_momentum
+module Co  = Coordinate
+module Cs  = Contracted_shell
+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.
@@ -13,29 +18,28 @@ let chop f g =
 
 
 (** Horizontal and Vertical Recurrence Relations (HVRR) *)
-let hvrr_one_e (angMom_a, angMom_b) 
-    zero_m_array (expo_b) (expo_inv_p) (center_ab, center_pa, center_pc)
-    map
-  =
+let hvrr_one_e  angMom_a  angMom_b  
+    zero_m_array expo_b expo_inv_p
+    center_ab center_pa center_pc
+    map =
 
-  let maxm = angMom_a.tot + angMom_b.tot  in
+  let maxm = angMom_a.Po.tot + angMom_b.Po.tot  in
   let maxsze = maxm+1 in
-  let empty = Array.make maxsze 0. in
 
   let get_xyz angMom =
     match angMom with
-    | { y=0 ; z=0 ; _ } -> X
-    | {       z=0 ; _ } -> Y
-    | _                 -> Z
+    | { Po.y=0 ; z=0 ; _ } -> Co.X
+    | {          z=0 ; _ } -> Co.Y
+    | _                    -> Co.Z
   in
 
 
   (** Vertical recurrence relations *)
   let rec vrr angMom_a =
-    let { x=ax ; y=ay ; z=az } = angMom_a in
+    let { Po.x=ax ; y=ay ; z=az } = angMom_a in
     if ax < 0 || ay < 0 || az < 0 then raise Exit
     else
-      match angMom_a.tot  with
+      match angMom_a.Po.tot  with
       | 0 -> zero_m_array
       | _ ->
         let key = Zkey.of_powers (Zkey.Three angMom_a) in
@@ -43,35 +47,26 @@ let hvrr_one_e (angMom_a, angMom_b)
         try Zmap.find map key with
         | Not_found -> 
           let result = 
-            let xyz = get_xyz angMom_a in
-            let am = Powers.decr xyz angMom_a in
-            let amxyz = Powers.get xyz am in
-            if amxyz < 0 then empty else
-              let f1 = Coordinate.get xyz center_pa 
-              and f2 = expo_inv_p *.  (Coordinate.get xyz center_pc)
-              in
-              if amxyz < 1 then
-                let v1 =
-                  vrr am
-                in
-                Array.init maxsze (fun m ->
-                    if m = maxm then (f1 *. v1.(m) ) else
-                    (f1 *. v1.(m) ) -. f2 *. v1.(m+1) )
-              else
+            let xyz   = get_xyz angMom_a in
+            let am    = Po.decr xyz angMom_a in
+            let amxyz = Po.get xyz am in
+            let f1    = Co.get xyz center_pa  in
+            let f2    = expo_inv_p *. Co.get xyz center_pc in
+            if amxyz < 1 then
+                let v1 = vrr am in
+                Array.init (maxsze - angMom_a.Po.tot)
+                  (fun m -> f1 *. v1.(m) -. f2 *. v1.(m+1))
+            else
                 let v3 =
-                  let amm = Powers.decr xyz am in
+                  let amm = Po.decr xyz am in
                   vrr amm
                 in
-                let v1 =
-                  vrr am
-                in
-                let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
-                Array.init maxsze (fun m -> f1 *. v1.(m) -.
-                    (if m = maxm then 0. else
-                        f2 *. v1.(m+1) )
-                    +. f3 *. (v3.(m) -.  if m = maxm then 0. else
-                                expo_inv_p *. v3.(m+1))
-                  )
+                let v1 = vrr am in
+                let f3 = float_of_int amxyz *. expo_inv_p *. 0.5 in
+                Array.init (maxsze - angMom_a.Po.tot)
+                  (fun m -> f1 *. v1.(m) -.  f2 *. v1.(m+1) +.
+                            f3 *. (v3.(m) -. expo_inv_p *. v3.(m+1))
+                )
           in Zmap.add map key result;
           result
 
@@ -79,28 +74,19 @@ let hvrr_one_e (angMom_a, angMom_b)
   (** Horizontal recurrence relations *)
   and hrr angMom_a angMom_b =
 
-    let { x=bx ; y=by ; z=bz } = angMom_b in
-    if bx < 0 || by < 0 || bz < 0 then raise Exit
-    else
-      match angMom_b.tot  with
-      | 0 -> (vrr angMom_a).(0)
-      | _ ->
-        let xyz = get_xyz angMom_b in
-        let bxyz = Powers.get xyz angMom_b in
-        if (bxyz < 1) then 0. else
-          let ap = Powers.incr xyz angMom_a in
-          let bm = Powers.decr xyz angMom_b in
-          let h1 =
-            hrr ap bm
-          in
-          let f2 =
-            Coordinate.get xyz center_ab
-          in
+    match angMom_b.Po.tot  with
+    | 0 -> (vrr angMom_a).(0)
+    | _ ->
+      let xyz = get_xyz angMom_b in
+      let bxyz = Po.get xyz angMom_b in
+      if (bxyz < 1) then 0. else
+          let ap = Po.incr xyz angMom_a in
+          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
-            let h2 =
-              hrr angMom_a bm
-            in
-            h1 +. f2 *. h2
+              let h2 = hrr angMom_a bm in
+              h1 +. f2 *. h2
 
   in
   hrr angMom_a angMom_b
@@ -116,19 +102,16 @@ let hvrr_one_e (angMom_a, angMom_b)
 (** Computes all the one-electron integrals of the contracted shell pair *)
 let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
 
-  let shell_a = shell_p.ContractedShellPair.shell_a
-  and shell_b = shell_p.ContractedShellPair.shell_b
+  let shell_a = shell_p.Csp.shell_a
+  and shell_b = shell_p.Csp.shell_b
   in
   let maxm = 
-    let open Angular_momentum in
-    (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) 
+    (Am.to_int @@ shell_a.totAngMom) + (Am.to_int @@ shell_b.totAngMom) 
   in
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
-    Angular_momentum.zkey_array
-      (Angular_momentum.Doublet
-         Contracted_shell.(totAngMom shell_a, totAngMom shell_b))
+    Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom))
   in
 
   let contracted_class = 
@@ -137,52 +120,50 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
-  let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale 
+  let norm_coef_scale_p = shell_p.Csp.norm_coef_scale 
   in
-  for ab=0 to (Array.length shell_p.ContractedShellPair.shell_pairs - 1)
+  for ab=0 to (Array.length shell_p.Csp.shell_pairs - 1)
   do
-    let b = shell_p.ContractedShellPair.shell_pairs.(ab).ShellPair.j in
+    let b = shell_p.Csp.shell_pairs.(ab).Sp.j in
     try
       begin
-        let coef_prod = shell_p.ContractedShellPair.coef.(ab) in
+        let coef_prod = shell_p.Csp.coef.(ab) in
 
         (** Screening on the product of coefficients *)
-        if (abs_float coef_prod) < 1.e-4*.cutoff then
+        if abs_float coef_prod < 1.e-3 *. cutoff then
           raise NullPair;
 
 
         let expo_pq_inv =
-          shell_p.ContractedShellPair.expo_inv.(ab) 
+          shell_p.Csp.expo_inv.(ab) 
         in
 
         let center_ab = 
-          shell_p.ContractedShellPair.center_ab
+          shell_p.Csp.center_ab
         in
         let center_p = 
-          shell_p.ContractedShellPair.shell_pairs.(ab).ShellPair.center
+          shell_p.Csp.shell_pairs.(ab).Sp.center
         in
         let center_pa = 
-          Coordinate.(center_p |- Contracted_shell.center shell_a) 
+          Co.(center_p |- shell_a.center) 
         in
 
         for c=0 to Array.length geometry - 1 do
           let element, nucl_coord = geometry.(c) in
           let charge = Element.to_charge element |> Charge.to_float in
           let center_pc =
-            Coordinate.(center_p |- nucl_coord ) 
+            Co.(center_p |- nucl_coord ) 
           in
           let norm_pq_sq = 
-            Coordinate.dot center_pc center_pc
+            Co.dot center_pc center_pc
           in
 
           let zero_m_array = 
             zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
           in
-          match Contracted_shell.(totAngMom shell_a, totAngMom shell_b) with
-          | Angular_momentum.(S,S) -> 
-            let integral =
-              zero_m_array.(0)
-            in
+          match (shell_a.totAngMom, shell_b.totAngMom) with
+          | Am.(S,S) -> 
+            let integral = zero_m_array.(0) in
             contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge
           | _ ->
             let map = Zmap.create (2*maxm) in
@@ -198,11 +179,12 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
                 let norm = norm_coef_scale.(i) in
                 let coef_prod = coef_prod *. norm in
                 let integral = 
-                  hvrr_one_e (angMomA, angMomB)
-                     zero_m_array
-                    (Contracted_shell.expo shell_b b)
-                    (shell_p.ContractedShellPair.expo_inv.(ab))
-                    (center_ab, center_pa, center_pc)
+                  hvrr_one_e
+                    angMomA  angMomB
+                    zero_m_array
+                    shell_b.expo.(b)
+                    shell_p.Csp.expo_inv.(ab)
+                    center_ab center_pa center_pc
                     map
                 in
                 contracted_class.(i) <- contracted_class.(i) -. coef_prod *. integral *. charge
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index cb363b8..b312db3 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -1,22 +1,26 @@
 open Util
 open Constants
 open Lacaml.D
-open Coordinate
 
 type t = Mat.t
 
+module Am  = Angular_momentum
+module Bs  = Basis
+module Co  = Coordinate
+module Cs  = Contracted_shell
+module Csp = ContractedShellPair
+module Sp  = ShellPair
 
 (** Computes all the overlap integrals of the contracted shell pair *)
 let contracted_class shell_a shell_b : float Zmap.t = 
 
   let shell_p =
-    ContractedShellPair.create shell_a shell_b
+    Csp.create shell_a shell_b
   in
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
-    Angular_momentum.zkey_array (Angular_momentum.Doublet
-                                   Contracted_shell.(totAngMom shell_a, totAngMom shell_b))
+    Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom))
   in
 
   let contracted_class = 
@@ -24,13 +28,13 @@ let contracted_class shell_a shell_b : float Zmap.t =
   in
 
   let sp =
-    shell_p.ContractedShellPair.shell_pairs
+    shell_p.Csp.shell_pairs
   in
   let center_ab = 
-    shell_p.ContractedShellPair.center_ab
+    shell_p.Csp.center_ab
   in
   let norm_coef_scale =
-    shell_p.ContractedShellPair.norm_coef_scale
+    shell_p.Csp.norm_coef_scale
   in
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
@@ -38,16 +42,16 @@ let contracted_class shell_a shell_b : float Zmap.t =
   for ab=0 to (Array.length sp - 1)
   do
     let coef_prod =
-      shell_p.ContractedShellPair.coef.(ab) 
+      shell_p.Csp.coef.(ab) 
     in
     (** Screening on thr product of coefficients *)
     if (abs_float coef_prod) > 1.e-4*.cutoff then
       begin
         let expo_inv = 
-          shell_p.ContractedShellPair.expo_inv.(ab)
+          shell_p.Csp.expo_inv.(ab)
         in
         let center_pa = 
-          sp.(ab).ShellPair.center_a
+          sp.(ab).Sp.center_a
         in
 
         Array.iteri (fun i key ->
@@ -58,14 +62,14 @@ let contracted_class shell_a shell_b : float Zmap.t =
             in
             let f k = 
               let xyz = match k with
-              | 0 -> X
-              | 1 -> Y
-              | _ -> Z
+              | 0 -> Co.X
+              | 1 -> Co.Y
+              | _ -> Co.Z
               in
               Overlap_primitives.hvrr (angMomA.(k), angMomB.(k)) 
                 expo_inv
-                (Coordinate.get xyz center_ab,
-                 Coordinate.get xyz center_pa) 
+                (Co.get xyz center_ab,
+                 Co.get xyz center_pa) 
             in
             let norm =  norm_coef_scale.(i) in
             let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in 
@@ -89,10 +93,8 @@ let of_basis basis =
     | _ -> assert false
   in
 
-  let n =
-    Basis.size basis
-  and shell =
-    Basis.contracted_shells basis
+  let n     = basis.Bs.size 
+  and shell = basis.Bs.contracted_shells 
   in
 
   let result = Mat.create n n in
@@ -104,10 +106,10 @@ let of_basis basis =
       in
 
       Array.iteri (fun j_c powers_j ->
-          let j_c = Contracted_shell.index shell.(j) + j_c + 1 in
+          let j_c = shell.(j).index + j_c + 1 in
           let xj = to_powers powers_j in
           Array.iteri (fun i_c powers_i ->
-              let i_c = Contracted_shell.index shell.(i) + i_c + 1 in
+              let i_c = shell.(i).index + i_c + 1 in
               let xi = to_powers powers_i in
               let key = 
                 Zkey.of_powers (Zkey.Six (xi,xj))
@@ -118,8 +120,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) (Contracted_shell.powers shell.(i));
-        ) (Contracted_shell.powers shell.(j))
+            ) shell.(i).powers
+        ) shell.(j).powers
     done;
   done;
   Mat.detri result;
diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml
index edbe122..3fdd502 100644
--- a/Basis/TwoElectronRR.ml
+++ b/Basis/TwoElectronRR.ml
@@ -1,7 +1,12 @@
 open Util
 open Constants
-open Powers
-open Coordinate
+
+module Am  = Angular_momentum
+module Co  = Coordinate
+module Cs  = Contracted_shell
+module Csp = ContractedShellPair
+module Sp  = ShellPair
+module Po  = Powers
 
 let debug=false
 
@@ -10,43 +15,44 @@ let cutoff2 = cutoff *. cutoff
 exception NullQuartet
 
 (** Horizontal and Vertical Recurrence Relations (HVRR) *)
-let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
+let rec hvrr_two_e
+    angMom_a angMom_b angMom_c angMom_d
     zero_m_array
-    (expo_b, expo_d)
-    (expo_inv_p, expo_inv_q)
-    (center_ab, center_cd, center_pq)
-    (center_pa, center_qc)
-    map_1d map_2d 
-  =
+    expo_b expo_d
+    expo_inv_p expo_inv_q
+    center_ab center_cd center_pq
+    center_pa center_qc
+    map_1d map_2d =
 
   (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *)
-  if  angMom_a.tot + angMom_b.tot < angMom_c.tot + angMom_d.tot then
-    hvrr_two_e (angMom_c, angMom_d, angMom_a, angMom_b)
+  if  angMom_a.Po.tot + angMom_b.Po.tot < angMom_c.Po.tot + angMom_d.Po.tot then
+    hvrr_two_e
+      angMom_c angMom_d angMom_a angMom_b
       zero_m_array
-      (expo_d, expo_b)
-      (expo_inv_q, expo_inv_p)
-      (center_cd, center_ab, (Coordinate.neg center_pq) )
-      (center_qc, center_pa)
+      expo_d expo_b
+      expo_inv_q expo_inv_p
+      center_cd center_ab (Co.neg center_pq) 
+      center_qc center_pa
       map_1d map_2d
 
   else
 
-    let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in
+    let maxm = angMom_a.Po.tot + angMom_b.Po.tot + angMom_c.Po.tot + angMom_d.Po.tot in
     let maxsze = maxm+1 in
 
 
     let get_xyz angMom = 
       match angMom with
-      | { y=0 ; z=0 ; _ } -> X
-      | {       z=0 ; _ } -> Y
-      | _                 -> Z
+      | { Po.y=0 ; z=0 ; _ } -> Co.X
+      | {          z=0 ; _ } -> Co.Y
+      | _                    -> Co.Z
     in
 
 
     (** Vertical recurrence relations *)
     let rec vrr0 angMom_a = 
 
-      match angMom_a.tot with 
+      match angMom_a.Po.tot with 
       | 0 -> zero_m_array
       | _ ->
         let key = Zkey.of_powers_three angMom_a in
@@ -55,13 +61,13 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
         | Not_found ->
           let result = 
             let xyz = get_xyz angMom_a in
-            let am    = Powers.decr xyz angMom_a in
-            let amxyz = Powers.get xyz am in
+            let am    = Po.decr xyz angMom_a in
+            let amxyz = Po.get xyz am in
 
-            let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) 
-            and f2 = expo_b *. expo_inv_p *.  (Coordinate.get xyz center_ab)
+            let f1 = expo_inv_p *. Co.get xyz center_pq
+            and f2 = expo_b *. expo_inv_p *. Co.get xyz center_ab
             in
-            let result = Array.create_float (maxsze-angMom_a.tot) in
+            let result = Array.create_float (maxsze - angMom_a.Po.tot) in
             if amxyz = 0 then
                 begin
                   let v1 = vrr0 am in
@@ -70,7 +76,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
                 end
             else
                 begin
-                  let amm = Powers.decr xyz am in
+                  let amm = Po.decr xyz am in
                   let v3 = vrr0 amm in
                   let v1 = vrr0 am in
                   let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
@@ -85,7 +91,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
 
     and vrr angMom_a angMom_c =
 
-      match angMom_a.tot, angMom_c.tot with
+      match angMom_a.Po.tot, angMom_c.Po.tot with
       | (i,0) -> if (i>0) then vrr0 angMom_a
                  else zero_m_array 
       | (_,_) ->
@@ -94,21 +100,21 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
         try Zmap.find map_2d key with
         | Not_found -> 
           let result = 
-            (* angMom_c.tot > 0  so cm.tot >= 0 *)
+            (* angMom_c.Po.tot > 0  so cm.Po.tot >= 0 *)
             let xyz   = get_xyz angMom_c in
-            let cm    = Powers.decr xyz angMom_c in
-            let cmxyz = Powers.get xyz cm in
-            let axyz  = Powers.get xyz angMom_a in
+            let cm    = Po.decr xyz angMom_c in
+            let cmxyz = Po.get xyz cm in
+            let axyz  = Po.get xyz angMom_a in
 
             let f1 =
-              -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) 
+              -. expo_d *. expo_inv_q *. Co.get xyz center_cd
             and f2 = 
-              expo_inv_q *. (Coordinate.get xyz center_pq)
+              expo_inv_q *. Co.get xyz center_pq
             in
-            let result = Array.make (maxsze - angMom_a.tot - angMom_c.tot) 0. in
+            let result = Array.make (maxsze - angMom_a.Po.tot - angMom_c.Po.tot) 0. in
             if axyz > 0 then 
                 begin
-                  let am = Powers.decr xyz angMom_a in
+                  let am = Po.decr xyz angMom_a in
                   let f5 = 
                     (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 
                   in
@@ -128,7 +134,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
                       (abs_float (f3 *. expo_inv_q) > cutoff) then 
                     begin
                       let v3 = 
-                        let cmm = Powers.decr xyz cm in
+                        let cmm = Po.decr xyz cm in
                         vrr angMom_a cmm 
                       in
                       Array.iteri (fun m _ ->
@@ -151,7 +157,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
     (*
     and trr angMom_a angMom_c =
 
-      match (angMom_a.tot, angMom_c.tot) with
+      match (angMom_a.Po.tot, angMom_c.Po.tot) with
       | (i,0) -> if (i>0) then (vrr0 angMom_a).(0)
                  else zero_m_array.(0)
       | (_,_) ->
@@ -161,14 +167,14 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
         | Not_found -> 
           let result = 
             let xyz   = get_xyz angMom_c in
-            let axyz  = Powers.get  xyz angMom_a in
-            let cm    = Powers.decr xyz angMom_c in
-            let cmxyz = Powers.get  xyz cm       in
+            let axyz  = Po.get  xyz angMom_a in
+            let cm    = Po.decr xyz angMom_c in
+            let cmxyz = Po.get  xyz cm       in
 
             let expo_inv_q_over_p = expo_inv_q /. expo_inv_p in
             let f =
-              Coordinate.get xyz center_qc +. expo_inv_q_over_p *.
-              (Coordinate.get xyz center_pa)
+              Co.get xyz center_qc +. expo_inv_q_over_p *.
+              Co.get xyz center_pa
             in
             let result = 0. in
 
@@ -176,7 +182,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
               if cmxyz < 1 then result else
                   let f = 0.5 *. (float_of_int cmxyz) *. expo_inv_q in
                   if abs_float f < cutoff then 0.  else
-                      let cmm = Powers.decr xyz cm in
+                      let cmm = Po.decr xyz cm in
                       let v3 = trr angMom_a cmm in
                       result +. f *. v3
             in
@@ -188,7 +194,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
             let result = 
               if cmxyz < 0 then result else
                   let f = -. expo_inv_q_over_p in
-                  let ap = Powers.incr xyz angMom_a in
+                  let ap = Po.incr xyz angMom_a in
                   let v4 = trr ap cm in
                   result +. v4 *. f
             in
@@ -196,7 +202,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
               if axyz < 1 then result else
                   let f = 0.5 *. (float_of_int axyz) *. expo_inv_q in
                   if abs_float f < cutoff then result else
-                      let am = Powers.decr xyz angMom_a in
+                      let am = Po.decr xyz angMom_a in
                       let v2 = trr am cm in
                       result +. f *. v2
             in
@@ -220,24 +226,24 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
     (** Horizontal recurrence relations *)
     let rec hrr0 angMom_a angMom_b angMom_c =
 
-      match angMom_b.tot with
+      match angMom_b.Po.tot with
       | 1 ->
         let xyz  = get_xyz angMom_b in
-        let ap = Powers.incr xyz angMom_a in
+        let ap = Po.incr xyz angMom_a in
         let v1 = vrr ap angMom_c in
-        let f2 = Coordinate.get xyz center_ab in
+        let f2 = Co.get xyz center_ab in
         if (abs_float f2 < cutoff) then v1 else
             let v2 = vrr angMom_a angMom_c in
             v1 +. f2 *. v2
       | 0 -> vrr angMom_a angMom_c
       | _ ->
         let xyz  = get_xyz angMom_b in
-        let bxyz = Powers.get xyz angMom_b in
+        let bxyz = Po.get xyz angMom_b in
         if bxyz > 0 then 
-            let ap = Powers.incr xyz angMom_a in
-            let bm = Powers.decr xyz angMom_b in
+            let ap = Po.incr xyz angMom_a in
+            let bm = Po.decr xyz angMom_b in
             let h1 = hrr0 ap bm angMom_c in
-            let f2 = Coordinate.get xyz center_ab in
+            let f2 = Co.get xyz center_ab in
             if abs_float f2 < cutoff then h1 else
                 let h2 = hrr0 angMom_a bm angMom_c in
                 h1 +. f2 *. h2
@@ -246,18 +252,18 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
 
     and hrr angMom_a angMom_b angMom_c angMom_d =
 
-      match (angMom_b.tot, angMom_d.tot) with
+      match (angMom_b.Po.tot, angMom_d.Po.tot) with
       | (_,0) ->
-        if (angMom_b.tot = 0) then 
+        if (angMom_b.Po.tot = 0) then 
             vrr angMom_a angMom_c
         else
             hrr0 angMom_a angMom_b angMom_c
       | (_,_) ->
         let xyz  = get_xyz angMom_d in
-        let cp = Powers.incr xyz angMom_c in
-        let dm = Powers.decr xyz angMom_d in 
+        let cp = Po.incr xyz angMom_c in
+        let dm = Po.decr xyz angMom_d in 
         let h1 = hrr angMom_a angMom_b cp dm in
-        let f2 = Coordinate.get xyz center_cd in
+        let f2 = Co.get xyz center_cd in
         if abs_float f2 < cutoff then h1 else
             let h2 = hrr angMom_a angMom_b angMom_c dm in
             h1 +. f2 *. h2
@@ -270,24 +276,20 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
 
 let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
 
-  let shell_a = shell_p.ContractedShellPair.shell_a
-  and shell_b = shell_p.ContractedShellPair.shell_b
-  and shell_c = shell_q.ContractedShellPair.shell_a
-  and shell_d = shell_q.ContractedShellPair.shell_b
-  and sp = shell_p.ContractedShellPair.shell_pairs
-  and sq = shell_q.ContractedShellPair.shell_pairs
-  in
-  let maxm = 
-    shell_p.ContractedShellPair.totAngMomInt +
-    shell_q.ContractedShellPair.totAngMomInt 
+  let shell_a = shell_p.Csp.shell_a
+  and shell_b = shell_p.Csp.shell_b
+  and shell_c = shell_q.Csp.shell_a
+  and shell_d = shell_q.Csp.shell_b
+  and sp      = shell_p.Csp.shell_pairs
+  and sq      = shell_q.Csp.shell_pairs
   in
+  let maxm = shell_p.Csp.totAngMomInt + shell_q.Csp.totAngMomInt in
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
-    Angular_momentum.zkey_array
-      (Angular_momentum.Quartet
-         Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
-                           totAngMom shell_c, totAngMom shell_d))
+    Am.zkey_array (Am.Quartet
+         (shell_a.totAngMom, shell_b.totAngMom,
+          shell_c.totAngMom, shell_d.totAngMom ))
   in
 
   let contracted_class = 
@@ -295,53 +297,53 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   in
 
   let monocentric =
-      shell_p.ContractedShellPair.monocentric &&
-      shell_q.ContractedShellPair.monocentric &&
-      Contracted_shell.center shell_p.ContractedShellPair.shell_a = 
-      Contracted_shell.center shell_q.ContractedShellPair.shell_a
+      shell_p.Csp.monocentric &&
+      shell_q.Csp.monocentric &&
+      shell_p.Csp.shell_a.Cs.center  = 
+      shell_q.Csp.shell_a.Cs.center 
   in
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
-  let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
-  let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in
+  let norm_coef_scale_p = shell_p.Csp.norm_coef_scale in
+  let norm_coef_scale_q = shell_q.Csp.norm_coef_scale in
   for ab=0 to (Array.length sp - 1) do
-    let cab = shell_p.ContractedShellPair.coef.(ab)  in
-    let b = sp.(ab).ShellPair.j in
-    for cd=0 to (Array.length shell_q.ContractedShellPair.shell_pairs - 1) do
+    let cab = shell_p.Csp.coef.(ab)  in
+    let b = sp.(ab).Sp.j in
+    for cd=0 to (Array.length shell_q.Csp.shell_pairs - 1) do
       let coef_prod =
-        cab *. shell_q.ContractedShellPair.coef.(cd)
+        cab *. shell_q.Csp.coef.(cd)
       in
       (** Screening on the product of coefficients *)
       try
-        if (abs_float coef_prod) < 1.e-3*.cutoff then
+        if (abs_float coef_prod) < 1.e-3 *. cutoff then
           raise NullQuartet;
 
         let center_pq = 
-          sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center
+          Co.(sp.(ab).Sp.center |- sq.(cd).Sp.center)
         in
         let norm_pq_sq =
-          Coordinate.dot center_pq center_pq
+          Co.dot center_pq center_pq
         in
 
         let expo_pq_inv =
-          shell_p.ContractedShellPair.expo_inv.(ab) +.
-          shell_q.ContractedShellPair.expo_inv.(cd) 
+          shell_p.Csp.expo_inv.(ab) +.
+          shell_q.Csp.expo_inv.(cd) 
         in
 
         let zero_m_array = 
           zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
         in
         begin
-          match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
-                                  totAngMom shell_c, totAngMom shell_d) with
-          | Angular_momentum.(S,S,S,S) ->
+          match (shell_a.totAngMom, shell_b.totAngMom,
+                 shell_c.totAngMom, shell_d.totAngMom) with
+          | Am.(S,S,S,S) ->
             let integral =
               zero_m_array.(0)
             in
             contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
           | _ -> 
-            let d = shell_q.ContractedShellPair.shell_pairs.(cd).ShellPair.j in
+            let d = shell_q.Csp.shell_pairs.(cd).Sp.j in
             let map_1d = Zmap.create (4*maxm) in
             let map_2d = Zmap.create (Array.length class_indices) in
             let norm_coef_scale =
@@ -362,20 +364,20 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                 try
                   if monocentric then
                     begin
-                      if (  ((1 land angMom_a.x+angMom_b.x+angMom_c.x+angMom_d.x)=1) ||
-                            ((1 land angMom_a.y+angMom_b.y+angMom_c.y+angMom_d.y)=1) ||
-                            ((1 land angMom_a.z+angMom_b.z+angMom_c.z+angMom_d.z)=1) 
+                      if (  ((1 land angMom_a.Po.x + angMom_b.Po.x + angMom_c.Po.x + angMom_d.Po.x)=1) ||
+                            ((1 land angMom_a.Po.y + angMom_b.Po.y + angMom_c.Po.y + angMom_d.Po.y)=1) ||
+                            ((1 land angMom_a.Po.z + angMom_b.Po.z + angMom_c.Po.z + angMom_d.Po.z)=1) 
                           ) then
                         raise NullQuartet
                     end;
 
-                  (*
                   (* Schwartz screening *)
-                  if (maxm > 2) then
+                  (*
+                  if (maxm > 8) then
                   (
                     let schwartz_p = 
-                      let key = Zkey.of_int_tuple (Zkey.Twelve
-                      (angMomA, angMomB, angMomA, angMomB) )
+                      let key =
+                        Zkey.of_powers_twelve angMom_a angMom_b angMom_a angMom_b
                       in
                       match schwartz_p with
                       | None -> 1.
@@ -383,15 +385,15 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                     in
                     if schwartz_p < cutoff then raise NullQuartet;
                     let schwartz_q = 
-                      let key = Zkey.of_int_tuple (Zkey.Twelve
-                      (angMomC, angMomD, angMomC, angMomD) ) 
+                      let key =
+                        Zkey.of_powers_twelve angMom_c angMom_d angMom_c angMom_d
                       in
                       match schwartz_q with
                       | None -> 1.
                       | Some schwartz_q -> Zmap.find schwartz_q key
                     in
                     if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet;
-                    );
+                  );
                   *)
 
 
@@ -399,13 +401,14 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                   let coef_prod = coef_prod *. norm in
 
                   let integral = 
-                    hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
+                    hvrr_two_e
+                      angMom_a angMom_b angMom_c angMom_d
                       zero_m_array
-                      (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
-                      (shell_p.ContractedShellPair.expo_inv.(ab), 
-                       shell_q.ContractedShellPair.expo_inv.(cd) )
-                      (sp.(ab).ShellPair.center_ab, sq.(cd).ShellPair.center_ab, center_pq)
-                      (sp.(ab).ShellPair.center_a , sq.(cd).ShellPair.center_a)
+                      shell_b.expo.(b) shell_d.expo.(d)
+                      shell_p.Csp.expo_inv.(ab)
+                      shell_q.Csp.expo_inv.(cd)
+                      sp.(ab).Sp.center_ab sq.(cd).Sp.center_ab center_pq
+                      sp.(ab).Sp.center_a  sq.(cd).Sp.center_a
                       map_1d map_2d 
                   in
                   contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
@@ -427,8 +430,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 (** Computes all the two-electron integrals of the contracted shell quartet *)
 let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
 
-  let shell_p = ContractedShellPair.create ~cutoff shell_a shell_b
-  and shell_q = ContractedShellPair.create ~cutoff shell_c shell_d
+  let shell_p = Csp.create ~cutoff shell_a shell_b
+  and shell_q = Csp.create ~cutoff shell_c shell_d
   in
   contracted_class_shell_pairs ~zero_m shell_p shell_q 
 
diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 1159627..f3225d0 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -3,6 +3,7 @@ open Lacaml.D
 open Bigarray
 open Powers
 open Coordinate
+open Contracted_shell_type
 
 let cutoff = Constants.cutoff 
 let cutoff2 = cutoff *. cutoff
@@ -565,8 +566,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   let class_indices = 
     Angular_momentum.zkey_array
       (Angular_momentum.Quartet
-         Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
-                           totAngMom shell_c, totAngMom shell_d))
+         (shell_a.totAngMom, shell_b.totAngMom,
+          shell_c.totAngMom, shell_d.totAngMom))
   in
 
   let contracted_class = 
@@ -576,8 +577,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   let monocentric =
       shell_p.ContractedShellPair.monocentric &&
       shell_q.ContractedShellPair.monocentric &&
-      Contracted_shell.center shell_p.ContractedShellPair.shell_a =
-      Contracted_shell.center shell_q.ContractedShellPair.shell_a
+      shell_p.ContractedShellPair.shell_a.center =
+      shell_q.ContractedShellPair.shell_a.center 
   in
 
   (** Screening on the product of coefficients *)
@@ -620,8 +621,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
   begin
-    match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
-                            totAngMom shell_c, totAngMom shell_d) with
+    match (shell_a.totAngMom, shell_b.totAngMom,
+           shell_c.totAngMom, shell_d.totAngMom) with
     | Angular_momentum.(S,S,S,S) ->
       contracted_class.(0) <-
         begin
@@ -681,9 +682,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
       in
 
       let expo_b =
-        Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp
+        Array.map (fun shell_ab -> shell_b.expo.(shell_ab.ShellPair.j)) sp
       and expo_d =
-        Array.map (fun shell_cd -> Contracted_shell.expo shell_d shell_cd.ShellPair.j) sq
+        Array.map (fun shell_cd -> shell_d.expo.(shell_cd.ShellPair.j)) sq
       in
       let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
 
@@ -716,7 +717,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           Array.init np (fun ab ->
             let shell_ab = sp.(ab) in
                 let cpa = 
-                  shell_ab.ShellPair.center |- Contracted_shell.center shell_a
+                  shell_ab.ShellPair.center |- shell_a.center
                 in
                 match xyz with
                 | 0 -> Coordinate.get X cpa;
@@ -735,7 +736,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           Array.init nq (fun cd ->
             let shell_cd = sq.(cd) in
                 let cqc = 
-                  shell_cd.ShellPair.center |- Contracted_shell.center shell_c
+                  shell_cd.ShellPair.center |- shell_c.center
                 in
                 match xyz with
                 | 0 -> Coordinate.get X cqc;
diff --git a/Utils/Util.ml b/Utils/Util.ml
index 7a66a05..5feb872 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -3,8 +3,8 @@ external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxe
 external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed]  [@@noalloc]
 external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed]  [@@noalloc]
 
-open Lacaml.D
 open Constants
+open Lacaml.D
 
 let factmax = 150