Primitive Shell

Basis/Basis.ml

@@ -6,6 +6,7 @@ type t =
 
 module Cs = ContractedShell
 module Gb = GeneralBasis
+module Ps = PrimitiveShell
 
 
 (** Returns an array of the basis set per atom *)
@@ -14,12 +15,11 @@ let of_nuclei_and_general_basis n b =
     Array.map (fun (e, center) ->
         List.assoc e b
         |> Array.map (fun (totAngMom, shell) ->
-            let expo = Array.map (fun Gb.{exponent ; coefficient} ->
-                exponent) shell
-            and coef = Array.map (fun Gb.{exponent ; coefficient} ->
-                coefficient) shell
+            let lc = 
+              Array.map (fun Gb.{exponent ; coefficient} ->
+                coefficient, Ps.make totAngMom center exponent) shell
             in
-            Cs.make ~expo ~coef ~totAngMom ~center ~index:0)
+            Cs.make lc)
       ) n
     |> Array.to_list
     |> Array.concat

Basis/ContractedShell.ml

@@ -7,59 +7,46 @@ type t = {
   coef      : float array;       (** Array of contraction coefficients {% $d_i$ %} *)
   center    : Coordinate.t;      (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %} *)
   totAngMom : AngularMomentum.t; (** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *)
-  size      : int;               (** Number of contracted functions, {% $m$ %} in the formula *)
   norm_coef : float array;       (** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *)
   norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *)
   index    : int;                (** Index in the basis set, represented as an array of contracted shells. *)
 }
 
 module Am = AngularMomentum
+module Ps = PrimitiveShell
 
 
-let compute_norm_coef expo totAngMom =
-  let atot =
-    Am.to_int totAngMom
-  in
-  let factor int_array = 
-       let dfa = Array.map (fun j ->
-           ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
-         ) int_array
-       in
-       sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
-  in
-  let expo = 
-    if atot mod 2 = 0 then
-      Array.map (fun alpha ->
-        let alpha_2 = alpha +. alpha in
-        (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
-      ) expo
-    else
-      Array.map (fun alpha ->
-        let alpha_2 = alpha +. alpha in
-        (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
-      ) expo
-  in
-  Array.map (fun x -> let f a = x *. (factor a) in f) expo
+let make ?(index=0) lc = 
+  assert (Array.length lc > 0);
 
+  let coef = Array.map fst lc
+  and prim = Array.map snd lc
+  in
 
-let make ~index ~expo ~coef ~center ~totAngMom =
-  assert (Array.length expo = Array.length coef);
-  assert (Array.length expo > 0);
-  let norm_coef_func =
-    compute_norm_coef expo totAngMom
+  let center = Ps.center prim.(0) in
+  let rec unique_center = function
+  | 0 -> true
+  | i -> if Ps.center prim.(i) = center then unique_center (i-1) else false
   in
-  let powers =
-    Am.zkey_array (Am.Singlet totAngMom) 
+  if not (unique_center (Array.length prim - 1)) then
+    invalid_arg "ContractedShell.make Coordinate.t differ";
+    
+  let totAngMom = Ps.totAngMom prim.(0) in
+  let rec unique_angmom = function
+  | 0 -> true
+  | i -> if Ps.totAngMom prim.(i) = totAngMom then unique_angmom (i-1) else false
   in
+  if not (unique_angmom (Array.length prim - 1)) then
+    invalid_arg "ContractedShell.make: AngularMomentum.t differ";
+    
+  let expo = Array.map Ps.expo prim in
+
   let norm_coef =
-    Array.map (fun f -> f [| Am.to_int totAngMom ; 0 ; 0 |])  norm_coef_func
+    Array.map Ps.norm_coef prim
   in
-  let norm_coef_scale = 
-    Array.map (fun a ->
-      (norm_coef_func.(0) (Zkey.to_int_array a)) /. norm_coef.(0)
-    ) powers
+  let norm_coef_scale = Ps.norm_coef_scale prim.(0)
   in
-  { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
+  { index ; expo ; coef ; center ; totAngMom ; norm_coef ;
     norm_coef_scale }
 
 
@@ -80,38 +67,6 @@ let to_string s =
    |> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) 
 
 
-(** Normalization coefficient of contracted function i, which depends on the
-    exponent and the angular momentum. Two conventions can be chosen : a single
-    normalisation factor for all functions of the class, or a coefficient which
-    depends on the powers of x,y and z.
-    Returns, for each contracted function, an array of functions taking as
-    argument the [|x;y;z|] powers.
-*)
-let compute_norm_coef expo totAngMom =
-  let atot =
-    Am.to_int totAngMom
-  in
-  let factor int_array = 
-       let dfa = Array.map (fun j ->
-           (float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
-         ) int_array
-       in
-       sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
-  in
-  let expo = 
-    if atot mod 2 = 0 then
-      Array.map (fun alpha ->
-        let alpha_2 = alpha +. alpha in
-        (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
-      ) expo
-    else
-      Array.map (fun alpha ->
-        let alpha_2 = alpha +. alpha in
-        (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
-      ) expo
-  in
-  Array.map (fun x -> let f a = x *. factor a in f) expo
-
 
 let expo x = x.expo
 

Basis/ContractedShell.mli

@@ -33,12 +33,8 @@ type t
 val to_string : t -> string
 (** Pretty-printing of the contracted shell in a string *)
 
-val make :
-  index:int ->
-  expo:float array ->
-  coef:float array ->
-  center:Coordinate.t -> totAngMom:AngularMomentum.t -> t
-(** Creates a contracted shell *)
+val make : ?index:int -> (float * PrimitiveShell.t) array -> t 
+(** 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 *)

Basis/ContractedShellPair.mli

@@ -61,8 +61,10 @@ val totAngMomInt : t -> int
 val monocentric : t -> bool
 (** If true, the two contracted shells have the same center. *)
 
-val hash : 'a array -> int array
-val cmp : 'a array -> 'a array -> int
-val equivalent : 'a array -> 'a array -> bool
+(*
+val hash : t -> int array
+val cmp : t -> t -> int
+val equivalent : t -> t -> bool
 val unique : 'a array array array -> 'a array list
 val indices : 'a array array array -> (int array, int * int) Hashtbl.t
+*)

Basis/PrimitiveShell.ml

@@ -0,0 +1,82 @@
+open Util
+open Constants
+open Coordinate
+
+type t = {
+  expo      : float;
+  norm_coef : float;
+  norm_coef_scale : float array lazy_t;
+  center    : Coordinate.t;
+  totAngMom : AngularMomentum.t;
+}
+
+module Am = AngularMomentum
+
+
+let compute_norm_coef alpha totAngMom =
+  let atot =
+    Am.to_int totAngMom
+  in
+  let factor int_array = 
+       let dfa = Array.map (fun j ->
+           ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
+         ) int_array
+       in
+       sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
+  in
+  let expo = 
+    if atot mod 2 = 0 then
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
+    else
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
+  in
+  let f a =
+    expo *. (factor a)
+  in f
+
+
+let make totAngMom center expo =
+  let norm_coef_func =
+    compute_norm_coef expo totAngMom
+  in
+  let norm_coef =
+    norm_coef_func [| Am.to_int totAngMom ; 0 ; 0 |]  
+  in
+  let powers =
+    Am.zkey_array (Am.Singlet totAngMom) 
+  in
+  let norm_coef_scale = lazy (
+    Array.map (fun a ->
+      (norm_coef_func (Zkey.to_int_array a)) /. norm_coef
+    ) powers )
+  in
+  { expo ; norm_coef ; norm_coef_scale ; center ; totAngMom }
+
+
+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)
+      (get X coord) (get Y coord) (get Z coord) s.expo
+
+
+(** Normalization coefficient of contracted function i, which depends on the
+    exponent and the angular momentum. Two conventions can be chosen : a single
+    normalisation factor for all functions of the class, or a coefficient which
+    depends on the powers of x,y and z.
+    Returns, for each contracted function, an array of functions taking as
+    argument the [|x;y;z|] powers.
+*)
+
+let expo x = x.expo
+
+let center x = x.center
+
+let totAngMom x = x.totAngMom
+
+let norm_coef x = x.norm_coef
+
+let norm_coef_scale x = Lazy.force x.norm_coef_scale
+
+let size_of_shell x = Array.length (norm_coef_scale x)

Basis/PrimitiveShell.mli

@@ -0,0 +1,54 @@
+(** Set of Gaussians with a given {!AngularMomentum.t}
+
+{% \\[
+g(r) = (x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} \exp \left( -\alpha |r-R_A|^2 \right)
+   \\] %}
+
+where:
+
+- {% $n_x + n_y + n_z = l$ %}, the total angular momentum
+
+- {% $\alpha$ %} is the exponent
+
+*)
+
+type t 
+
+val to_string : t -> string
+(** Pretty-printing of the primitive shell in a string. *)
+
+val make : AngularMomentum.t -> Coordinate.t -> float -> t
+(** Creates a primitive shell from the total angular momentum, the coordinates of the
+    center and the exponent. *)
+
+val expo : t -> float 
+(** Returns the exponent {% $\alpha$ %}. *)
+
+val center : t -> Coordinate.t
+(** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *)
+
+val totAngMom : t -> AngularMomentum.t
+(** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *)
+
+val norm_coef : t -> float 
+(** Normalization coefficient of the shell:
+
+    {% \\[
+    \mathcal{N} = \sqrt{\iiint \left[ (x-X_A)^{l}
+                  \exp (-\alpha |r-R_A|^2) \right]^2 \, dx\, dy\, dz}
+      \\] %}
+*)
+
+val norm_coef_scale : t -> float array
+(** Scaling factors adjusting the normalization coefficient for the.
+  particular powers of {% $x,y,z$ %}.  They are given in the same order as
+  [AngularMomentum.zkey_array totAngMom]: 
+
+  {% \\[
+  f = \frac{1}{\mathcal{N}} \sqrt{\iiint [g(r)]^2 \, d^3r}
+    \\] %}
+*)
+
+val size_of_shell : t -> int
+(** Number of functions in the shell. *)
+