diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 6123a8c..917b025 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -42,11 +42,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
         | _       -> 2 
       in
       let f = expo_b *. (Coordinate.coord center_ab xyz) in
-      Array.mapi (fun k c -> c *.  expo_inv_p *.
+      Array.init ncoef (fun k -> coef_prod.(k) *.  expo_inv_p *.
         ( (Coordinate.coord center_pq.(k) xyz) *. zero_m_array.(m+1).(k)
-          -. f *. zero_m_array.(m).(k) ) ) coef_prod
-
-    | 0 -> Array.mapi (fun k c -> c *. zero_m_array.(m).(k)) coef_prod
+          -. f *. zero_m_array.(m).(k) ) ) 
+    | 0 -> Array.map2 ( *. ) zero_m_array.(m) coef_prod
     | totAngMom_a ->
       let key = Zkey.of_int_tuple (Zkey.Three angMom_a)
       in
@@ -66,7 +65,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 		    -. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)
 		  in
 		  if (abs_float f < cutoff) then empty else
-		    Array.mapi (fun k v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) )
+		    Array.map (fun v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) )
 		in
 		let p1 = 
 		      
@@ -82,8 +81,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 		      if (abs_float (f *. expo_inv_p)) < cutoff then empty else
 		      vrr0_v (m+1) amm (totAngMom_a-2) 
 		    in
-		    Array.mapi (fun k _ -> p1.(k) +.
-			f *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) coef_prod
+		    Array.init ncoef (fun k -> p1.(k) +.
+			f *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) 
 		in Zmap.add map_1d.(m) key result;
 	    result
 
@@ -93,7 +92,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
     | (i,0) -> if (i>0) then
                  vrr0_v m angMom_a totAngMom_a
                else
-                 Array.mapi (fun k c -> c *. zero_m_array.(m).(k)) coef_prod
+                 Array.map2 ( *. ) zero_m_array.(m) coef_prod
     | (_,_) ->
 
       let key =  Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c))
@@ -125,12 +124,13 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 		in
 		if cxyz < 1 then empty else
 		  let f1 =
-		    Array.mapi (fun k _ -> 
-			expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) ) expo_inv_q
+		    Array.init ncoef (fun k -> 
+			expo_d.(k) *. expo_inv_q.(k) *.
+                        (Coordinate.coord center_cd.(k) xyz) ) 
 		  in
 		  let f2 =
-		    Array.mapi (fun k _ -> 
-			expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz) ) expo_inv_q
+		    Array.init ncoef (fun k -> 
+			expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz) )
 		  in
 		  let v1 = 
 		    if (at_least_one_valid f1) then
@@ -142,7 +142,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 		    else empty
 		  in
 		  let p1 = 
-		    Array.mapi (fun k _ -> -. v1.(k) *. f1.(k) -. v2.(k) *. f2.(k)) coef_prod
+		    Array.init ncoef (fun k -> -. v1.(k) *. f1.(k) -. v2.(k) *. f2.(k)) 
 		  in
 		  let p2 = 
 		    if cxyz < 2 then p1 else
@@ -150,10 +150,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 			(float_of_int (cxyz-1)) *. 0.5
 		      in
 		      let f1 =
-			Array.mapi (fun k _ -> fcm *. expo_inv_q.(k) ) coef_prod
+			Array.map (fun e -> fcm *. e) expo_inv_q
 		      in
 		      let f2 =
-			Array.mapi (fun k _ -> f1.(k) *. expo_inv_q.(k) ) coef_prod
+			Array.map2 ( *. ) f1 expo_inv_q 
 		      in
 		      let v1 = 
 			if (at_least_one_valid f1) then
@@ -165,20 +165,20 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 			  vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) 
 			else empty
 		      in
-		      Array.mapi (fun k _ -> p1.(k) +.  f1.(k) *. v1.(k) +. f2.(k) *. v2.(k)) coef_prod
+		      Array.init ncoef (fun k -> p1.(k) +.  f1.(k) *. v1.(k) +. f2.(k) *. v2.(k)) 
 		    in
 		    if (axyz < 1) || (cxyz < 1) then p2 else
 		      let fa =
 			(float_of_int axyz) *. expo_inv_p *. 0.5
 		      in
 		      let f1 =
-			Array.mapi (fun k _ -> fa *. expo_inv_q.(k) ) coef_prod
+			Array.map (fun e -> fa *. e ) expo_inv_q
 		      in
 		      if (at_least_one_valid f1) then
 			let v = 
 			  vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) 
 			in
-			Array.mapi (fun k _ -> p2.(k) -. f1.(k) *. v.(k)) coef_prod
+			Array.init ncoef (fun k -> p2.(k) -. f1.(k) *. v.(k)) 
 		      else p2
               end
 	      in Zmap.add map_2d.(m) key result;
@@ -195,7 +195,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
     | 0 ->
       begin
         match (totAngMom_a, totAngMom_c) with
-        | (0,0) -> Array.mapi (fun k c -> c *. zero_m_array.(0).(k)) coef_prod
+        | (0,0) -> Array.map2 ( *. )  zero_m_array.(0) coef_prod
         | (_,0) -> vrr0_v 0 angMom_a totAngMom_a
         | (_,_) -> vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c
       end