diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
new file mode 100644
index 0000000..18dff5c
--- /dev/null
+++ b/Basis/KinInt.ml
@@ -0,0 +1,141 @@
+open Util
+open Constants
+
+(** Computes all the kinetic integrals of the contracted shell pair *)
+let contracted_class shell_a shell_b : float Zmap.t = 
+
+  let shell_p =
+    Shell_pair.create_array 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))
+  in
+
+  let contracted_class = 
+    Array.make (Array.length class_indices) 0.
+  in
+
+  (* Compute all integrals in the shell for each pair of significant shell pairs *)
+
+  for ab=0 to (Array.length shell_p - 1)
+  do
+    let coef_prod =
+      shell_p.(ab).Shell_pair.coef 
+    in
+    (** Screening on thr product of coefficients *)
+    if (abs_float coef_prod) > 1.e-4*.cutoff then
+      begin
+        let center_ab = 
+          shell_p.(ab).Shell_pair.center_ab 
+        in
+        let center_pa = 
+          shell_p.(ab).Shell_pair.center_a
+        in
+        let expo_inv = 
+          shell_p.(ab).Shell_pair.expo_inv
+        in
+        let norm_coef_scale =
+          shell_p.(ab).Shell_pair.norm_coef_scale
+        in
+        let i, j = 
+          shell_p.(ab).Shell_pair.i, shell_p.(ab).Shell_pair.j 
+        in
+        let expo_a = 
+          Contracted_shell.expo shell_p.(ab).Shell_pair.shell_a i
+        and expo_b = 
+          Contracted_shell.expo shell_p.(ab).Shell_pair.shell_b j
+        in
+
+        Array.iteri (fun i key ->
+            let (angMomA,angMomB) =
+              let a = Zkey.to_int_array Zkey.Kind_6 key in
+              ( [| a.(0) ;  a.(1) ;  a.(2) |], 
+                [| a.(3) ;  a.(4) ;  a.(5) |] )
+            in
+            let ov a b k = 
+              Overlap_primitives.hvrr (a, b)
+                expo_inv
+                (Coordinate.coord center_ab k,
+                 Coordinate.coord center_pa k) 
+            in
+            let f k = 
+              ov angMomA.(k) angMomB.(k) k
+            and g k =
+               let s1 = ov (angMomA.(k)-1) (angMomB.(k)-1) k
+               and s2 = ov (angMomA.(k)+1) (angMomB.(k)-1) k
+               and s3 = ov (angMomA.(k)-1) (angMomB.(k)+1) k
+               and s4 = ov (angMomA.(k)+1) (angMomB.(k)+1) k
+               and a = float_of_int angMomA.(k) 
+               and b = float_of_int angMomB.(k)
+               in
+               0.5 *. a *. b *. s1 -. expo_a *. b *. s2 -. expo_b *. a *. s3 +. 
+                 2.0 *. expo_a *. expo_b *. s4
+            in
+            let s = Array.init 3 f
+            and k = Array.init 3 g
+            in
+            let norm =  norm_coef_scale.(i) in
+            let integral = chop norm (fun () ->
+                k.(0)*.s.(1)*.s.(2) +.
+                s.(0)*.k.(1)*.s.(2) +.
+                s.(0)*.s.(1)*.k.(2) 
+              ) in 
+            contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+          ) class_indices
+      end
+  done;
+  let result = 
+    Zmap.create (Array.length contracted_class)
+  in
+  Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
+  result
+
+
+
+
+
+(** Write all kinetic integrals to a file *)
+let to_file ~filename basis =
+  let to_int_tuple x = 
+    let open Zkey in
+    match to_int_tuple Kind_3 x with
+    | Three x -> x
+    | _ -> assert false
+  in
+
+  let oc = open_out filename in
+  for i=0 to (Array.length basis) - 1 do
+    print_int basis.(i).Contracted_shell.indice ; print_newline ();
+    for j=0 to i do
+      (* Compute all the integrals of the class *)
+      let cls =
+        contracted_class basis.(i) basis.(j) 
+      in
+
+      (* 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 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 xj = to_int_tuple powers_j in
+              let key = 
+                Zkey.of_int_tuple (Zkey.Six (xi,xj))
+              in
+              let value = 
+                try
+                  Zmap.find cls key 
+                with Not_found -> failwith "Bug in Kinetic integrals"
+              in
+              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;
+    done;
+  done;
+  close_out oc
+
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index e3d7da7..4c31890 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -5,13 +5,13 @@ open Constants
 let contracted_class shell_a shell_b : float Zmap.t = 
 
   let shell_p =
-     Shell_pair.create_array shell_a shell_b
+    Shell_pair.create_array 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))
+                                   Contracted_shell.(totAngMom shell_a, totAngMom shell_b))
   in
 
   let contracted_class = 
@@ -22,42 +22,42 @@ let contracted_class shell_a shell_b : float Zmap.t =
 
   for ab=0 to (Array.length shell_p - 1)
   do
-      let coef_prod =
-         shell_p.(ab).Shell_pair.coef 
-      in
-      (** Screening on thr product of coefficients *)
-      if (abs_float coef_prod) > 1.e-4*.cutoff then
-        begin
-            let center_ab = 
-              shell_p.(ab).Shell_pair.center_ab 
-            in
-            let center_pa = 
-              shell_p.(ab).Shell_pair.center_a
-            in
-            let expo_inv = 
-              shell_p.(ab).Shell_pair.expo_inv
-            in
-            let norm_coef_scale =
-	      shell_p.(ab).Shell_pair.norm_coef_scale
-	    in
+    let coef_prod =
+      shell_p.(ab).Shell_pair.coef 
+    in
+    (** Screening on thr product of coefficients *)
+    if (abs_float coef_prod) > 1.e-4*.cutoff then
+      begin
+        let center_ab = 
+          shell_p.(ab).Shell_pair.center_ab 
+        in
+        let center_pa = 
+          shell_p.(ab).Shell_pair.center_a
+        in
+        let expo_inv = 
+          shell_p.(ab).Shell_pair.expo_inv
+        in
+        let norm_coef_scale =
+          shell_p.(ab).Shell_pair.norm_coef_scale
+        in
 
-            Array.iteri (fun i key ->
-                      let (angMomA,angMomB) =
-                        let a = Zkey.to_int_array Zkey.Kind_6 key in
-                        ( [| a.(0) ;  a.(1) ;  a.(2) |], 
-                          [| a.(3) ;  a.(4) ;  a.(5) |] )
-                      in
-                      let f k = 
-                        Overlap_primitives.hvrr (angMomA.(k), angMomB.(k)) 
-                        expo_inv
-                        (Coordinate.coord center_ab k,
-                         Coordinate.coord center_pa k) 
-                      in
-		      let norm =  norm_coef_scale.(i) in
-                      let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in 
-                      contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
-                ) class_indices
-        end
+        Array.iteri (fun i key ->
+            let (angMomA,angMomB) =
+              let a = Zkey.to_int_array Zkey.Kind_6 key in
+              ( [| a.(0) ;  a.(1) ;  a.(2) |], 
+                [| a.(3) ;  a.(4) ;  a.(5) |] )
+            in
+            let f k = 
+              Overlap_primitives.hvrr (angMomA.(k), angMomB.(k)) 
+                expo_inv
+                (Coordinate.coord center_ab k,
+                 Coordinate.coord center_pa k) 
+            in
+            let norm =  norm_coef_scale.(i) in
+            let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in 
+            contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+          ) class_indices
+      end
   done;
   let result = 
     Zmap.create (Array.length contracted_class)
diff --git a/README.rst b/README.rst
new file mode 100644
index 0000000..09f23ea
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,10 @@
+Requirements
+------------
+
+* gsl: GNU Scientific Library
+* gmp : GNU Multiple Precision arithmetic library
+* zarith : Arbitrary-precision integers
+* BLAS/LAPACK : Linear algebra
+* LaCaml : LAPACK OCaml interface
+* SklMl : Parallel skeletons for OCaml
+
diff --git a/run_integrals.ml b/run_integrals.ml
index 7c5aa87..f7bdc6e 100644
--- a/run_integrals.ml
+++ b/run_integrals.ml
@@ -29,6 +29,7 @@ let run ~out =
 
   Overlap.to_file ~filename:(out_file^".overlap") basis;
   NucInt.to_file ~filename:(out_file^".nuc") basis nuclei;
+  KinInt.to_file ~filename:(out_file^".kin") basis;
   ERI.to_file ~filename:(out_file^".eri") basis