Introduced shell pair couples

Basis/AtomicShell.ml

@@ -5,7 +5,7 @@ type t = {
   expo      : float array array;
   coef      : float array array;
   center    : Coordinate.t;
-  totAngMom : AngularMomentum.t;
+  ang_mom   : AngularMomentum.t;
   norm_coef : float array array;
   norm_coef_scale : float array;
   index     : int;
@@ -32,10 +32,10 @@ let make ?(index=0) contr =
   if not (unique_center (Array.length contr - 1)) then
     invalid_arg "ContractedAtomicShell.make Coordinate.t differ";
     
-  let totAngMom = Cs.totAngMom contr.(0) in
+  let ang_mom = Cs.ang_mom contr.(0) in
   let rec unique_angmom = function
   | 0 -> true
-  | i -> if Cs.totAngMom contr.(i) = totAngMom then unique_angmom (i-1) else false
+  | i -> if Cs.ang_mom contr.(i) = ang_mom then unique_angmom (i-1) else false
   in
   if not (unique_angmom (Array.length contr - 1)) then
     invalid_arg "ContractedShell.make: AngularMomentum.t differ";
@@ -45,7 +45,7 @@ let make ?(index=0) contr =
   in
   let norm_coef_scale = Cs.norm_scales contr.(0)
   in
-  { index ; expo ; coef ; center ; totAngMom ; norm_coef ;
+  { index ; expo ; coef ; center ; ang_mom ; norm_coef ;
     norm_coef_scale ; contr  }
 
 
@@ -59,7 +59,7 @@ let coefficients x = x.coef
 
 let center x = x.center
 
-let totAngMom x = x.totAngMom
+let ang_mom x = x.ang_mom
 
 let size x = Array.length x.contr
 
@@ -83,7 +83,7 @@ let pp_debug ppf x =
   fprintf ppf "@[<2>expo =@ %a ;@]@ " pp_float_2darray_size x.expo;
   fprintf ppf "@[<2>coef =@ %a ;@]@ " pp_float_2darray_size x.coef;
   fprintf ppf "@[<2>center =@ %a ;@]@ " Co.pp_angstrom x.center;
-  fprintf ppf "@[<2>totAngMom =@ %a ;@]@ " Am.pp_string x.totAngMom;
+  fprintf ppf "@[<2>ang_mom =@ %a ;@]@ " Am.pp_string x.ang_mom;
   fprintf ppf "@[<2>norm_coef =@ %a ;@]@ " pp_float_2darray_size x.norm_coef;
   fprintf ppf "@[<2>norm_coef_scale =@ %a ;@]@ " pp_float_array_size x.norm_coef_scale;
   fprintf ppf "@[<2>index =@ %d ;@]@ " x.index;
@@ -91,7 +91,7 @@ let pp_debug ppf x =
 
 let pp ppf s =
   fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+ (size_of_shell s)*(size s));
-  fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.totAngMom Co.pp s.center;
+  fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.ang_mom Co.pp s.center;
   Array.iter2 (fun e_arr c_arr -> fprintf ppf "@[<v>";
     Array.iter2 (fun e c -> fprintf ppf "@[%16.8e  %16.8e@]@;" e c)
     e_arr c_arr;

Basis/AtomicShell.mli

@@ -42,7 +42,7 @@ val index : t -> int
 val center : t -> Coordinate.t
 (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *)
 
-val totAngMom : t -> AngularMomentum.t
+val ang_mom : t -> AngularMomentum.t
 (** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *)
 
 val size : t -> int
@@ -68,7 +68,7 @@ val normalizations : t -> float array array
 
 val norm_scales : t -> float array
 (** Scaling factors {% $f(n_x,n_y,n_z)$ %}, given in the same order as
-    [AngularMomentum.zkey_array totAngMom]. *)
+    [AngularMomentum.zkey_array ang_mom]. *)
 
 val size_of_shell : t -> int
 (** Number of contracted functions in the shell: length of {!norm_coef_scale}. *)

Basis/AtomicShellPair.mli

@@ -43,7 +43,7 @@ val norm_sq : t -> float
 val norm_scales : t -> float array
   (* norm_coef.(i) / norm_coef.(0) *)
 
-val totAngMom : t -> AngularMomentum.t
+val ang_mom : t -> AngularMomentum.t
   (* Total angular Momentum *)
 
 val monocentric : t -> bool

Basis/Basis.ml

@@ -17,10 +17,10 @@ let of_nuclei_and_general_basis nucl bas =
   let contracted_shells = 
     Array.map (fun (e, center) ->
         List.assoc e bas
-        |> Array.map (fun (totAngMom, shell) ->
+        |> Array.map (fun (ang_mom, shell) ->
             let lc = 
               Array.map (fun Gb.{exponent ; coefficient} ->
-                coefficient, Ps.make totAngMom center exponent) shell
+                coefficient, Ps.make ang_mom center exponent) shell
             in
             let result = Cs.make ~index:!index_ lc in
             index_ := !index_ + Cs.size_of_shell result;
@@ -32,16 +32,16 @@ let of_nuclei_and_general_basis nucl bas =
   in
   let atomic_shells = lazy(
     let uniq_center_angmom =
-      Array.map (fun x -> Cs.center x, Cs.totAngMom x) contracted_shells
+      Array.map (fun x -> Cs.center x, Cs.ang_mom x) contracted_shells
       |> Array.to_list
       |> List.sort_uniq compare
     in
     let csl = 
       Array.to_list contracted_shells
     in
-    List.map (fun (center, totAngMom) ->
+    List.map (fun (center, ang_mom) ->
       let a = 
-        List.filter (fun x -> Cs.center x = center && Cs.totAngMom x = totAngMom) csl
+        List.filter (fun x -> Cs.center x = center && Cs.ang_mom x = ang_mom) csl
         |> Array.of_list
       in
       As.make ~index:(Cs.index a.(0)) a

Basis/ContractedShell.ml

@@ -5,7 +5,7 @@ type t = {
   expo      : float array;
   coef      : float array;
   center    : Coordinate.t;
-  totAngMom : AngularMomentum.t;
+  ang_mom   : AngularMomentum.t;
   norm_coef : float array;
   norm_coef_scale : float array;
   index     : int;
@@ -32,10 +32,10 @@ let make ?(index=0) lc =
   if not (unique_center (Array.length prim - 1)) then
     invalid_arg "ContractedShell.make Coordinate.t differ";
     
-  let totAngMom = Ps.totAngMom prim.(0) in
+  let ang_mom = Ps.ang_mom prim.(0) in
   let rec unique_angmom = function
   | 0 -> true
-  | i -> if Ps.totAngMom prim.(i) = totAngMom then unique_angmom (i-1) else false
+  | i -> if Ps.ang_mom prim.(i) = ang_mom then unique_angmom (i-1) else false
   in
   if not (unique_angmom (Array.length prim - 1)) then
     invalid_arg "ContractedShell.make: AngularMomentum.t differ";
@@ -47,7 +47,7 @@ let make ?(index=0) lc =
   in
   let norm_coef_scale = Ps.norm_scales prim.(0)
   in
-  { index ; expo ; coef ; center ; totAngMom ; norm_coef ;
+  { index ; expo ; coef ; center ; ang_mom ; norm_coef ;
     norm_coef_scale ; prim }
 
 
@@ -61,7 +61,7 @@ let coefficients x = x.coef
 
 let center x = x.center
 
-let totAngMom x = x.totAngMom
+let ang_mom x = x.ang_mom
 
 let size x = Array.length x.prim
 
@@ -85,18 +85,18 @@ let pp_debug ppf x =
   fprintf ppf "@[<2>expo =@ %a ;@]@ " pp_float_array_size x.expo;
   fprintf ppf "@[<2>coef =@ %a ;@]@ " pp_float_array_size x.coef;
   fprintf ppf "@[<2>center =@ %a ;@]@ " Co.pp_angstrom x.center;
-  fprintf ppf "@[<2>totAngMom =@ %a ;@]@ " Am.pp_string x.totAngMom;
+  fprintf ppf "@[<2>ang_mom =@ %a ;@]@ " Am.pp_string x.ang_mom;
   fprintf ppf "@[<2>norm_coef =@ %a ;@]@ " pp_float_array_size x.norm_coef;
   fprintf ppf "@[<2>norm_coef_scale =@ %a ;@]@ " pp_float_array_size x.norm_coef_scale;
   fprintf ppf "@[<2>index =@ %d ;@]@ " x.index;
   fprintf ppf "}@,@]"
 
 let pp ppf s =
-  (match s.totAngMom with
+  (match s.ang_mom with
    | Am.S -> fprintf ppf "@[%3d@]    " (s.index+1)
    | _    -> fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+(Array.length s.norm_coef_scale))
   );
-  fprintf ppf "@[%a@ %a@]@[" Am.pp_string s.totAngMom Co.pp s.center;
+  fprintf ppf "@[%a@ %a@]@[" Am.pp_string s.ang_mom Co.pp s.center;
   Array.iter2 (fun e c -> fprintf ppf "@[%16.8e  %16.8e@]@;" e c) s.expo s.coef;
   fprintf ppf "@]"
 

Basis/ContractedShell.mli

@@ -40,7 +40,7 @@ val index : t -> int
 val center : t -> Coordinate.t
 (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *)
 
-val totAngMom : t -> AngularMomentum.t
+val ang_mom : t -> AngularMomentum.t
 (** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *)
 
 val size : t -> int
@@ -60,7 +60,7 @@ val normalizations : t -> float array
 
 val norm_scales : t -> float array
 (** Scaling factors {% $f(n_x,n_y,n_z)$ %}, given in the same order as
-    [AngularMomentum.zkey_array totAngMom]. *)
+    [AngularMomentum.zkey_array ang_mom]. *)
 
 val size_of_shell : t -> int
 (** Number of contracted functions in the shell: length of {!norm_coef_scale}. *)

Basis/ContractedShellPair.ml

@@ -47,6 +47,7 @@ let make ?(cutoff=1.e-32) s_a s_b =
   | coefs_and_shell_pairs -> Some { shell_a = s_a ; shell_b = s_b ; coefs_and_shell_pairs }
 
 
+
 let shell_a x = x.shell_a
 let shell_b x = x.shell_b
 let coefs_and_shell_pairs x = x.coefs_and_shell_pairs
@@ -63,20 +64,20 @@ let exponents_inv x =
   List.map (fun (_,sp) -> Psp.exponent_inv sp) x.coefs_and_shell_pairs
   |> Array.of_list
 
-let center_ab x =
+let a_minus_b x =
   match x.coefs_and_shell_pairs with
   | [] -> assert false
   | (_,sp)::_ -> Psp.a_minus_b sp
 
-let norm_sq x = 
+let a_minus_b_sq x = 
   match x.coefs_and_shell_pairs with
   | [] -> assert false
   | (_,sp)::_ -> Psp.a_minus_b_sq sp
 
-let totAngMom x = 
+let ang_mom x = 
   match x.coefs_and_shell_pairs with
   | [] -> assert false
-  | (_,sp)::_ -> Psp.totAngMom sp
+  | (_,sp)::_ -> Psp.ang_mom sp
 
 let norm_scales x = 
   match x.coefs_and_shell_pairs with
@@ -89,6 +90,7 @@ let monocentric x =
   | (_,sp)::_ -> Psp.monocentric  sp
 
 
+
 (** Returns an integer characteristic of a contracted shell pair *)
 let hash a = 
   Array.map Hashtbl.hash a
@@ -151,22 +153,9 @@ let unique sp =
   aux [] sp
 
 
-(** A map from a shell pair hash to the list of indices in the array
+let zkey_array x =
-    of shell pairs.
+  Am.zkey_array (Am.Doublet 
-*)
+    Cs.(ang_mom x.shell_a, ang_mom x.shell_b)
-let indices sp =
+  )
-  let map =
-    Hashtbl.create 129
-  in
-  Array.iteri (fun i s ->
-    Array.iteri (fun j shell_p ->
-      let key =
-        hash shell_p
-      in
-      Hashtbl.add map key (i,j); ) s
-  ) sp;
-  map
-      
-
 
 

Basis/ContractedShellPair.mli

@@ -26,6 +26,7 @@ val make : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option
     [None].
 *)
 
+
 val of_contracted_shell_array : ContractedShell.t array -> t option array array
 (** Creates all possible contracted shell pairs from a list of contracted shells.
     If a shell pair is not significant, sets the value to [None]:
@@ -60,21 +61,25 @@ val coefficients : t -> float array
 val exponents_inv : t -> float array
 
 
-val center_ab : t -> Coordinate.t
+val a_minus_b : t -> Coordinate.t
   (* A-B *)
 
-val norm_sq : t -> float
+val a_minus_b_sq : t -> float
   (* |A-B|^2 *)
 
 val norm_scales : t -> float array
   (* normalizations.(i) / normalizations.(0) *)
 
-val totAngMom : t -> AngularMomentum.t
+val ang_mom : t -> AngularMomentum.t
   (* Total angular Momentum *)
 
 val monocentric : t -> bool
 (** If true, the two contracted shells have the same center. *)
 
+val zkey_array : t -> Zkey.t array
+(** Returns the array of Zkeys associated with the contracted shell pair. *)
+
+
 (*
 val hash : t -> int array
 val cmp : t -> t -> int

Basis/ContractedShellPairCouple.ml

@@ -0,0 +1,75 @@
+type t = 
+{
+  shell_pair_p: ContractedShellPair.t ;
+  shell_pair_q: ContractedShellPair.t ;
+  shell_a     : ContractedShell.t ;
+  shell_b     : ContractedShell.t ;
+  shell_c     : ContractedShell.t ;
+  shell_d     : ContractedShell.t ;
+  ang_mom     : AngularMomentum.t ;
+  coefs_and_shell_pair_couples : (float * PrimitiveShellPairCouple.t) list ;
+}
+
+module Am   = AngularMomentum
+module Co   = Coordinate
+module Cs   = ContractedShell
+module Csp  = ContractedShellPair
+module Pspc = PrimitiveShellPairCouple
+
+let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q = 
+  let ang_mom =
+    Am.(Csp.ang_mom shell_pair_p + Csp.ang_mom shell_pair_q)
+  in
+  let shell_a = Csp.shell_a shell_pair_p
+  and shell_b = Csp.shell_b shell_pair_p
+  and shell_c = Csp.shell_a shell_pair_q
+  and shell_d = Csp.shell_b shell_pair_q
+  in
+  let cutoff = 1.e-3 *. cutoff in
+  let coefs_and_shell_pair_couples =
+    List.map (fun (c_ab, sp_ab) -> 
+      List.map (fun (c_cd, sp_cd) -> 
+        let coef_prod = c_ab *. c_cd in
+        if abs_float coef_prod < cutoff then None
+        else Some (coef_prod, Pspc.make sp_ab sp_cd)
+      ) (Csp.coefs_and_shell_pairs shell_pair_q)
+    ) (Csp.coefs_and_shell_pairs shell_pair_p)
+    |> List.concat
+    |> List.filter (function None -> false | _ -> true)
+    |> List.map (function None -> assert false | Some x -> x)
+  in
+  match coefs_and_shell_pair_couples with
+  | [] -> None 
+  | _  -> Some {  shell_pair_p ; shell_pair_q ; ang_mom ;
+                  shell_a ; shell_b ; shell_c ; shell_d ;
+                  coefs_and_shell_pair_couples ;
+                }
+
+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
+
+
+let ang_mom t = t.ang_mom
+
+let shell_pair_p t = t.shell_pair_p
+let shell_pair_q t = t.shell_pair_q
+
+let shell_a t = t.shell_a
+let shell_b t = t.shell_b
+let shell_c t = t.shell_c
+let shell_d t = t.shell_d
+
+
+let zkey_array t =
+  match t.coefs_and_shell_pair_couples with
+  | (_,f)::_ -> Pspc.zkey_array f
+  |  _ -> invalid_arg "ContractedShellPairCouple.zkey_array"
+
+let norm_scales t = 
+  match t.coefs_and_shell_pair_couples with
+  | (_,f)::_ -> Pspc.norm_scales f
+  |  _ -> invalid_arg "ContractedShellPairCouple.norm_scales"
+

Basis/ContractedShellPairCouple.mli

@@ -0,0 +1,61 @@
+(** Data structure describing a couple of contracted shells pairs.
+
+A contracted 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 contracted 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 -> ContractedShellPair.t -> ContractedShellPair.t -> t option
+(** Creates a contracted shell pair couple using two contracted 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 shell_a : t -> ContractedShell.t
+(** Returns the first contracted shell of the first shell pair. *)
+
+val shell_b : t -> ContractedShell.t
+(** Returns the second contracted shell of the first shell pair. *)
+
+val shell_c : t -> ContractedShell.t
+(** Returns the first contracted shell of the second shell pair. *)
+
+val shell_d : t -> ContractedShell.t
+(** Returns the second contracted shell of the second shell pair. *)
+
+
+val shell_pair_p : t -> ContractedShellPair.t
+(** Returns the first contracted shell pair that was used to build the shell pair. *)
+
+val shell_pair_q : t -> ContractedShellPair.t
+(** Returns the second contracted shell pair that was used to build the shell pair. *)
+
+val monocentric : t -> bool
+(** True if all four contracted shells have the same center. *)
+
+val coefs_and_shell_pair_couples : t -> (float * PrimitiveShellPairCouple.t) list
+(** Returns the list of significant 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}.
+*)
+
+

Basis/ERI.ml

@@ -221,10 +221,10 @@ let of_basis basis =
                               )
                             else
                               out := !out + 1;
-                          ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(l))))
+                          ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(l))))
-                      ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(k))))
+                      ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(k))))
-                  ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(j))))
+                  ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(j))))
-              ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(i))))
+              ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(i))))
             with NullIntegral ->  ()
           done;
         done;
@@ -280,7 +280,7 @@ let to_file ~filename basis =
   let oc = open_out filename in
   let zkey = Array.map (fun b -> 
       let result =
-        Angular_momentum.(zkey_array (Kind_1 b.totAngMom))
+        Angular_momentum.(zkey_array (Kind_1 b.ang_mom))
       in
       { n=Array.length result ;  cls=result }
     ) basis
@@ -378,7 +378,7 @@ let to_file ~filename basis =
 let xto_file ~filename basis =
   let zkey = Array.map (fun b -> 
       let result =
-        Angular_momentum.(zkey_array (Kind_1 b.totAngMom))
+        Angular_momentum.(zkey_array (Kind_1 b.ang_mom))
       in
       { n=Array.length result ;  cls=result }
     ) basis

Basis/KinInt.ml

@@ -28,9 +28,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
   | Some shell_p ->
     begin
       (* Pre-computation of integral class indices *)
-      let class_indices = 
+      let class_indices =  Csp.zkey_array shell_p in
-        Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b))
-      in
 
       let contracted_class = 
         Array.make (Array.length class_indices) 0.
@@ -39,8 +37,8 @@ let contracted_class shell_a shell_b : float Zmap.t =
       (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
       let sp = Csp.shell_pairs shell_p in
-      let center_ab = 
+      let a_minus_b = 
-        Csp.center_ab shell_p
+        Csp.a_minus_b shell_p
       in
       let norm_coef_scales =
         Csp.norm_scales shell_p
@@ -78,7 +76,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
                   let xyz = xyz_of_int k in
                   Overlap_primitives.hvrr (a, b)
                     expo_inv
-                    (Co.get xyz center_ab,
+                    (Co.get xyz a_minus_b,
                     Co.get xyz center_pa) 
                 in
                 let f k = 
@@ -154,8 +152,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
+            ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i))))
-        ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) 
+        ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) 
     done;
   done;
   Mat.detri result;

Basis/NucInt.ml

@@ -87,8 +87,8 @@ let of_basis_nuclei basis nuclei =
                 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.ang_mom shell.(j))))
-          ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
+          ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i))))
     done;
   done;
   Mat.detri eni_array;

Basis/OneElectronRR.ml

@@ -101,15 +101,10 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
   let shell_a = Csp.shell_a shell_p
   and shell_b = Csp.shell_b shell_p
   in
-  let maxm = 
+  let maxm = Am.to_int (Csp.ang_mom shell_p) in
-    Am.(Cs.totAngMom shell_a + Cs.totAngMom shell_b) 
-    |> Am.to_int
-  in
 
   (* Pre-computation of integral class indices *)
-  let class_indices = 
+  let class_indices = Csp.zkey_array shell_p in
-    Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b))
-  in
 
   let contracted_class = 
     Array.make (Array.length class_indices) 0.;
@@ -119,7 +114,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
 
   let norm_scales_p = Csp.norm_scales shell_p in
 
-  let center_ab = Csp.center_ab shell_p  in
+  let center_ab = Csp.a_minus_b shell_p  in
 
   List.iter (fun (coef_prod, psp) ->
     try
@@ -146,7 +141,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
           let zero_m_array = 
             zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
           in
-          match Cs.(totAngMom shell_a, totAngMom shell_b) with
+          match Cs.(ang_mom shell_a, ang_mom shell_b) with
           | Am.(S,S) -> 
             let integral = zero_m_array.(0) in
             contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge

Basis/Overlap.ml

@@ -28,16 +28,14 @@ let contracted_class shell_a shell_b : float Zmap.t =
     begin
 
   (* Pre-computation of integral class indices *)
-  let class_indices = 
+  let class_indices = Csp.zkey_array shell_p in
-    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 center_ab = 
+  let a_minus_b = 
-    Csp.center_ab shell_p
+    Csp.a_minus_b shell_p
   in
   let norm_coef_scales =
     Csp.norm_scales shell_p
@@ -66,7 +64,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
                 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 a_minus_b,
                   Co.get xyz center_pa) 
               in
               let norm =  norm_coef_scales.(i) in
@@ -120,8 +118,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
+            ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i))))
-        ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j))))
+        ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j))))
     done;
   done;
   Mat.detri result;

Basis/PrimitiveShell.ml

@@ -7,15 +7,15 @@ type t = {
   normalization : float;
   norm_scales   : float array lazy_t;
   center    : Coordinate.t;
-  totAngMom : AngularMomentum.t;
+  ang_mom   : AngularMomentum.t;
 }
 
 module Am = AngularMomentum
 
 
-let compute_norm_coef alpha totAngMom =
+let compute_norm_coef alpha ang_mom =
   let atot =
-    Am.to_int totAngMom
+    Am.to_int ang_mom
   in
   let factor int_array = 
        let dfa = Array.map (fun j ->
@@ -37,15 +37,15 @@ let compute_norm_coef alpha totAngMom =
   in f
 
 
-let make totAngMom center exponent =
+let make ang_mom center exponent =
   let norm_coef_func =
-    compute_norm_coef exponent totAngMom
+    compute_norm_coef exponent ang_mom
   in
   let norm =
-    1. /. norm_coef_func [| Am.to_int totAngMom ; 0 ; 0 |] 
+    1. /. norm_coef_func [| Am.to_int ang_mom ; 0 ; 0 |] 
   in
   let powers =
-    Am.zkey_array (Am.Singlet totAngMom) 
+    Am.zkey_array (Am.Singlet ang_mom) 
   in
   let norm_scales = lazy (
     Array.map (fun a ->
@@ -53,12 +53,12 @@ let make totAngMom center exponent =
     ) powers )
   in
   let normalization = 1. /. norm in
-  { exponent ; normalization ; norm_scales ; center ; totAngMom }
+  { exponent ; normalization ; norm_scales ; center ; ang_mom }
 
 
 let to_string s =
   let coord = s.center in
-  Printf.sprintf "%1s %8.3f %8.3f %8.3f  %16.8e" (Am.to_string s.totAngMom)
+  Printf.sprintf "%1s %8.3f %8.3f %8.3f  %16.8e" (Am.to_string s.ang_mom)
       (get X coord) (get Y coord) (get Z coord) s.exponent
 
 
@@ -74,7 +74,7 @@ let exponent x = x.exponent
 
 let center x = x.center
 
-let totAngMom x = x.totAngMom
+let ang_mom x = x.ang_mom
 
 let norm x = 1. /. x.normalization
 

Basis/PrimitiveShell.mli

@@ -33,7 +33,7 @@ val exponent : t -> float
 val center : t -> Coordinate.t
 (** Coordinate {% $\mathbf{A}$ %}.of the center. *)
 
-val totAngMom : t -> AngularMomentum.t
+val ang_mom : t -> AngularMomentum.t
 (** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *)
 
 val norm : t -> float 
@@ -54,7 +54,7 @@ val norm_scales : t -> float array
 (** Scaling factors {% $f(n_x,n_y,n_z)$ %} adjusting the normalization coefficient
     for the powers of {% $x,y,z$ %}. The normalization coefficients of the
     functions of the shell are given by {% $\mathcal{N}\times f$ %}.  They are
-    given in the same order as [AngularMomentum.zkey_array totAngMom]: 
+    given in the same order as [AngularMomentum.zkey_array ang_mom]: 
     {% \\[ f(n_x,n_y,n_z) = \frac{|| g_{l,0,0}(\mathbf{r}) ||}{|| g_{n_x,n_y,n_z}(\mathbf{r}) ||} \\] %}
 *)
 

Basis/PrimitiveShellPair.ml

@@ -12,7 +12,7 @@ type t = {
   norm_scales     : float array lazy_t;
   normalization   : float;              (* norm_coef_a * norm_coef_b * g, with
                      g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *)
-  totAngMom       : AngularMomentum.t; 
+  ang_mom         : AngularMomentum.t; 
   shell_a         : PrimitiveShell.t;
   shell_b         : PrimitiveShell.t;
 }
@@ -30,7 +30,7 @@ let hash a =
 
 let equivalent a b =
   a.exponent = b.exponent &&
-  a.totAngMom = b.totAngMom &&
+  a.ang_mom = b.ang_mom &&
   a.normalization = b.normalization &&
   a.center = b.center &&
   a.center_minus_a = b.center_minus_a &&
@@ -59,8 +59,8 @@ let create_make_of p_a p_b =
     |> Array.concat
   ) in
 
-  let totAngMom =
+  let ang_mom =
-    Am.( Ps.totAngMom p_a + Ps.totAngMom p_b )
+    Am.( Ps.ang_mom p_a + Ps.ang_mom p_b )
   in
 
   function p_a -> 
@@ -109,7 +109,7 @@ let create_make_of p_a p_b =
                   in
 
                   Some {
-                    totAngMom ; 
+                    ang_mom ; 
                     exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ; 
                     a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a;
                     shell_b = p_b }
@@ -136,7 +136,7 @@ let monocentric x =
   Ps.center x.shell_a = Ps.center x.shell_b
 
 
-let totAngMom x = x.totAngMom
+let ang_mom x = x.ang_mom
 
 let a_minus_b x = x.a_minus_b
 
@@ -154,3 +154,11 @@ let shell_a x = x.shell_a
 
 let shell_b x = x.shell_b
 
+
+let zkey_array x =
+  Am.zkey_array (Am.Doublet 
+    Ps.(ang_mom x.shell_a, ang_mom x.shell_b)
+  )
+
+
+

Basis/PrimitiveShellPair.mli

@@ -1,7 +1,7 @@
 (** Data structure describing a pair of primitive shells.
 
 A primitive shell pair is the cartesian product between two sets of functions, each
-set containing all the functions of a primitive shell.
+set containing all the functions of a primitive shell. These are one-electron functions.
 
 {% \\[
   \left\\{ p_{k_x,k_y,k_z}(\mathbf{r}) \right\\} =
@@ -24,7 +24,7 @@ where
 
 
 {!a_minus_b}, {!a_minus_b_sq} and {!norm_coef_scale} depend only on the
-centering of the two shells, and {!totAngMom} only depends on the angular
+centering of the two shells, and {!ang_mom} only depends on the angular
 momenta of the two shells. Hence, these quantities need to be computed only
 once when a {!ContractedShellPair.t} is built. Hence, there is the
 {!create_make_of} function which creates a [make] function which is suitable
@@ -52,7 +52,7 @@ val create_make_of : PrimitiveShell.t -> PrimitiveShell.t ->
 
 *)
 
-val totAngMom : t -> AngularMomentum.t
+val ang_mom : t -> AngularMomentum.t
 (** Total angular momentum of the shell pair: sum of the angular momenta of
     the shells.  *)
 
@@ -104,3 +104,6 @@ val hash : t -> int
 val cmp : t -> t -> int
 (** Arbitray comparison function for sorting. *)
 
+val zkey_array : t -> Zkey.t array
+(** Returns the array of Zkeys associated with the shell pair. *)
+

Basis/PrimitiveShellPairCouple.ml

@@ -0,0 +1,71 @@
+exception NullContribution
+
+type t = 
+{
+  shell_pair_p: PrimitiveShellPair.t ;
+  shell_pair_q: PrimitiveShellPair.t ;
+  shell_a     : PrimitiveShell.t ;
+  shell_b     : PrimitiveShell.t ;
+  shell_c     : PrimitiveShell.t ;
+  shell_d     : PrimitiveShell.t ;
+  ang_mom     : AngularMomentum.t ;
+  zkey_array  : Zkey.t array lazy_t;
+}
+
+module Am  = AngularMomentum
+module Co  = Coordinate
+module Ps  = PrimitiveShell
+module Psp = PrimitiveShellPair
+
+let make shell_pair_p shell_pair_q = 
+  let ang_mom =
+    Am.(Psp.ang_mom shell_pair_p + Psp.ang_mom shell_pair_q)
+  in
+  let shell_a = Psp.shell_a shell_pair_p
+  and shell_b = Psp.shell_b shell_pair_p
+  and shell_c = Psp.shell_a shell_pair_q
+  and shell_d = Psp.shell_b shell_pair_q
+  in
+  let zkey_array = lazy (
+    Am.zkey_array (Am.Quartet
+      Ps.(ang_mom shell_a, ang_mom shell_b,
+          ang_mom shell_c, ang_mom shell_d)
+    )
+  )
+  in
+  { shell_pair_p ; shell_pair_q ; ang_mom ;
+    shell_a ; shell_b ; shell_c ; shell_d ;
+    zkey_array ;
+  }
+
+
+
+let monocentric t =
+  Psp.monocentric t.shell_pair_p && Psp.monocentric t.shell_pair_q &&
+  Psp.center t.shell_pair_p = Psp.center t.shell_pair_q
+
+
+let ang_mom t = t.ang_mom
+
+let shell_pair_p t = t.shell_pair_p
+let shell_pair_q t = t.shell_pair_q
+
+let shell_a t = t.shell_a
+let shell_b t = t.shell_b
+let shell_c t = t.shell_c
+let shell_d t = t.shell_d
+
+let p_minus_q t = 
+  let p = Psp.center t.shell_pair_p
+  and q = Psp.center t.shell_pair_q
+  in Co.(p |- q)
+
+let zkey_array t = Lazy.force t.zkey_array
+
+let norm_scales t = 
+  let norm_coef_scale_p_list = Array.to_list (Psp.norm_scales t.shell_pair_p) in
+  let norm_coef_scale_q = Psp.norm_scales t.shell_pair_q in
+  List.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q)
+    norm_coef_scale_p_list
+  |> Array.concat
+

Basis/PrimitiveShellPairCouple.mli

@@ -0,0 +1,62 @@
+(** Data structure describing a couple of primitive shells pairs.
+
+A primitive 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 primitive 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 : PrimitiveShellPair.t -> PrimitiveShellPair.t -> t 
+(** Creates a primitive shell pair couple using two primitive 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 shell_a : t -> PrimitiveShell.t
+(** Returns the first primitive shell of the first shell pair. *)
+
+val shell_b : t -> PrimitiveShell.t
+(** Returns the second primitive shell of the first shell pair. *)
+
+val shell_c : t -> PrimitiveShell.t
+(** Returns the first primitive shell of the second shell pair. *)
+
+val shell_d : t -> PrimitiveShell.t
+(** Returns the second primitive shell of the second shell pair. *)
+
+
+val shell_pair_p : t -> PrimitiveShellPair.t
+(** Returns the first primitive shell pair that was used to build the shell pair. *)
+
+val shell_pair_q : t -> PrimitiveShellPair.t
+(** Returns the second primitive shell pair that was used to build the shell pair. *)
+
+val p_minus_q : t -> Coordinate.t
+(** {% $\mathbf{P-Q}$ %}. *)
+
+
+val monocentric : t -> bool
+(** True if all four primitive shells have the same center. *)
+
+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}.
+*)
+
+

Basis/TwoElectronRR.ml

@@ -4,12 +4,15 @@ module Am  = AngularMomentum
 module Co   = Coordinate
 module Cs   = ContractedShell
 module Csp  = ContractedShellPair
+module Cspc = ContractedShellPairCouple
 module Po   = Powers
 module Psp  = PrimitiveShellPair
+module Pspc = PrimitiveShellPairCouple
 module Ps   = PrimitiveShell
 
 let cutoff = Constants.integrals_cutoff 
 let cutoff2 = cutoff *. cutoff
+let empty = Zmap.create 0
 
 exception NullQuartet
 
@@ -275,50 +278,39 @@ let rec hvrr_two_e
 
 let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
 
-  let shell_a = Csp.shell_a     shell_p
+  match Cspc.make ~cutoff shell_p shell_q with
-  and shell_b = Csp.shell_b     shell_p
+  | None -> empty
-  and shell_c = Csp.shell_a     shell_q
+  | Some shell_pair_couple ->
-  and shell_d = Csp.shell_b     shell_q
+
-  in
+    let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in
-  let maxm = Am.(Csp.totAngMom shell_p + Csp.totAngMom shell_q |> to_int) in
 
     (* Pre-computation of integral class indices *)
-  let class_indices = 
+    let class_indices = Cspc.zkey_array shell_pair_couple in
-    Am.zkey_array (Am.Quartet
-         Cs.(totAngMom shell_a, totAngMom shell_b,
-             totAngMom shell_c, totAngMom shell_d ))
-  in
 
     let contracted_class = 
       Array.make (Array.length class_indices) 0.;
     in
 
     let monocentric =
-      Csp.monocentric shell_p && Csp.monocentric shell_q &&
+      Cspc.monocentric shell_pair_couple
-      Cs.center (Csp.shell_a shell_p) = Cs.center (Csp.shell_a shell_q) 
     in
 
     (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
-  let norm_coef_scale_p_list = Array.to_list (Csp.norm_scales shell_p) in
+    let center_ab = Csp.a_minus_b shell_p
-  let norm_coef_scale_q = Csp.norm_scales shell_q in
+    and center_cd = Csp.a_minus_b shell_q
-
-  let center_ab = Csp.center_ab shell_p in
-  List.iter (fun (c_ab, sp_ab) ->
-
-    let expo_b = Ps.exponent (Psp.shell_b sp_ab) 
-    and expo_inv_p = Psp.exponent_inv sp_ab
-    and center_pa = Psp.center_minus_a sp_ab
     in
 
-    List.iter (fun (c_cd, sp_cd) ->
+    let norm_scales = Cspc.norm_scales shell_pair_couple in
 
-      let coef_prod = c_ab *. c_cd in
+    List.iter (fun (coef_prod, spc) ->
 
-      (** Screening on the product of coefficients *)
+        let sp_ab = Pspc.shell_pair_p spc
-      try
+        and sp_cd = Pspc.shell_pair_q spc
-        if (abs_float coef_prod) < 1.e-3 *. cutoff then
+        in
-          raise NullQuartet;
+
+        let expo_inv_p = Psp.exponent_inv sp_ab
+        in
 
         let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in
         let norm_pq_sq = Co.dot center_pq center_pq in
@@ -330,23 +322,19 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
         in
 
         begin
-          match Cs.(totAngMom shell_a, totAngMom shell_b,
+          match Cspc.ang_mom shell_pair_couple with
-                    totAngMom shell_c, totAngMom shell_d) with
+          | Am.S ->
-          | Am.(S,S,S,S) ->
+            let integral = zero_m_array.(0) in
-            let integral =
-              zero_m_array.(0)
-            in
             contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
           | _ -> 
-            let expo_d = Ps.exponent (Psp.shell_b sp_cd) in
+            let expo_b = Ps.exponent (Psp.shell_b sp_ab) 
-            let map_1d = Zmap.create (4*maxm) in
+            and expo_d = Ps.exponent (Psp.shell_b sp_cd)
-            let map_2d = Zmap.create (Array.length class_indices) in
+            and center_pa = Psp.center_minus_a sp_ab
-            let center_cd = Psp.a_minus_b sp_cd in
+            in
-            let center_qc = Psp.center_minus_a sp_cd in
+            let map_1d = Zmap.create (4*maxm) 
-            let norm_coef_scale =
+            and map_2d = Zmap.create (Array.length class_indices)
-              List.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) 
+            in
-                norm_coef_scale_p_list
+            let center_qc = Psp.center_minus_a sp_cd
-              |> Array.concat
             in
 
             (* Compute the integral class from the primitive shell quartet *)
@@ -393,7 +381,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                   *)
 
 
-                  let norm =  norm_coef_scale.(i) in
+                  let norm =  norm_scales.(i) in
                   let coef_prod = coef_prod *. norm in
 
                   let integral = 
@@ -410,9 +398,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                 with NullQuartet -> ()
               )
         end
-      with NullQuartet -> ()
+      ) (Cspc.coefs_and_shell_pair_couples shell_pair_couple);
-    )  (Csp.coefs_and_shell_pairs shell_q)
-  ) (Csp.coefs_and_shell_pairs shell_p);
 
     let result = 
       Zmap.create (Array.length contracted_class)
@@ -431,6 +417,6 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
   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
+  | _ -> empty
 
 

Basis/TwoElectronRRVectorized.ml

@@ -6,6 +6,7 @@ module Am  = AngularMomentum
 module Co   = Coordinate
 module Cs   = ContractedShell
 module Csp  = ContractedShellPair
+module Cspc = ContractedShellPairCouple
 module Po   = Powers
 module Psp  = PrimitiveShellPair
 module Ps   = PrimitiveShell
@@ -15,6 +16,7 @@ exception Found
 
 let cutoff = Constants.integrals_cutoff 
 let cutoff2 = cutoff *. cutoff
+let empty = Zmap.create 0
 
 let at_least_one_valid arr =
   try
@@ -555,57 +557,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 
 let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
 
-  let shell_a = Csp.shell_a     shell_p
+  let sp = Csp.shell_pairs shell_p
-  and shell_b = Csp.shell_b     shell_p
+  and sq = Csp.shell_pairs shell_q
-  and shell_c = Csp.shell_a     shell_q
+  and cp = Csp.coefficients shell_p
-  and shell_d = Csp.shell_b     shell_q
+  and cq = Csp.coefficients shell_q
-  in
-  let maxm = 
-    Am.(Csp.totAngMom shell_p + Csp.totAngMom shell_q |> to_int)
-  in
-
-  (* Pre-computation of integral class indices *)
-  let class_indices = 
-    Am.zkey_array (Am.Quartet
-      Cs.(totAngMom shell_a, totAngMom shell_b,
-          totAngMom shell_c, totAngMom shell_d))
-  in
-
-  let contracted_class = 
-    Array.make (Array.length class_indices) 0.;
-  in
-
-  let monocentric =
-      Csp.monocentric shell_p && Csp.monocentric shell_q &&
-      Cs.center (Csp.shell_a shell_p) =
-      Cs.center (Csp.shell_a shell_q)
-  in
-
-  (** Screening on the product of coefficients *)
-  let coef_max_p =
-    List.fold_left (fun accu (x,_) ->
-      if (abs_float x) > accu then (abs_float x) else accu)
-      0. (Csp.coefs_and_shell_pairs shell_p)
-  and coef_max_q =
-    List.fold_left (fun accu (x,_) ->
-      if (abs_float x) > accu then (abs_float x) else accu)
-      0. (Csp.coefs_and_shell_pairs shell_q)
-  in
-
-  let filter_p =
-    let cutoff_p = cutoff /. coef_max_p in
-    Csp.coefs_and_shell_pairs shell_p
-    |> List.filter (fun (ck,f) -> abs_float ck > cutoff_p) 
-  and filter_q =
-    let cutoff_q = cutoff /. coef_max_q in
-    Csp.coefs_and_shell_pairs shell_q
-    |> List.filter (fun (ck,f) -> abs_float ck > cutoff_q)
-  in
-
-  let sp = List.map snd filter_p |> Array.of_list
-  and sq = List.map snd filter_q |> Array.of_list
-  and cp = List.map fst filter_p |> Array.of_list
-  and cq = List.map fst filter_q |> Array.of_list
   in
 
   let np, nq = 
@@ -613,13 +568,34 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
     Array.length sq
   in
 
+  try
+    match Cspc.make ~cutoff shell_p shell_q with
+    | None -> raise NullQuartet
+    | Some shell_pair_couple -> 
+
+      let shell_a = Cspc.shell_a shell_pair_couple
+      and shell_c = Cspc.shell_c shell_pair_couple
+      in
+
+      let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in
+
+
+      (* Pre-computation of integral class indices *)
+      let class_indices = Cspc.zkey_array shell_pair_couple in
+
+      let contracted_class = 
+        Array.make (Array.length class_indices) 0.;
+      in
+
+      let monocentric =
+        Cspc.monocentric shell_pair_couple
+      in
 
       (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
       begin
-    match Cs.(totAngMom shell_a, totAngMom shell_b,
+        match Cspc.ang_mom shell_pair_couple with
-              totAngMom shell_c, totAngMom shell_d) with
+        | Am.S ->
-    | Am.(S,S,S,S) ->
           contracted_class.(0) <-
             begin
               try
@@ -665,6 +641,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
             Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) )
           in
 
+          let norm = Cspc.norm_scales shell_pair_couple in
+
           let expo_inv_p = 
             Array.map (fun shell_ab -> Psp.exponent_inv shell_ab) sp
           and expo_inv_q = 
@@ -676,7 +654,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           and expo_d =
             Array.map (fun shell_cd -> Ps.exponent (Psp.shell_b shell_cd) ) sq
           in
-      let norm_coef_scale_p = Csp.norm_scales shell_p in
 
           let center_pq =
             let result = 
@@ -778,16 +755,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
             result
           in
 
-        let norm = 
-          let norm_coef_scale_q =
-            Csp.norm_scales shell_q
-          in
-          Array.to_list norm_coef_scale_p
-          |> List.map (fun v1 ->
-              Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q)
-          |> Array.concat
-        in
-
           let map_1d = Zmap.create (4*maxm)
           and map_2d = Array.init (maxm+1) (fun _ -> Zmap.create (Array.length class_indices))
           in
@@ -836,8 +803,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                     zero_m_array
                     (expo_b, expo_d)
                     (expo_inv_p, expo_inv_q)
-                      (Csp.center_ab shell_p,
+                    (Csp.a_minus_b shell_p,
-                       Csp.center_ab shell_q, center_pq)
+                     Csp.a_minus_b shell_q, center_pq)
                     (center_pa, center_qc)
                     map_1d map_2d np nq
                 in
@@ -852,6 +819,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
       in
       Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
       result
+  with NullQuartet -> empty
 
 
 
@@ -864,5 +832,5 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
   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
+  | _ -> empty