diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index 1edacad..23a2fb1 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -13,6 +13,7 @@ open Bigarray
     norm_pq_sq  : square of the distance between the centers of $\phi_p$
                 and $\phi_q$
 *)
+
 let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq =
   let exp_pq = 1. /. expo_pq_inv in
   let t = norm_pq_sq *. exp_pq in
diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml
index 7ad24c5..5a564d8 100644
--- a/Basis/TwoElectronRR.ml
+++ b/Basis/TwoElectronRR.ml
@@ -325,6 +325,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
+  let zero_m_cache =
+    Hashtbl.create 129
+  in
   for ab=0 to (Array.length shell_p - 1) do
     let cab = shell_p.(ab).Shell_pair.coef  in
     let b = shell_p.(ab).Shell_pair.j in
@@ -349,7 +352,20 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
         in
 
         let zero_m_array = 
-          zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
+          let key = String.concat " " 
+            [ string_of_int maxm ; string_of_float expo_pq_inv ; string_of_float norm_pq_sq ]
+          in
+          try
+            let result = 
+              Hashtbl.find zero_m_cache key
+            in
+            result
+          with
+          | Not_found ->
+            let result = 
+              zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
+            in
+            (Hashtbl.add zero_m_cache key result ; result)
         in
         begin
           match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 91b7a76..0ad0df0 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -306,6 +306,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
     Array.make (Array.length class_indices) 0.;
   in
 
+  let zero_m_cache =
+     Hashtbl.create 129
+  in
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
   begin
@@ -334,8 +337,21 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                   Coordinate.dot center_pq center_pq
                 in
 
-                let zero_m_array =
-                  zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq
+                let zero_m_array = 
+                  let key = String.concat " " 
+                    [ string_of_int 0 ; string_of_float expo_pq_inv ; string_of_float norm_pq_sq ]
+                  in
+                  try
+                    let result = 
+                      Hashtbl.find zero_m_cache key
+                    in
+                    result
+                  with
+                  | Not_found ->
+                    let result = 
+                      zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq
+                    in
+                    (Hashtbl.add zero_m_cache key result ; result)
                 in
 
                 accu +. coef_prod *. zero_m_array.(0)
@@ -362,8 +378,21 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                   Coordinate.dot center_pq center_pq
                 in
 
-                let zero_m_array =
-                  zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
+                let zero_m_array = 
+                  let key = String.concat " " 
+                    [ string_of_int maxm ; string_of_float expo_pq_inv ; string_of_float norm_pq_sq ]
+                  in
+                  try
+                    let result = 
+                      Hashtbl.find zero_m_cache key
+                    in
+                    result
+                  with
+                  | Not_found ->
+                    let result = 
+                      zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
+                    in
+                    (Hashtbl.add zero_m_cache key result ; result)
                 in
 
                 let d = shell_cd.Shell_pair.j in