diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 426fd8a..e555d14 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -458,12 +458,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
           hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1)
         in h1 +. f *. h2
   in
-  hrr_v 
-      (angMom_a.(0),angMom_a.(1),angMom_a.(2))
-      (angMom_b.(0),angMom_b.(1),angMom_b.(2))
-      (angMom_c.(0),angMom_c.(1),angMom_c.(2))
-      (angMom_d.(0),angMom_d.(1),angMom_d.(2))
-      totAngMom_a totAngMom_b totAngMom_c totAngMom_d 
+  hrr_v angMom_a angMom_b angMom_c angMom_d
+    totAngMom_a totAngMom_b totAngMom_c totAngMom_d 
 
 
 
@@ -497,6 +493,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
     Array.make (Array.length class_indices) 0.;
   in
 
+  let monocentric =
+      shell_p.ContractedShellPair.monocentric &&
+      shell_q.ContractedShellPair.monocentric
+  in
+
   (** Screening on the product of coefficients *)
   let coef_max_p =
     Array.fold_left (fun accu x ->
@@ -671,25 +672,47 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
         in
           (* Compute the integral class from the primitive shell quartet *)
             Array.iteri (fun i key ->
-                let a = Zkey.to_int_array Zkey.Kind_12 key in
                 let (angMomA,angMomB,angMomC,angMomD) =
-                  ( [| a.(0) ;  a.(1) ;  a.(2) |], 
-                    [| a.(3) ;  a.(4) ;  a.(5) |], 
-                    [| a.(6) ;  a.(7) ;  a.(8) |], 
-                    [| a.(9) ; a.(10) ; a.(11) |] )
+                  match Zkey.to_int_tuple ~kind:Zkey.Kind_12 key with
+                  | Zkey.Twelve x -> x
+                  | _             -> assert false
                 in
-                let integral = 
-                  hvrr_two_e_vector (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)
-                    (expo_b, expo_d)
-                    (expo_inv_p, expo_inv_q)
-                    (shell_p.ContractedShellPair.center_ab,
-                     shell_q.ContractedShellPair.center_ab, center_pq)
-                     map_1d map_2d np nq
-                in
-                contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i)
+                try
+                    if monocentric then
+                      begin
+                        let ax,ay,az = angMomA
+                        and bx,by,bz = angMomB
+                        and cx,cy,cz = angMomC
+                        and dx,dy,dz = angMomD
+                        in
+                        if ( ((1 land ax+bx+cx+dx)=1) ||
+                             ((1 land ay+by+cy+dy)=1) ||
+                             ((1 land az+bz+cz+dz)=1)
+                           ) then
+                          raise NullQuartet
+                      end;
+  (*
+                  let a = Zkey.to_int_array Zkey.Kind_12 key in
+                  let (angMomA,angMomB,angMomC,angMomD) =
+                    ( [| 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 = 
+                    hvrr_two_e_vector (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)
+                      (expo_b, expo_d)
+                      (expo_inv_p, expo_inv_q)
+                      (shell_p.ContractedShellPair.center_ab,
+                      shell_q.ContractedShellPair.center_ab, center_pq)
+                      map_1d map_2d np nq
+                  in
+                  contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i)
+                with NullQuartet -> ()
               ) class_indices
 
   end;