diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index b6f490a..110a419 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -121,7 +121,7 @@ let to_file ~filename basis geometry =
     for j_c=1 to (Array2.dim2 eni_array) do
       let value = eni_array.{(i_c-1),(j_c-1)} in
       if (value <> 0.) then
-        Printf.fprintf oc " %5d %5d%20.15f\n" i_c j_c value;
+        Printf.fprintf oc " %5d %5d %20.15f\n" i_c j_c value;
     done;
   done;
   Printf.printf "In: %d  Out:%d\n" !inn !out ;
diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml
index f4e4f66..2be685f 100644
--- a/Basis/OneElectronRR.ml
+++ b/Basis/OneElectronRR.ml
@@ -13,7 +13,7 @@ let chop f g =
 (** Horizontal and Vertical Recurrence Relations (HVRR) *)
 let hvrr_one_e
     (angMom_a, angMom_b) (totAngMom_a, totAngMom_b)
-    (maxm, zero_m_array) (expo_inv_p) (center_ab, center_pa, center_pc)
+    (maxm, zero_m_array) (expo_b) (expo_inv_p) (center_ab, center_pa, center_pc)
     map
   =
 
@@ -45,7 +45,7 @@ let hvrr_one_e
               | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2
             in
             if amxyz < 0 then empty else
-              let f1 = Coordinate.coord center_pa xyz
+              let f1 = Coordinate.coord center_pa xyz 
               and f2 = expo_inv_p *.  (Coordinate.coord center_pc xyz)
               in
               if amxyz < 1 then
@@ -53,7 +53,8 @@ let hvrr_one_e
                   vrr am (totAngMom_a-1)
                 in
                 Array.init maxsze (fun m ->
-                    if m = maxm then 0. else (f1 *. v1.(m) ) -. f2 *. v1.(m+1) )
+                    if m = maxm then (f1 *. v1.(m) ) else
+                    (f1 *. v1.(m) ) -. f2 *. v1.(m+1) )
               else
                 let v3 =
                   vrr amm (totAngMom_a-2)
@@ -62,10 +63,10 @@ let hvrr_one_e
                   vrr am (totAngMom_a-1)
                 in
                 let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
-                Array.init maxsze (fun m ->
+                Array.init maxsze (fun m -> f1 *. v1.(m) -.
                     (if m = maxm then 0. else
-                       (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
-                    +. f3 *. (v3.(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))
                   )
           in Zmap.add map key result;
@@ -148,6 +149,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
 
   for ab=0 to (Array.length shell_p - 1)
   do
+    let b = shell_p.(ab).Shell_pair.j in
     try
       begin
         let coef_prod = shell_p.(ab).Shell_pair.coef in
@@ -165,15 +167,18 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
         let center_ab = 
           shell_p.(ab).Shell_pair.center_ab 
         in
+        let center_p = 
+          shell_p.(ab).Shell_pair.center
+        in
         let center_pa = 
-          Coordinate.(center_ab |- shell_a.Contracted_shell.center) 
+          Coordinate.(center_p |- shell_a.Contracted_shell.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.(shell_p.(ab).Shell_pair.center |- nucl_coord ) 
+            Coordinate.(center_p |- nucl_coord ) 
           in
           let norm_pq_sq = 
             Coordinate.dot center_pc center_pc
@@ -204,7 +209,9 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
                 let integral = 
                   hvrr_one_e (angMomA, angMomB)
                     (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b)
-                    (maxm, zero_m_array) shell_p.(ab).Shell_pair.expo_inv
+                    (maxm, zero_m_array)
+                    (Contracted_shell.expo shell_b b)
+                    (shell_p.(ab).Shell_pair.expo_inv)
                     (center_ab, center_pa, center_pc)
                     map
                 in