diff --git a/Basis/AtomicShellPairCouple.ml b/Basis/AtomicShellPairCouple.ml
new file mode 100644
index 0000000..6d90bae
--- /dev/null
+++ b/Basis/AtomicShellPairCouple.ml
@@ -0,0 +1,72 @@
+type t = 
+{
+  atomic_shell_pair_p: AtomicShellPair.t ;
+  atomic_shell_pair_q: AtomicShellPair.t ;
+  atomic_shell_a     : AtomicShell.t ;
+  atomic_shell_b     : AtomicShell.t ;
+  atomic_shell_c     : AtomicShell.t ;
+  atomic_shell_d     : AtomicShell.t ;
+  ang_mom     : AngularMomentum.t ;
+  contracted_shell_pair_couples : ContractedShellPairCouple.t list ;
+}
+
+module Am   = AngularMomentum
+module Co   = Coordinate
+module As   = AtomicShell
+module Asp  = AtomicShellPair
+module Cspc = ContractedShellPairCouple
+
+let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q = 
+  let ang_mom =
+    Am.(Asp.ang_mom atomic_shell_pair_p + Asp.ang_mom atomic_shell_pair_q)
+  in
+  let atomic_shell_a = Asp.atomic_shell_a atomic_shell_pair_p
+  and atomic_shell_b = Asp.atomic_shell_b atomic_shell_pair_p
+  and atomic_shell_c = Asp.atomic_shell_a atomic_shell_pair_q
+  and atomic_shell_d = Asp.atomic_shell_b atomic_shell_pair_q
+  in
+  let contracted_shell_pair_couples =
+    List.map (fun ap_ab -> 
+      List.map (fun ap_cd -> 
+        ContractedShellPairCouple.make ~cutoff ap_ab ap_cd
+      ) (Asp.contracted_shell_pairs atomic_shell_pair_q)
+    ) (Asp.contracted_shell_pairs atomic_shell_pair_p)
+    |> List.concat
+    |> List.filter (function None -> false | _ -> true)
+    |> List.map (function None -> assert false | Some x -> x)
+  in
+  match contracted_shell_pair_couples with
+  | [] -> None 
+  | _  -> Some  { atomic_shell_pair_p ; atomic_shell_pair_q ; ang_mom ;
+                  atomic_shell_a ; atomic_shell_b ; atomic_shell_c ; atomic_shell_d ;
+                  contracted_shell_pair_couples ;
+                }
+
+let contracted_shell_pair_couples t = t.contracted_shell_pair_couples
+
+let monocentric t =
+  Asp.monocentric t.atomic_shell_pair_p && Asp.monocentric t.atomic_shell_pair_q &&
+  As.center (Asp.atomic_shell_a t.atomic_shell_pair_p) = As.center (Asp.atomic_shell_a t.atomic_shell_pair_q)
+
+
+let ang_mom t = t.ang_mom
+
+let atomic_shell_pair_p t = t.atomic_shell_pair_p
+let atomic_shell_pair_q t = t.atomic_shell_pair_q
+
+let atomic_shell_a t = t.atomic_shell_a
+let atomic_shell_b t = t.atomic_shell_b
+let atomic_shell_c t = t.atomic_shell_c
+let atomic_shell_d t = t.atomic_shell_d
+
+
+let zkey_array t =
+  match t.contracted_shell_pair_couples with
+  | f::_ -> Cspc.zkey_array f
+  |  _ -> invalid_arg "AtomicShellPairCouple.zkey_array"
+
+let norm_scales t = 
+  match t.contracted_shell_pair_couples with
+  | f::_ -> Cspc.norm_scales f
+  |  _ -> invalid_arg "AtomicShellPairCouple.norm_scales"
+
diff --git a/Basis/AtomicShellPairCouple.mli b/Basis/AtomicShellPairCouple.mli
new file mode 100644
index 0000000..84ec73c
--- /dev/null
+++ b/Basis/AtomicShellPairCouple.mli
@@ -0,0 +1,61 @@
+(** Data structure describing a couple of atomic shells pairs.
+
+An atomic shell pair couple is the cartesian product between two sets of functions, one
+set over electron one and one set over electron two. Both sets are atomic shell
+pairs.
+
+These are usually called {e shell quartets} in the literature, but we prefer to use
+{e pair} for two functions with the same electron, and {e couple} for two functions
+acting on different electrons, since they will be coupled by a two-electron operator.
+
+
+*)
+
+type t 
+
+
+val make : ?cutoff:float -> AtomicShellPair.t -> AtomicShellPair.t -> t option
+(** Creates an atomic shell pair couple using two atomic shell pairs.
+*)
+
+val ang_mom : t -> AngularMomentum.t
+(** Total angular momentum of the shell pair couple: sum of the angular momenta of
+    all the shells.  *)
+
+
+val atomic_shell_a : t -> AtomicShell.t
+(** Returns the first atomic shell of the first shell pair. *)
+
+val atomic_shell_b : t -> AtomicShell.t
+(** Returns the second atomic shell of the first shell pair. *)
+
+val atomic_shell_c : t -> AtomicShell.t
+(** Returns the first atomic shell of the second shell pair. *)
+
+val atomic_shell_d : t -> AtomicShell.t
+(** Returns the second atomic shell of the second shell pair. *)
+
+
+val atomic_shell_pair_p : t -> AtomicShellPair.t
+(** Returns the first atomic shell pair that was used to build the shell pair. *)
+
+val atomic_shell_pair_q : t -> AtomicShellPair.t
+(** Returns the second atomic shell pair that was used to build the shell pair. *)
+
+val monocentric : t -> bool
+(** True if all four atomic shells have the same center. *)
+
+val contracted_shell_pair_couples : t -> ContractedShellPairCouple.t list
+(** Returns the list of significant contracted shell pair couples. *)
+
+val zkey_array : t -> Zkey.t array
+(** Returns the array of {!Zkey.t} relative to the four shells of the
+    shell pair couple.
+*)
+
+val norm_scales : t -> float array
+(** Scaling factors of normalization coefficients inside the shell. The
+    ordering is the same as {!zkey_array}.
+*)
+
+
diff --git a/Basis/ContractedShellPairCouple.ml b/Basis/ContractedShellPairCouple.ml
index 3c6a7ce..b3c1b5f 100644
--- a/Basis/ContractedShellPairCouple.ml
+++ b/Basis/ContractedShellPairCouple.ml
@@ -49,7 +49,7 @@ let coefs_and_shell_pair_couples t = t.coefs_and_shell_pair_couples
 
 let monocentric t =
   Csp.monocentric t.shell_pair_p && Csp.monocentric t.shell_pair_q &&
-  Csp.a_minus_b t.shell_pair_p = Csp.a_minus_b t.shell_pair_q
+  Cs.center (Csp.shell_a t.shell_pair_p) = Cs.center (Csp.shell_a t.shell_pair_q)
 
 
 let ang_mom t = t.ang_mom