diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 4babfc8..d148bc0 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -53,31 +53,32 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d)
                         in
                         am.(xyz)  <- am.(xyz) - 1;
                         amm.(xyz) <- amm.(xyz) - 2;
-                        if am.(xyz) < 0 then
-                          empty 
-                        else
+                        if am.(xyz) < 0 then empty else
                           let v1 = 
                             let f = 
                               -. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)
                             in
-                            Array.mapi (fun k v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) )
+                            if (abs_float f < cutoff) then empty else
+                              Array.mapi (fun k v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) )
                           in
                           let p1 = 
                                
                             Array.mapi (fun k v2k -> v1.(k) +. expo_inv_p *. (Coordinate.coord center_pq.(k) xyz) *. v2k) (vrr0_v (m+1) am (totAngMom_a-1))
                           in
                           if amm.(xyz) < 0 then p1 else
-                            let v1 = vrr0_v m amm (totAngMom_a-2)  
-                            and v2 = vrr0_v (m+1) amm (totAngMom_a-2) 
-                            and f = (float_of_int am.(xyz)) *. expo_inv_p *. 0.5
+                            let f = (float_of_int am.(xyz)) *. expo_inv_p *. 0.5
                             in
-                            Array.mapi (fun k _ -> p1.(k) +.
-                                f *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) coef_prod
-      )
-      in
-      if not found then
-        Zmap.add map.(m) key result;
-      result
+                            if (abs_float f < cutoff) then empty else
+                              let v1 = vrr0_v m amm (totAngMom_a-2)  
+                              and v2 = 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
+                          )
+                          in
+                          if not found then
+                            Zmap.add map.(m) key result;
+                          result
 
   and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c =
 
@@ -135,14 +136,18 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d)
                                 ) coef_prod
                             in
                             if (am.(xyz) < 0) || (cm.(xyz) < 0) then p2 else
-                              let v = 
-                                vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) 
-                              and fa =
+                              let fa =
                                 (float_of_int angMom_a.(xyz)) *. expo_inv_p *. 0.5
                               in
-                              Array.mapi (fun k _ -> 
-                                  p2.(k) -. fa *. expo_inv_q.(k) *. v.(k) 
-                              ) coef_prod
+                              (*
+                              if (abs_float fa < cutoff) then empty else
+                              *)
+                                let v = 
+                                  vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) 
+                                in
+                                Array.mapi (fun k _ -> 
+                                    p2.(k) -. fa *. expo_inv_q.(k) *. v.(k) 
+                                ) coef_prod
       )
       in
       if not found then
@@ -167,13 +172,15 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d)
     | 1 -> let xyz = if angMom_b.(0) = 1 then 0 else if angMom_b.(1) = 1 then 1 else 2 in
       let ap  = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in
       ap.(xyz) <- ap.(xyz) + 1;
+      let f = Coordinate.coord center_ab xyz in
       let v1 = 
-        vrr_v m ap angMom_c (totAngMom_a+1) totAngMom_c
-      and v2 = 
-        vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c
-      and f = Coordinate.coord center_ab xyz
+          vrr_v m ap angMom_c (totAngMom_a+1) totAngMom_c
       in
-      Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2
+      if (abs_float f < cutoff) then v1 else
+        let v2 = 
+          vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c
+        in
+        Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2
     | _ ->
       let ap  = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
       and bm  = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |]
@@ -188,9 +195,11 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d)
       if (bm.(xyz) < 0) then empty  else
         let h1 = 
           hrr0_v m ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c 
-        and h2 = 
+        in
+        let f = (Coordinate.coord center_ab xyz) in
+        if (abs_float f < cutoff) then h1 else
+        let h2 = 
           hrr0_v m angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c 
-        and f = (Coordinate.coord center_ab xyz)
         in Array.map2 (fun h1 h2 -> h1 +. h2 *. f) h1 h2
 
   and hrr_v m angMom_a angMom_b angMom_c angMom_d