diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml
index 4861e15..0de6114 100644
--- a/Basis/TwoElectronRR.ml
+++ b/Basis/TwoElectronRR.ml
@@ -178,7 +178,7 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
   and shell_q = Shell_pair.create_array shell_c shell_d
   and maxm = 
     let open Angular_momentum in
-      (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) 
+    (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) 
     + (to_int @@ Contracted_shell.totAngMom shell_c) + (to_int @@ Contracted_shell.totAngMom shell_d)
   in
 
@@ -186,8 +186,8 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
   let class_indices = 
     Angular_momentum.zkey_array
       (Angular_momentum.Quartet
-        Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
-                          totAngMom shell_c, totAngMom shell_d))
+         Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
+                           totAngMom shell_c, totAngMom shell_d))
   in
 
   let contracted_class = 
@@ -202,51 +202,66 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
 
     for cd=0 to (Array.length shell_q - 1)
     do
-
-      let d = shell_q.(cd).Shell_pair.j in
-
-      let expo_pq_inv =
-        shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv
-      in
-      let center_pq = 
-        Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center)
-      in
-      let norm_pq_sq =
-        Coordinate.dot center_pq center_pq
+      let coef_prod =
+         shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef
       in
+      (** Screening on thr product of coefficients *)
+      if (abs_float coef_prod) > cutoff then
+        begin
 
-      let zero_m_array = 
-        zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
-      in
-      let map = Array.init maxm (fun _ -> Zmap.create 129) in
-      (* Compute the integral class from the primitive shell quartet *)
-      Array.iteri (fun i key ->
-          let (angMomA,angMomB,angMomC,angMomD) =
-            let a = Zkey.to_int_array Zkey.Kind_12 key in
-              ( [| a.(0) ;  a.(1) ;  a.(2) |], 
-                [| a.(3) ;  a.(4) ;  a.(5) |], 
-                [| a.(6) ;  a.(7) ;  a.(8) |], 
-                [| a.(9) ; a.(10) ; a.(11) |] )
-          in
-          let integral = 
-            ghvrr 0 (angMomA, angMomB, angMomC, angMomD)
-              (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,
-               Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)
-              (maxm, zero_m_array)
-              (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
-              (shell_p.(ab).Shell_pair.expo_inv, shell_q.(cd).Shell_pair.expo_inv)
-              (shell_p.(ab).Shell_pair.center_ab, shell_q.(cd).Shell_pair.center_ab, center_pq)
-              map 
-          in
-          let norm = 
-            shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD 
-          in
-          let coef_prod =
-            shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef *. norm
-          in
-          contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
-        ) class_indices
+            let expo_pq_inv =
+              shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv
+            in
+            let center_pq = 
+              Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center)
+            in
+            let norm_pq_sq =
+              Coordinate.dot center_pq center_pq
+            in
 
+            let zero_m_array = 
+              zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
+            in
+
+              match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
+                                      totAngMom shell_c, totAngMom shell_d) with
+              | Angular_momentum.(S,S,S,S) -> Array.iteri (fun i key ->
+                  let coef_prod =
+                    shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef 
+                  in
+                  let integral =
+                    zero_m_array.(0)
+                  in
+                  contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+                ) class_indices
+              | _ ->
+                  let d = shell_q.(cd).Shell_pair.j in
+                  let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in
+                  (* Compute the integral class from the primitive shell quartet *)
+                  Array.iteri (fun i key ->
+                      let (angMomA,angMomB,angMomC,angMomD) =
+                        let a = Zkey.to_int_array Zkey.Kind_12 key in
+                        ( [| a.(0) ;  a.(1) ;  a.(2) |], 
+                          [| a.(3) ;  a.(4) ;  a.(5) |], 
+                          [| a.(6) ;  a.(7) ;  a.(8) |], 
+                          [| a.(9) ; a.(10) ; a.(11) |] )
+                      in
+                      let norm = 
+                        shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD 
+                      in
+                      let integral = chop norm (fun () ->
+                        ghvrr 0 (angMomA, angMomB, angMomC, angMomD)
+                          (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,
+                          Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)
+                          (maxm, zero_m_array)
+                          (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
+                          (shell_p.(ab).Shell_pair.expo_inv, shell_q.(cd).Shell_pair.expo_inv)
+                          (shell_p.(ab).Shell_pair.center_ab, shell_q.(cd).Shell_pair.center_ab, center_pq)
+                          map )
+                      in
+                      contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+                ) class_indices
+        end
     done
   done;
   let result = 
diff --git a/Utils/Constants.ml b/Utils/Constants.ml
index c2d1085..556587a 100644
--- a/Utils/Constants.ml
+++ b/Utils/Constants.ml
@@ -1,4 +1,4 @@
-let cutoff = 1.e-16
+let cutoff = 1.e-15
 
 (** Constants *)
 let pi = acos (-1.)
diff --git a/Utils/Zmap.ml b/Utils/Zmap.ml
index e3c5a85..96696b1 100644
--- a/Utils/Zmap.ml
+++ b/Utils/Zmap.ml
@@ -2,3 +2,4 @@
 
 module Zmap = Hashtbl.Make(Zkey)
 include Zmap
+