diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli
index 2fed7c5..993b17b 100644
--- a/Basis/ContractedShell.mli
+++ b/Basis/ContractedShell.mli
@@ -32,7 +32,7 @@ val to_string : t -> string
 (** Pretty-printing of the contracted shell in a string. *)
 
 val make : ?index:int -> (float * PrimitiveShell.t) array -> t 
-(** Creates a contracted shell from a list of coefficients and primitives. *)
+(** Creates a contracted shell from a list of coefficients and primitives.  *)
 
 val with_index  : t -> int -> t
 (** Returns a copy of the contracted shell with a modified index. *)
diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml
index 2a6a7a8..1005821 100644
--- a/Basis/ContractedShellPair.ml
+++ b/Basis/ContractedShellPair.ml
@@ -59,14 +59,17 @@ Format.printf "@[<2>shell_b :@ %a@]@;" Cs.pp s_b;
   and expo_inv = Array.map (fun (_,y) -> Psp.expo_inv y) shell_pairs
   in
   let shell_pairs = Array.map snd shell_pairs in
-  let root = shell_pairs.(0) in
-  {
-    shell_a = s_a ; shell_b = s_b ; coef ; expo_inv ; shell_pairs ;
-    center_ab = Psp.a_minus_b root;
-    norm_coef_scale = Psp.norm_coef_scale root;
-    norm_sq=Psp.a_minus_b_sq root;
-    totAngMomInt = Psp.totAngMom root |> Am.to_int;
-  }
+  if Array.length shell_pairs = 0 then
+    None
+  else
+    let root = shell_pairs.(0) in
+    Some {
+      shell_a = s_a ; shell_b = s_b ; coef ; expo_inv ; shell_pairs ;
+      center_ab = Psp.a_minus_b root;
+      norm_coef_scale = Psp.norm_coef_scale root;
+      norm_sq=Psp.a_minus_b_sq root;
+      totAngMomInt = Psp.totAngMom root |> Am.to_int;
+    }
 
 
 let shell_a         x = x.shell_a
diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli
index 15faa95..a457789 100644
--- a/Basis/ContractedShellPair.mli
+++ b/Basis/ContractedShellPair.mli
@@ -14,18 +14,22 @@ A contracted shell pair is a product of two {!ContractedShell.t}:
 type t
 
 
-val create : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t
+val create : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option
 (** Creates an contracted shell pair {% $\varphi_{ab}$ %} from a contracted
     shell {% $\chi_a$ %} (first argument) and a contracted shell {% $\chi_b$ %}
     (second argument).
 
-    The contracted shell pair contains the pairs of primitives for which
-    the norm is greater than [cutoff], and for which ... TODO ....
-    All other pairs are discarded.
+    The contracted shell pair contains the only pairs of primitives for which
+    the norm is greater than [cutoff]. 
+
+    If all the primitive shell pairs are not significant, the function returns
+    [None].
 *)
 
-val of_basis : ContractedShell.t array -> t array array
-(** Creates all possible contracted shell pairs from the basis set. *)
+val of_basis : ContractedShell.t array -> t option array array
+(** Creates all possible contracted shell pairs from the basis set.
+    If the shell pair is not significant, sets the value to [None].
+*)
 
 val shell_a : t -> ContractedShell.t
 (** Returns the first {!ContractedShell.t}  {% $\chi_a$ %} which was used to
diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index 0dc0cdd..4a7fca9 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -98,11 +98,13 @@ let of_basis basis =
   let t0 = Unix.gettimeofday () in
 
   let schwartz = 
-    Array.map (fun pair_array -> Array.map (fun pair ->
-        let cls = 
-          contracted_class_shell_pairs pair pair
-        in
-        (cls, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
+    Array.map (fun pair_array -> Array.map (function 
+        | None -> (Zmap.create 0, 0.)
+        | Some pair ->
+          let cls = 
+            contracted_class_shell_pairs pair pair
+          in
+          (cls, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
       ) pair_array ) shell_pairs
   in
 
@@ -138,8 +140,10 @@ let of_basis basis =
       try
         if (schwartz_p_max < cutoff) then raise NullIntegral;
 
-        let
-          shell_p = shell_pairs.(i).(j)
+        let shell_p = 
+          match shell_pairs.(i).(j) with
+          | None -> raise NullIntegral
+          | Some x -> x
         in
 
         let sp =
@@ -152,9 +156,13 @@ let of_basis basis =
             try
             if schwartz_p_max *. schwartz_q_max < cutoff2 then
               raise NullIntegral;
-            let
-              shell_q = shell_pairs.(k).(l)
+
+            let shell_q =
+              match shell_pairs.(k).(l) with
+              | None -> raise NullIntegral
+              | Some x -> x
             in
+
             let sq =
               Csp.shell_pairs shell_q
             in
diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
index 02fdcc2..1b6d6cc 100644
--- a/Basis/KinInt.ml
+++ b/Basis/KinInt.ml
@@ -24,97 +24,98 @@ let to_powers x =
 (** Computes all the kinetic integrals of the contracted shell pair *)
 let contracted_class shell_a shell_b : float Zmap.t = 
 
-  let shell_p =
-    Csp.create shell_a shell_b
-  in
+  match Csp.create shell_a shell_b with
+  | Some shell_p ->
+    begin
+      (* Pre-computation of integral class indices *)
+      let class_indices = 
+        Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b))
+      in
 
-  (* Pre-computation of integral class indices *)
-  let class_indices = 
-    Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b))
-  in
+      let contracted_class = 
+        Array.make (Array.length class_indices) 0.
+      in
 
-  let contracted_class = 
-    Array.make (Array.length class_indices) 0.
-  in
+      (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
-  (* Compute all integrals in the shell for each pair of significant shell pairs *)
+      let sp = Csp.shell_pairs shell_p in
+      let center_ab = 
+        Csp.center_ab shell_p
+      in
+      let norm_coef_scale =
+        Csp.norm_coef_scale shell_p
+      in
 
-  let sp = Csp.shell_pairs shell_p in
-  let center_ab = 
-    Csp.center_ab shell_p
-  in
-  let norm_coef_scale =
-    Csp.norm_coef_scale shell_p
-  in
-
-  for ab=0 to (Array.length sp - 1)
-  do
-    let coef_prod =
-      (Csp.coef shell_p).(ab) 
-    in
-    (** Screening on thr product of coefficients *)
-    if (abs_float coef_prod) > 1.e-4*.cutoff then
-      begin
-        let center_pa = 
-          Psp.center_minus_a sp.(ab)
+      for ab=0 to (Array.length sp - 1)
+      do
+        let coef_prod =
+          (Csp.coef shell_p).(ab) 
         in
-        let expo_inv = 
-          (Csp.expo_inv shell_p).(ab)
-        in
-        let expo_a = 
-          Ps.expo (Psp.shell_a sp.(ab))
-        and expo_b = 
-          Ps.expo (Psp.shell_b sp.(ab))
-        in
-
-        let xyz_of_int k =
-          match k with
-          | 0 -> Co.X
-          | 1 -> Co.Y
-          | _ -> Co.Z
-        in
-        Array.iteri (fun i key ->
-            let (angMomA,angMomB) = to_powers key in
-            let ov a b k = 
-              let xyz = xyz_of_int k in
-              Overlap_primitives.hvrr (a, b)
-                expo_inv
-                (Co.get xyz center_ab,
-                 Co.get xyz center_pa) 
+        (** Screening on thr product of coefficients *)
+        if (abs_float coef_prod) > 1.e-4*.cutoff then
+          begin
+            let center_pa = 
+              Psp.center_minus_a sp.(ab)
             in
-            let f k = 
-               let xyz = xyz_of_int k in 
-               ov (Po.get xyz angMomA) (Po.get xyz angMomB) k
-            and g k =
-               let xyz = xyz_of_int k in 
-               let s1 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB - 1) k
-               and s2 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB - 1) k
-               and s3 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB + 1) k
-               and s4 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB + 1) k
-               and a = float_of_int (Po.get xyz angMomA)
-               and b = float_of_int (Po.get xyz angMomB) 
-               in
-               0.5 *. a *. b *. s1 -. expo_a *. b *. s2 -. expo_b *. a *. s3 +. 
-                 2.0 *. expo_a *. expo_b *. s4
+            let expo_inv = 
+              (Csp.expo_inv shell_p).(ab)
             in
-            let s = Array.init 3 f
-            and k = Array.init 3 g
+            let expo_a = 
+              Ps.expo (Psp.shell_a sp.(ab))
+            and expo_b = 
+              Ps.expo (Psp.shell_b sp.(ab))
             in
-            let norm =  norm_coef_scale.(i) in
-            let integral = chop norm (fun () ->
-                k.(0)*.s.(1)*.s.(2) +.
-                s.(0)*.k.(1)*.s.(2) +.
-                s.(0)*.s.(1)*.k.(2) 
-              ) in 
-            contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
-          ) class_indices
-      end
-  done;
-  let result = 
-    Zmap.create (Array.length contracted_class)
-  in
-  Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
-  result
+
+            let xyz_of_int k =
+              match k with
+              | 0 -> Co.X
+              | 1 -> Co.Y
+              | _ -> Co.Z
+            in
+            Array.iteri (fun i key ->
+                let (angMomA,angMomB) = to_powers key in
+                let ov a b k = 
+                  let xyz = xyz_of_int k in
+                  Overlap_primitives.hvrr (a, b)
+                    expo_inv
+                    (Co.get xyz center_ab,
+                    Co.get xyz center_pa) 
+                in
+                let f k = 
+                  let xyz = xyz_of_int k in 
+                  ov (Po.get xyz angMomA) (Po.get xyz angMomB) k
+                and g k =
+                  let xyz = xyz_of_int k in 
+                  let s1 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB - 1) k
+                  and s2 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB - 1) k
+                  and s3 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB + 1) k
+                  and s4 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB + 1) k
+                  and a = float_of_int (Po.get xyz angMomA)
+                  and b = float_of_int (Po.get xyz angMomB) 
+                  in
+                  0.5 *. a *. b *. s1 -. expo_a *. b *. s2 -. expo_b *. a *. s3 +. 
+                    2.0 *. expo_a *. expo_b *. s4
+                in
+                let s = Array.init 3 f
+                and k = Array.init 3 g
+                in
+                let norm =  norm_coef_scale.(i) in
+                let integral = chop norm (fun () ->
+                    k.(0)*.s.(1)*.s.(2) +.
+                    s.(0)*.k.(1)*.s.(2) +.
+                    s.(0)*.s.(1)*.k.(2) 
+                  ) in 
+                contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+              ) class_indices
+          end
+      done;
+      let result = 
+        Zmap.create (Array.length contracted_class)
+      in
+      Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
+      result
+    end
+  | None -> Zmap.create 0
 
 
 (** Create kinetic energy matrix *)
@@ -149,7 +150,7 @@ let of_basis basis =
               in
               let value =
                 try Zmap.find cls key
-                with Not_found -> failwith "Bug in kinetic integrals"
+                with Not_found -> 0.
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index e86b2f8..b845948 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -64,32 +64,31 @@ let of_basis_nuclei basis nuclei =
   (* Compute Integrals *)
   for i=0 to (Array.length shell) - 1 do
     for j=0 to i do
-      let
-        shell_p = shell_pairs.(i).(j)
-      in
+      match shell_pairs.(i).(j) with
+      | None -> ()
+      | Some shell_p ->
+        (* Compute all the integrals of the class *)
+        let cls =
+          contracted_class_shell_pair shell_p nuclei
+        in
 
-      (* Compute all the integrals of the class *)
-      let cls =
-        contracted_class_shell_pair shell_p nuclei
-      in
-
-      (* Write the data in the output file *)
-      Array.iteri (fun i_c powers_i ->
-          let i_c = Cs.index shell.(i) + i_c + 1 in
-          let xi = to_powers powers_i in
-          Array.iteri (fun j_c powers_j ->
-              let j_c = Cs.index shell.(j) + j_c + 1 in
-              let xj = to_powers powers_j in
-              let key = 
-                Zkey.of_powers_six xi xj
-              in
-              let value = 
-                Zmap.find cls key 
-              in
-              eni_array.{j_c,i_c} <- value;
-              eni_array.{i_c,j_c} <- value;
-            ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j))))
-        ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
+        (* Write the data in the output file *)
+        Array.iteri (fun i_c powers_i ->
+            let i_c = Cs.index shell.(i) + i_c + 1 in
+            let xi = to_powers powers_i in
+            Array.iteri (fun j_c powers_j ->
+                let j_c = Cs.index shell.(j) + j_c + 1 in
+                let xj = to_powers powers_j in
+                let key = 
+                  Zkey.of_powers_six xi xj
+                in
+                let value = 
+                  Zmap.find cls key 
+                in
+                eni_array.{j_c,i_c} <- value;
+                eni_array.{i_c,j_c} <- value;
+              ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j))))
+          ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
     done;
   done;
   Mat.detri eni_array;
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index 2f5fa64..08a1ecb 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -23,9 +23,9 @@ let to_powers x =
 (** Computes all the overlap integrals of the contracted shell pair *)
 let contracted_class shell_a shell_b : float Zmap.t = 
 
-  let shell_p =
-    Csp.create shell_a shell_b
-  in
+  match Csp.create shell_a shell_b with
+  | Some shell_p -> 
+    begin
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
@@ -39,56 +39,59 @@ let contracted_class shell_a shell_b : float Zmap.t =
   let sp =
     Csp.shell_pairs shell_p
   in
-  let center_ab = 
-    Csp.center_ab shell_p
-  in
-  let norm_coef_scale =
-    Csp.norm_coef_scale shell_p
-  in
 
-  (* Compute all integrals in the shell for each pair of significant shell pairs *)
+      let center_ab = 
+        Csp.center_ab shell_p
+      in
+      let norm_coef_scale =
+        Csp.norm_coef_scale shell_p
+      in
 
-  let xyz_of_int k =
-    match k with
-    | 0 -> Co.X
-    | 1 -> Co.Y
-    | _ -> Co.Z
-  in
-  for ab=0 to (Array.length sp - 1)
-  do
-    let coef_prod =
-      (Csp.coef shell_p).(ab) 
-    in
-    (** Screening on thr product of coefficients *)
-    if (abs_float coef_prod) > 1.e-3*.cutoff then
-      begin
-        let expo_inv = 
-          (Csp.expo_inv shell_p).(ab)
+      (* Compute all integrals in the shell for each pair of significant shell pairs *)
+
+      let xyz_of_int k =
+        match k with
+        | 0 -> Co.X
+        | 1 -> Co.Y
+        | _ -> Co.Z
+      in
+      for ab=0 to (Array.length sp - 1)
+      do
+        let coef_prod =
+          (Csp.coef shell_p).(ab) 
         in
-        let center_pa = 
-          Psp.center_minus_a sp.(ab)
-        in
-
-        Array.iteri (fun i key ->
-            let (angMomA,angMomB) = to_powers key in
-            let f k = 
-              let xyz = xyz_of_int k in
-              Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)  
-                expo_inv
-                (Co.get xyz center_ab,
-                 Co.get xyz center_pa) 
+        (** Screening on thr product of coefficients *)
+        if (abs_float coef_prod) > 1.e-3*.cutoff then
+          begin
+            let expo_inv = 
+              (Csp.expo_inv shell_p).(ab)
             in
-            let norm =  norm_coef_scale.(i) in
-            let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in 
-            contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
-          ) class_indices
-      end
-  done;
-  let result = 
-    Zmap.create (Array.length contracted_class)
-  in
-  Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
-  result
+            let center_pa = 
+              Psp.center_minus_a sp.(ab)
+            in
+
+            Array.iteri (fun i key ->
+                let (angMomA,angMomB) = to_powers key in
+                let f k = 
+                  let xyz = xyz_of_int k in
+                  Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)  
+                    expo_inv
+                    (Co.get xyz center_ab,
+                    Co.get xyz center_pa) 
+                in
+                let norm =  norm_coef_scale.(i) in
+                let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in 
+                contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+              ) class_indices
+          end
+      done;
+      let result = 
+        Zmap.create (Array.length contracted_class)
+      in
+      Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
+      result
+    end
+  | None -> Zmap.create 0
 
 
 (** Create overlap matrix *)
@@ -123,7 +126,7 @@ let of_basis basis =
               in
               let value = 
                 try Zmap.find cls key 
-                with Not_found -> failwith "Bug in overlap integrals"
+                with Not_found -> 0.
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml
index fcdf0b0..751d981 100644
--- a/Basis/TwoElectronRR.ml
+++ b/Basis/TwoElectronRR.ml
@@ -432,6 +432,9 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
   let shell_p = Csp.create ~cutoff shell_a shell_b
   and shell_q = Csp.create ~cutoff shell_c shell_d
   in
-  contracted_class_shell_pairs ~zero_m shell_p shell_q 
+  match shell_p, shell_q with
+  | Some shell_p, Some shell_q -> 
+    contracted_class_shell_pairs ~zero_m shell_p shell_q 
+  | _ -> Zmap.create 0
 
 
diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml
index 745e0ae..bd6f9e7 100644
--- a/Basis/TwoElectronRRVectorized.ml
+++ b/Basis/TwoElectronRRVectorized.ml
@@ -873,5 +873,8 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
   let shell_p = Csp.create ~cutoff shell_a shell_b
   and shell_q = Csp.create ~cutoff shell_c shell_d
   in
-  contracted_class_shell_pairs ~zero_m shell_p shell_q 
+  match shell_p, shell_q with
+  | Some shell_p, Some shell_q ->
+    contracted_class_shell_pairs ~zero_m shell_p shell_q
+  | _ -> Zmap.create 0