Introducing AtomicShellPair.ml

Basis/AtomicShellPair.ml

@@ -0,0 +1,63 @@
+open Util
+open Constants
+
+exception Null_contribution
+
+type t =
+{
+  atomic_shell_a         : AtomicShell.t;
+  atomic_shell_b         : AtomicShell.t;
+  contracted_shell_pairs : ContractedShellPair.t list;
+}
+
+
+module Am  = AngularMomentum
+module As  = AtomicShell
+module Co  = Coordinate
+module Cs  = ContractedShell
+module Ps  = PrimitiveShell
+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 l_a = Array.to_list (As.contracted_shells atomic_shell_a)
+  and l_b = Array.to_list (As.contracted_shells atomic_shell_b)
+  in
+
+  let contracted_shell_pairs =
+    List.map (fun s_a -> 
+      List.map (fun s_b ->
+        Csp.make ~cutoff s_a s_b
+      ) l_b
+    ) l_a
+  in
+
+
+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 monocentric x = Csp.monocentric @@ List.hd x.contracted_shell_pairs
+
+let center_ab x = Csp.center_ab @@ List.hd x.contracted_shell_pairs
+
+let totAngMon x = Csp.totAngMon @@ 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 =
+  Array.mapi (fun i shell_a ->
+    Array.mapi (fun j shell_b -> 
+      make ~cutoff shell_a shell_b)
+    (Array.sub basis 0 (i+1)) 
+  ) basis
+
+
+

Basis/AtomicShellPair.mli

@@ -0,0 +1,51 @@
+(** A datastructure to represent pairs of atomic shells.  
+*)
+
+type t
+
+
+val make : ?cutoff:float -> AtomicShell.t -> AtomicShell.t -> t option
+(** Creates an atomic shell pair from two atomic shells.
+
+    The contracted shell pairs contains the only pairs of primitives for which
+    the norm is greater than [cutoff]. 
+
+    If all the contracted shell pairs are not significant, the function returns
+    [None].
+*)
+
+val of_atomic_shell_array : AtomicShell.t array -> t option list
+(** 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].
+*)
+
+val atomic_shell_a : t -> AtomicShell.t
+(** Returns the first {!AtomicShell.t} which was used to build the atomic
+    shell pair.
+*)
+
+val atomic_shell_b : t -> AtomicShell.t
+(** Returns the second {!AtomicShell.t} which was used to build the atomic
+    shell pair.
+*)
+
+val contracted_shell_pairs : t -> ContractedShellPair.t list
+(** Returns an array of {!ContractedShellPair.t}, containing all the pairs of
+    contracted functions used to build the atomic shell pair.
+*)
+
+val center_ab : t -> Coordinate.t
+  (* A-B *)
+
+val norm_sq : t -> float
+  (* |A-B|^2 *)
+
+val norm_scales : t -> float array
+  (* norm_coef.(i) / norm_coef.(0) *)
+
+val totAngMom : t -> AngularMomentum.t
+  (* Total angular Momentum *)
+
+val monocentric : t -> bool
+(** If true, the two atomic shells have the same center. *)
+

Basis/ContractedShellPair.ml

@@ -7,7 +7,7 @@ type t =
 {
   shell_a         : ContractedShell.t;
   shell_b         : ContractedShell.t;
-  shell_pairs     : PrimitiveShellPair.t array;
+  shell_pairs     : PrimitiveShellPair.t list;
   coefficients    : float array; 
   exponents_inv   : float array; 
   center_ab       : Coordinate.t;       (* A-B *)
@@ -46,18 +46,17 @@ let make ?(cutoff=1.e-32) s_a s_b =
     |> Array.to_list
     |> List.filter (function (_, Some _) -> true | _ -> false)
     |> List.map (function (c, Some x) -> (c,x) | _ -> assert false)
-    |> Array.of_list
   in
 
 
-  let coefficients = Array.map (fun (c,y) -> c *. Psp.normalization y) shell_pairs
+  match shell_pairs with
-  and exponents_inv = Array.map (fun (_,y) -> Psp.exponent_inv y) shell_pairs
+  | [] -> None
+  | head :: _ -> 
+    let coefficients = List.map (fun (c,y) -> c *. Psp.normalization y) shell_pairs |> Array.of_list
+    and exponents_inv = List.map (fun (_,y) -> Psp.exponent_inv y) shell_pairs |> Array.of_list
     in
-  let shell_pairs = Array.map snd shell_pairs in
+    let shell_pairs = List.map snd shell_pairs in
-  if Array.length shell_pairs = 0 then
+    let root = snd head in
-    None
-  else
-    let root = shell_pairs.(0) in
     Some {
       shell_a = s_a ; shell_b = s_b ; coefficients ; exponents_inv ; shell_pairs ;
       center_ab = Psp.a_minus_b root;
@@ -69,7 +68,7 @@ let make ?(cutoff=1.e-32) s_a s_b =
 
 let shell_a         x = x.shell_a
 let shell_b         x = x.shell_b
-let shell_pairs     x = x.shell_pairs
+let shell_pairs     x = Array.of_list x.shell_pairs
 let coefficients    x = x.coefficients
 let exponents_inv   x = x.exponents_inv
 let center_ab       x = x.center_ab
@@ -77,7 +76,7 @@ let norm_sq         x = x.norm_sq
 let totAngMom       x = x.totAngMom
 let norm_scales     x = x.norm_scales
 
-let monocentric x = Psp.monocentric x.shell_pairs.(0)
+let monocentric x = Psp.monocentric (List.hd x.shell_pairs)
 
 (** Returns an integer characteristic of a contracted shell pair *)
 let hash a =