diff --git a/Basis/Multipole.ml b/Basis/Multipole.ml
index 3b4f38f..cbc9d44 100644
--- a/Basis/Multipole.ml
+++ b/Basis/Multipole.ml
@@ -15,31 +15,15 @@ module Csp = ContractedShellPair
 module Po  = Powers
 module Psp = PrimitiveShellPair
 
-let multiply a b =
-  let n = Mat.dim1 a in
-  let c = Mat.create n n in
-  Mat.cpab c a b;
-  c
-
-let x0 t = t.(0)
-let y0 t = t.(1)
-let z0 t = t.(2)
-let x1 t = t.(3)
-let y1 t = t.(4)
-let z1 t = t.(5)
-let x2 t = t.(6)
-let y2 t = t.(7)
-let z2 t = t.(8)
-
-let matrix_x  t = multiply (x1 t) @@ multiply (y0 t) (z0 t)
-let matrix_y  t = multiply (x0 t) @@ multiply (y1 t) (z0 t)
-let matrix_z  t = multiply (x0 t) @@ multiply (y0 t) (z1 t)
-let matrix_x2 t = multiply (x2 t) @@ multiply (y0 t) (z0 t)
-let matrix_y2 t = multiply (x0 t) @@ multiply (y2 t) (z0 t)
-let matrix_z2 t = multiply (x0 t) @@ multiply (y0 t) (z2 t)
-let matrix_xy t = multiply (x1 t) @@ multiply (y1 t) (z0 t)
-let matrix_yz t = multiply (x0 t) @@ multiply (y1 t) (z1 t)
-let matrix_zx t = multiply (x1 t) @@ multiply (y0 t) (z1 t)
+let matrix_x  t =  t.(0)
+let matrix_y  t =  t.(1)
+let matrix_z  t =  t.(2)
+let matrix_x2 t =  t.(3)
+let matrix_y2 t =  t.(4)
+let matrix_z2 t =  t.(5)
+let matrix_xy t =  t.(6)
+let matrix_yz t =  t.(7)
+let matrix_zx t =  t.(8)
 
 let cutoff = integrals_cutoff
 
@@ -86,9 +70,9 @@ let contracted_class shell_a shell_b : float Zmap.t array =
             begin
               let expo_inv = Psp.exponent_inv psp
               and center_pa = Psp.center_minus_a psp
-              and xa = Co.get Co.X @@ Cs.center shell_a
-              and ya = Co.get Co.Y @@ Cs.center shell_a
-              and za = Co.get Co.Z @@ Cs.center shell_a
+              and xa = Co.(get X) @@ Cs.center shell_a
+              and ya = Co.(get Y) @@ Cs.center shell_a
+              and za = Co.(get Z) @@ Cs.center shell_a
               in
 
               Array.iteri (fun i key ->
@@ -126,15 +110,15 @@ let contracted_class shell_a shell_b : float Zmap.t array =
                   let z2 = h2 +. za *. (2. *. z -. za *. f2) in
                   let c = contracted_class in
                   let d = coef_prod *. norm in
-                  c.(0).(i) <- c.(0).(i) +. d *. f0 ;
-                  c.(1).(i) <- c.(1).(i) +. d *. f1 ;
-                  c.(2).(i) <- c.(2).(i) +. d *. f2 ;
-                  c.(3).(i) <- c.(3).(i) +. d *.  x ;
-                  c.(4).(i) <- c.(4).(i) +. d *.  y ;
-                  c.(5).(i) <- c.(5).(i) +. d *.  z ;
-                  c.(6).(i) <- c.(6).(i) +. d *. x2 ;
-                  c.(7).(i) <- c.(7).(i) +. d *. y2 ;
-                  c.(8).(i) <- c.(8).(i) +. d *. z2 ;
+                  c.(0).(i) <- c.(0).(i) +. d *.  x *. f1 *. f2;
+                  c.(1).(i) <- c.(1).(i) +. d *. f0 *.  y *. f2;
+                  c.(2).(i) <- c.(2).(i) +. d *. f0 *. f1 *.  z;
+                  c.(3).(i) <- c.(3).(i) +. d *. x2 *. f1 *. f2;
+                  c.(4).(i) <- c.(4).(i) +. d *. f0 *. y2 *. f2;
+                  c.(5).(i) <- c.(5).(i) +. d *. f0 *. f1 *. z2;
+                  c.(6).(i) <- c.(6).(i) +. d *.  x *.  y *. f2;
+                  c.(7).(i) <- c.(7).(i) +. d *. f0 *.  y *.  z;
+                  c.(8).(i) <- c.(8).(i) +. d *.  x *. f1 *.  z;
                 ) class_indices
             end
         ) (Csp.coefs_and_shell_pairs shell_p);