Introduced atomic_shell_pair_couples

Basis/AtomicShellPair.ml

@@ -15,12 +15,11 @@ module Am  = AngularMomentum
 module As  = AtomicShell
 module Co  = Coordinate
 module Cs  = ContractedShell
-module Ps  = PrimitiveShell
+module Csp = ContractedShellPair
-module Psp = PrimitiveShellPair
 
 (** Creates an atomic shell pair : an array of pairs of contracted shells.
 *)
-let make ?(cutoff=1.e-32) atomic_shell_a atomic_shell_b =
+let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b =
 
   let l_a = Array.to_list (As.contracted_shells atomic_shell_a)
   and l_b = Array.to_list (As.contracted_shells atomic_shell_b)
@@ -32,27 +31,33 @@ let make ?(cutoff=1.e-32) atomic_shell_a atomic_shell_b =
         Csp.make ~cutoff s_a s_b
       ) l_b
     ) l_a
+    |> List.concat
+    |> List.filter (function None -> false | _ -> true)
+    |> List.map (function None -> assert false | Some x -> x)
   in
+  match contracted_shell_pairs with
+  | [] -> None
+  | _  -> Some { atomic_shell_a ; atomic_shell_b ; contracted_shell_pairs }
 
 
 let atomic_shell_a  x = x.atomic_shell_a
 let atomic_shell_b  x = x.atomic_shell_b
-let contracted_shell_pairs = x.contracted_shell_pairs
+let contracted_shell_pairs x = x.contracted_shell_pairs
 
 let monocentric x = Csp.monocentric @@ List.hd x.contracted_shell_pairs
 
-let center_ab x = Csp.center_ab @@ List.hd x.contracted_shell_pairs
+let a_minus_b x = Csp.a_minus_b @@ List.hd x.contracted_shell_pairs
 
-let totAngMon x = Csp.totAngMon @@ List.hd x.contracted_shell_pairs
+let a_minus_b_sq x = Csp.a_minus_b_sq @@ List.hd x.contracted_shell_pairs
+
+let ang_mom x = Csp.ang_mom @@ List.hd x.contracted_shell_pairs
 
 let norm_scales x = Csp.norm_scales @@ List.hd x.contracted_shell_pairs
 
-let norm_sq x = Csp.norm_sq @@ List.hd x.contracted_shell_pairs
-
 (** The array of all shell pairs with their correspondance in the list
     of contracted shells.
 *)
-let of_atomic_shell_array basis =
+let of_atomic_shell_array ?(cutoff=Constants.epsilon) basis =
   Array.mapi (fun i shell_a ->
     Array.mapi (fun j shell_b -> 
       make ~cutoff shell_a shell_b)

Basis/AtomicShellPair.mli

@@ -14,7 +14,7 @@ val make : ?cutoff:float -> AtomicShell.t -> AtomicShell.t -> t option
     [None].
 *)
 
-val of_atomic_shell_array : AtomicShell.t array -> t option list
+val of_atomic_shell_array : ?cutoff:float -> AtomicShell.t array -> t option array array
 (** Creates all possible atomic shell pairs from an array of atomic shells.
     If an atomic shell pair is not significant, sets the value to [None].
 *)
@@ -34,10 +34,10 @@ val contracted_shell_pairs : t -> ContractedShellPair.t list
     contracted functions used to build the atomic shell pair.
 *)
 
-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

Basis/ContractedShellPair.ml

@@ -21,7 +21,7 @@ module Psp = PrimitiveShellPair
     A contracted shell with N functions combined with a contracted
     shell with M functions generates a NxM array of shell pairs.
 *)
-let make ?(cutoff=1.e-32) s_a s_b =
+let make ?(cutoff=Constants.epsilon) s_a s_b =
 
   let make = Psp.create_make_of (Cs.primitives s_a).(0) (Cs.primitives s_b).(0) in
 
@@ -119,10 +119,10 @@ let cmp a b =
 (** The array of all shell pairs with their correspondance in the list
     of contracted shells.
 *)
-let of_contracted_shell_array basis =
+let of_contracted_shell_array ?(cutoff=Constants.epsilon) basis =
   Array.mapi (fun i shell_a ->
     Array.mapi (fun j shell_b -> 
-      make shell_a shell_b)
+      make ~cutoff shell_a shell_b)
     (Array.sub basis 0 (i+1)) 
   ) basis
 

Basis/ContractedShellPair.mli

@@ -27,7 +27,7 @@ val make : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option
 *)
 
 
-val of_contracted_shell_array : ContractedShell.t array -> t option array array
+val of_contracted_shell_array : ?cutoff:float -> 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]:
 

Basis/ERI.ml

@@ -42,20 +42,12 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq =
 
 
 (** Compute all the integrals of a contracted class *)
-(*
-let contracted_class  shell_a shell_b shell_c shell_d : float Zmap.t =
-  TwoElectronRRVectorized.contracted_class ~zero_m shell_a shell_b shell_c shell_d 
-  *)
-let contracted_class  shell_a shell_b shell_c shell_d : float Zmap.t =
-  TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d 
-
-(** Compute all the integrals of a contracted class *)
-let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
-  TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
-
 let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
   TwoElectronRR.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
+let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
+  TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
+
 
 let cutoff2 = cutoff *. cutoff
 (*

Basis/TwoElectronRR.ml

@@ -1,6 +1,8 @@
 open Util
 
 module Am   = AngularMomentum
+module Asp  = AtomicShellPair
+module Aspc = AtomicShellPairCouple
 module Co   = Coordinate
 module Cs   = ContractedShell
 module Csp  = ContractedShellPair
@@ -408,15 +410,140 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
 
 
-(** Computes all the two-electron integrals of the contracted shell quartet *)
-let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
-
-  let shell_p = Csp.make ~cutoff shell_a shell_b
-  and shell_q = Csp.make ~cutoff shell_c shell_d
-  in
-  match shell_p, shell_q with
-  | Some shell_p, Some shell_q -> 
-    contracted_class_shell_pairs ~zero_m shell_p shell_q 
-  | _ -> empty
 
 
+
+
+
+let contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
+
+  match Aspc.make shell_p shell_q with
+  | None -> empty
+  | Some atomic_shell_pair_couple ->
+
+    let maxm = Am.to_int (Aspc.ang_mom atomic_shell_pair_couple) in
+
+    (* Pre-computation of integral class indices *)
+    let class_indices = Aspc.zkey_array atomic_shell_pair_couple in
+
+    let contracted_class = 
+      Array.make (Array.length class_indices) 0.;
+    in
+
+    let monocentric =
+      Aspc.monocentric atomic_shell_pair_couple
+    in
+
+    (* Compute all integrals in the shell for each pair of significant shell pairs *)
+
+    let center_ab = Asp.a_minus_b shell_p
+    and center_cd = Asp.a_minus_b shell_q
+    in
+
+    let norm_scales = Aspc.norm_scales atomic_shell_pair_couple in
+
+
+    List.iter (fun cspc ->
+      List.iter (fun (coef_prod, spc) ->
+        let sp_ab = Pspc.shell_pair_p spc
+        and sp_cd = Pspc.shell_pair_q spc
+        in
+
+        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
+        let expo_inv_q = Psp.exponent_inv sp_cd in
+        let expo_pq_inv = expo_inv_p +. expo_inv_q in
+
+        let zero_m_array = 
+          zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
+        in
+
+        begin
+          match Aspc.ang_mom atomic_shell_pair_couple with
+          | Am.S ->
+            let integral = zero_m_array.(0) in
+            contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
+          | _ -> 
+            let expo_b = Ps.exponent (Psp.shell_b sp_ab) 
+            and expo_d = Ps.exponent (Psp.shell_b sp_cd)
+            and center_pa = Psp.center_minus_a sp_ab
+            in
+            let map_1d = Zmap.create (4*maxm) 
+            and map_2d = Zmap.create (Array.length class_indices)
+            in
+            let center_qc = Psp.center_minus_a sp_cd
+            in
+
+            (* Compute the integral class from the primitive shell quartet *)
+            class_indices
+            |> Array.iteri (fun i key ->
+                let (angMom_a,angMom_b,angMom_c,angMom_d) =
+                  match Zkey.to_powers key with
+                  | Zkey.Twelve x -> x
+                  | _             -> assert false
+                in
+                try
+                  if monocentric then
+                    begin
+                      if (  ((1 land angMom_a.Po.x + angMom_b.Po.x + angMom_c.Po.x + angMom_d.Po.x)=1) ||
+                            ((1 land angMom_a.Po.y + angMom_b.Po.y + angMom_c.Po.y + angMom_d.Po.y)=1) ||
+                            ((1 land angMom_a.Po.z + angMom_b.Po.z + angMom_c.Po.z + angMom_d.Po.z)=1) 
+                         ) then
+                        raise NullQuartet
+                    end;
+
+                  (* Schwartz screening *)
+                  (*
+                  if (maxm > 8) then
+                  (
+                    let schwartz_p = 
+                      let key =
+                        Zkey.of_powers_twelve angMom_a angMom_b angMom_a angMom_b
+                      in
+                      match schwartz_p with
+                      | None -> 1.
+                      | Some schwartz_p -> Zmap.find schwartz_p key
+                    in
+                    if schwartz_p < cutoff then raise NullQuartet;
+                    let schwartz_q = 
+                      let key =
+                        Zkey.of_powers_twelve angMom_c angMom_d angMom_c angMom_d
+                      in
+                      match schwartz_q with
+                      | None -> 1.
+                      | Some schwartz_q -> Zmap.find schwartz_q key
+                    in
+                    if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet;
+                  );
+                  *)
+
+
+                  let norm =  norm_scales.(i) in
+                  let coef_prod = coef_prod *. norm in
+
+                  let integral = 
+                    hvrr_two_e
+                      angMom_a angMom_b angMom_c angMom_d
+                      zero_m_array
+                      expo_b  expo_d
+                      expo_inv_p expo_inv_q
+                      center_ab center_cd center_pq
+                      center_pa center_qc 
+                      map_1d map_2d 
+                  in
+                  contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+                with NullQuartet -> ()
+              )
+        end
+        ) (Cspc.coefs_and_shell_pair_couples cspc)
+      ) (Aspc.contracted_shell_pair_couples atomic_shell_pair_couple);
+
+    let result = 
+      Zmap.create (Array.length contracted_class)
+    in
+    Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
+    result
+

Basis/TwoElectronRRVectorized.ml

@@ -822,15 +822,3 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   with NullQuartet -> empty
 
 
-
-(** Computes all the two-electron integrals of the contracted shell quartet *)
-let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
-
-  let shell_p = Csp.make ~cutoff shell_a shell_b
-  and shell_q = Csp.make ~cutoff shell_c shell_d
-  in
-  match shell_p, shell_q with
-  | Some shell_p, Some shell_q ->
-    contracted_class_shell_pairs ~zero_m shell_p shell_q
-  | _ -> empty
-