diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml
index 82499bb..ca71c9e 100644
--- a/Basis/TwoElectronRR.ml
+++ b/Basis/TwoElectronRR.ml
@@ -17,26 +17,42 @@ let cutoff2 = cutoff *. cutoff
 
 exception NullQuartet
 
+type four_idx_intermediates =
+{
+  expo_b : float ;
+  expo_d : float ;
+  expo_inv_p : float ;
+  expo_inv_q : float ;
+  center_ab  : Co.t ;
+  center_cd  : Co.t ;
+  center_pq  : Co.t ;
+  center_pa  : Co.t ;
+  center_qc  : Co.t ;
+  zero_m_array : float array ;
+}
+
 (** Horizontal and Vertical Recurrence Relations (HVRR) *)
 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 =
+    abcd map_1d map_2d =
 
   (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *)
   if  angMom_a.Po.tot + angMom_b.Po.tot < angMom_c.Po.tot + angMom_d.Po.tot then
+    let abcd = {
+        expo_b       = abcd.expo_d ;
+        expo_d       = abcd.expo_b ;
+        expo_inv_p   = abcd.expo_inv_q ;
+        expo_inv_q   = abcd.expo_inv_p ;
+        center_ab    = abcd.center_cd ;
+        center_cd    = abcd.center_ab ;
+        center_pq    = Co.neg abcd.center_pq ;
+        center_pa    = abcd.center_qc ;
+        center_qc    = abcd.center_pa ;
+        zero_m_array = abcd.zero_m_array ;
+    } in
     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 (Co.neg center_pq) 
-      center_qc center_pa
-      map_1d map_2d
+      abcd map_1d map_2d
 
   else
 
@@ -51,12 +67,18 @@ let rec hvrr_two_e
       | _                    -> Co.Z
     in
 
+    let expo_inv_p = abcd.expo_inv_p
+    and expo_inv_q = abcd.expo_inv_q
+    and center_ab  = abcd.center_ab
+    and center_cd  = abcd.center_cd
+    and center_pq  = abcd.center_pq
+    in
 
     (** Vertical recurrence relations *)
     let rec vrr0 angMom_a = 
 
       match angMom_a.Po.tot with 
-      | 0 -> zero_m_array
+      | 0 -> abcd.zero_m_array
       | _ ->
         let key = Zkey.of_powers_three angMom_a in
 
@@ -68,7 +90,7 @@ let rec hvrr_two_e
             let amxyz = Po.get xyz am in
 
             let f1 = expo_inv_p *. Co.get xyz center_pq
-            and f2 = expo_b *. expo_inv_p *. Co.get xyz center_ab
+            and f2 = abcd.expo_b *. expo_inv_p *. Co.get xyz center_ab
             in
             let result = Array.create_float (maxsze - angMom_a.Po.tot) in
             if amxyz = 0 then
@@ -96,7 +118,7 @@ let rec hvrr_two_e
 
       match angMom_a.Po.tot, angMom_c.Po.tot with
       | (i,0) -> if (i>0) then vrr0 angMom_a
-                 else zero_m_array 
+                 else abcd.zero_m_array 
       | (_,_) ->
         let key = Zkey.of_powers_six angMom_a angMom_c in
 
@@ -110,7 +132,7 @@ let rec hvrr_two_e
             let axyz  = Po.get xyz angMom_a in
 
             let f1 =
-              -. expo_d *. expo_inv_q *. Co.get xyz center_cd
+              -. abcd.expo_d *. expo_inv_q *. Co.get xyz center_cd
             and f2 = 
               expo_inv_q *. Co.get xyz center_pq
             in
@@ -162,7 +184,7 @@ let rec hvrr_two_e
 
       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)
+                 else abcd.zero_m_array.(0)
       | (_,_) ->
         let key = Zkey.of_powers_six angMom_a angMom_c in
 
@@ -359,15 +381,15 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t
                   let norm =  norm_scales.(i) in
                   let coef_prod = coef_prod *. norm in
 
+                  let abcd = {
+                    expo_b ; expo_d ; expo_inv_p ; expo_inv_q ;
+                    center_ab ; center_cd ; center_pq ;
+                    center_pa ; center_qc ; zero_m_array ;
+                  } in
                   let integral = 
                     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 
+                      abcd map_1d map_2d 
                   in
                   contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
                 with NullQuartet -> ()
@@ -471,14 +493,15 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple :
                   let norm =  norm_scales.(i) in
                   let coef_prod = coef_prod *. norm in
 
+                  let abcd = {
+                    expo_b ; expo_d ; expo_inv_p ; expo_inv_q ;
+                    center_ab ; center_cd ; center_pq ;
+                    center_pa ; center_qc ; zero_m_array ;
+                  } in
                   let integral = 
                     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 
+                      abcd
                       map_1d map_2d 
                   in
                   contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral