diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml
index 6645c24..fd6710b 100644
--- a/Basis/TwoElectronRR.ml
+++ b/Basis/TwoElectronRR.ml
@@ -305,7 +305,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   let norm_coef_scale_q = Csp.norm_coef_scale shell_q in
   for ab=0 to (Array.length sp - 1) do
     let cab = (Csp.coef shell_p).(ab)  in
-    let b = sp.(ab).Sp.j in
+    let b = (Csp.shell_pairs shell_p).(ab).Sp.j in
+    let expo_b = (Cs.expo shell_b).(b) in
+    let expo_inv_p = (Csp.expo_inv shell_p).(ab) in
+    let center_ab = sp.(ab).Sp.center_ab in
+    let center_pa = sp.(ab).Sp.center_a in
     for cd=0 to (Array.length (Csp.shell_pairs shell_q) - 1) do
       let coef_prod =
         cab *. (Csp.coef shell_q).(cd)
@@ -322,10 +326,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           Co.dot center_pq center_pq
         in
 
-        let expo_pq_inv =
-         (Csp.expo_inv shell_p).(ab) +.
-         (Csp.expo_inv shell_q).(cd) 
-        in
+        let expo_inv_q = (Csp.expo_inv shell_q).(cd) in
+
+        let expo_pq_inv = expo_inv_p +. expo_inv_q in
 
         let zero_m_array = 
           zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
@@ -340,8 +343,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
             contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
           | _ -> 
             let d = (Csp.shell_pairs shell_q).(cd).Sp.j in
+            let expo_d = (Cs.expo shell_d).(d) in
             let map_1d = Zmap.create (4*maxm) in
             let map_2d = Zmap.create (Array.length class_indices) in
+            let center_cd = sq.(cd).Sp.center_ab in
+            let center_qc = sq.(cd).Sp.center_a in
             let norm_coef_scale =
               Array.to_list norm_coef_scale_p
               |> List.map (fun v1 ->
@@ -400,11 +406,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                     hvrr_two_e
                       angMom_a angMom_b angMom_c angMom_d
                       zero_m_array
-                      (Cs.expo shell_b).(b) (Cs.expo shell_d).(d)
-                      (Csp.expo_inv shell_p).(ab)
-                      (Csp.expo_inv shell_q).(cd)
-                      sp.(ab).Sp.center_ab sq.(cd).Sp.center_ab center_pq
-                      sp.(ab).Sp.center_a  sq.(cd).Sp.center_a
+                      expo_b  expo_d
+                      expo_inv_p expo_inv_q
+                      center_ab center_cd center_pq
+                      center_pa center_qc 
                       map_1d map_2d 
                   in
                   contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 816f3e4..3a7450e 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -777,9 +777,12 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
               ) sq
             in
             (* Transpose result *)
+              let coef_ab = coef.(ab) in
               for m=0 to maxm do
+                let result_m_ab = result.(m).(ab)
+                in
                 for cd=0 to nq-1 do
-                  result.(m).(ab).(cd) <- zero_m_array_tmp.(cd).(m) *. coef.(ab).(cd)
+                  result_m_ab.(cd) <- zero_m_array_tmp.(cd).(m) *. coef_ab.(cd)
                 done
               done
           ) sp;