ContractedShell has now a hidden type

Basis/Basis.ml

@@ -27,14 +27,14 @@ let of_nuclei_and_general_basis n b =
   Array.iteri (fun i x ->
       if (i > 0) then
         result.(i) <- Cs.with_index x (
-            result.(i-1).index  +
+            Cs.index result.(i-1)  +
-            (Array.length result.(i-1).norm_coef_scale )) 
+            Cs.size_of_shell result.(i-1) ) 
     ) result ;
 
   let size =
     let n = ref 0 in
     for i=0 to (Array.length result) - 1 do
-      n := !n + (Array.length (result.(i).norm_coef_scale))
+      n := !n + Cs.size_of_shell result.(i)
     done; !n
   in
   { contracted_shells = result ; size }

Basis/ContractedShell.ml

@@ -3,14 +3,14 @@ open Constants
 open Coordinate
 
 type t = {
-  expo      : float array;
+  expo      : float array;       (** Array of exponents {% $\alpha_i$ %} *)
-  coef      : float array;
+  coef      : float array;       (** Array of contraction coefficients {% $d_i$ %} *)
-  center    : Coordinate.t;
+  center    : Coordinate.t;      (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %} *)
-  totAngMom : AngularMomentum.t;
+  totAngMom : AngularMomentum.t; (** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *)
-  size      : int;
+  size      : int;               (** Number of contracted functions, {% $m$ %} in the formula *)
-  norm_coef : float array;
+  norm_coef : float array;       (** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *)
-  norm_coef_scale : float array;
+  norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *)
-  index    : int;
+  index    : int;                (** Index in the basis set, represented as an array of contracted shells. *)
 }
 
 module Am = AngularMomentum
@@ -113,4 +113,20 @@ let compute_norm_coef expo totAngMom =
   Array.map (fun x -> let f a = x *. factor a in f) expo
 
 
+let expo x = x.expo
 
+let coef x = x.coef
+
+let center x = x.center
+
+let totAngMom x = x.totAngMom
+
+let size x = Array.length x.coef
+
+let norm_coef x = x.norm_coef
+
+let norm_coef_scale x = x.norm_coef_scale
+
+let index x = x.index
+
+let size_of_shell x = Array.length x.norm_coef_scale

Basis/ContractedShell.mli

@@ -28,28 +28,46 @@ f_i = \frac{1}{\mathcal{N}_i}
 
 *)
 
-type t = private {
+type t 
-  expo      : float array;       (** Array of exponents {% $\alpha_i$ %} *)
-  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. *)
-}
 
-
-
-(** Pretty-printing of the contracted shell in a string *)
 val to_string : t -> string
+(** Pretty-printing of the contracted shell in a string *)
 
-(** Creates a contracted shell *)
 val make :
   index:int ->
   expo:float array ->
   coef:float array ->
   center:Coordinate.t -> totAngMom:AngularMomentum.t -> t
+(** Creates a contracted shell *)
 
-(** Returns a copy of the contracted shell with a modified index *)
 val with_index  : t -> int -> t
+(** Returns a copy of the contracted shell with a modified index *)
+
+
+val expo : t -> float array
+(** Array of exponents {% $\alpha_i$ %} *)
+
+val coef : t -> float array
+(** Array of contraction coefficients {% $d_i$ %} *)
+
+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 size : t -> int
+(** Number of contracted functions, {% $m$ %} in the formula *)
+
+val norm_coef : t -> float array
+(** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *)
+
+val norm_coef_scale : t -> float array
+(** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *)
+
+val index : t -> int
+(** Index in the basis set, represented as an array of contracted shells. *)
+
+val size_of_shell : t -> int
+(** Number of contracted functions in the shell *)
+

Basis/ContractedShellPair.ml

@@ -34,15 +34,15 @@ let create ?cutoff p_a p_b =
     | Some cutoff -> cutoff, -. (log cutoff)
   in
 
-  let center_ab = Co.( p_a.Cs.center |- p_b.Cs.center )
+  let center_ab = Co.( Cs.center p_a |- Cs.center p_b )
   in
   let norm_sq =
     Co.dot center_ab center_ab
   in
   let norm_coef_scale_a =
-    p_a.norm_coef_scale 
+    Cs. norm_coef_scale p_a
   and norm_coef_scale_b =
-    p_b.norm_coef_scale 
+    Cs. norm_coef_scale p_b 
   in
   let norm_coef_scale =
     Array.map (fun v1 ->
@@ -52,24 +52,24 @@ let create ?cutoff p_a p_b =
     |> Array.concat
   in
   let shell_pairs = 
-    Array.init p_a.size (fun i ->
+    Array.init (Cs.size p_a) (fun i ->
-        let p_a_expo_center = Co.(p_a.Cs.expo.(i) |. p_a.Cs.center) in
+        let p_a_expo_center = Co.( (Cs.expo p_a).(i) |. Cs.center p_a ) in
-        let norm_coef_a = p_a.norm_coef.(i)  in
+        let norm_coef_a = (Cs.norm_coef p_a).(i)  in
 
-        Array.init p_b.size (fun j ->
+        Array.init (Cs.size p_b) (fun j ->
             try
-              let norm_coef_b = p_b.norm_coef.(j) in
+              let norm_coef_b = (Cs.norm_coef p_b).(j) in
               let norm_coef   = norm_coef_a *. norm_coef_b
               in
               if norm_coef < cutoff  then
                 raise Null_contribution;
-              let p_b_expo_center = Co.(p_b.expo.(j) |. p_b.center) in
+              let p_b_expo_center = Co.( (Cs.expo p_b).(j) |. Cs.center p_b ) in
-              let expo = p_a.expo.(i) +. p_b.expo.(j) in
+              let expo = (Cs.expo p_a).(i) +. (Cs.expo p_b).(j) in
               let expo_inv = 1. /. expo in
               let center = Co.(expo_inv |. (p_a_expo_center |+ p_b_expo_center ) )
               in
               let argexpo =
-                p_a.Cs.expo.(i) *. p_b.Cs.expo.(j) *. norm_sq *. expo_inv
+                (Cs.expo p_a).(i) *. (Cs.expo p_b).(j) *. norm_sq *. expo_inv
               in
               if (argexpo > log_cutoff) then
                 raise Null_contribution;
@@ -77,19 +77,19 @@ let create ?cutoff p_a p_b =
                 (pi *. expo_inv)**(1.5) *. exp (-. argexpo)
               in
               let coef =
-                norm_coef *. p_a.coef.(i) *. p_b.coef.(j) *. g
+                norm_coef *. (Cs.coef p_a).(i) *. (Cs.coef p_b).(j) *. g
               in
               if abs_float coef  < cutoff then
                 raise Null_contribution;
               let center_a =
-                Co.(center |- p_a.center)
+                Co.(center |- Cs.center p_a)
               in
               let monocentric =
-                p_a.center  = p_b.center 
+                Cs.(center p_a = center p_b)
               in
               let totAngMomInt = 
-                Am.to_int p_a.totAngMom +
+                Am.to_int (Cs.totAngMom p_a) +
-                Am.to_int p_b.totAngMom 
+                Am.to_int (Cs.totAngMom p_b)
               in
               Some {
                 Sp.i ; j ;

Basis/ERI.ml

@@ -108,7 +108,7 @@ let of_basis basis =
 
   let icount = ref 0 in
   for i=0 to (Array.length shell) - 1 do
-    print_int shell.(i).index ; print_newline ();
+    print_int (Cs.index shell.(i)) ; print_newline ();
     for j=0 to i do
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
         if (schwartz_p_max >= cutoff) then
@@ -132,7 +132,7 @@ let of_basis basis =
   let inn = ref 0 and out = ref 0 in
 
   for i=0 to (Array.length shell) - 1 do
-    print_int shell.(i).index ; print_newline ();
+    print_int (Cs.index shell.(i)) ; print_newline ();
     for j=0 to i do
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
       try
@@ -180,16 +180,16 @@ let of_basis basis =
 
             (* Write the data in the output file *)
             Array.iteri (fun i_c powers_i ->
-                let i_c = shell.(i).index + i_c + 1 in
+                let i_c = Cs.index shell.(i) + i_c + 1 in
                 let xi = to_powers powers_i in
                 Array.iteri (fun j_c powers_j ->
-                    let j_c = shell.(j).index + j_c + 1 in
+                    let j_c = Cs.index shell.(j) + j_c + 1 in
                     let xj = to_powers powers_j in
                     Array.iteri (fun k_c powers_k ->
-                        let k_c = shell.(k).index + k_c + 1 in 
+                        let k_c = Cs.index shell.(k) + k_c + 1 in 
                         let xk = to_powers powers_k in
                         Array.iteri (fun l_c powers_l ->
-                            let l_c = shell.(l).index + l_c + 1 in
+                            let l_c = Cs.index shell.(l) + l_c + 1 in
                             let xl = to_powers powers_l in
                             let key = 
                               if swap then
@@ -213,10 +213,10 @@ let of_basis basis =
                               )
                             else
                               out := !out + 1;
-                          ) Am.(zkey_array (Singlet shell.(l).Cs.totAngMom))
+                          ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(l))))
-                      ) Am.(zkey_array (Singlet shell.(k).Cs.totAngMom))
+                      ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(k))))
-                  ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom))
+                  ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(j))))
-              ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom))
+              ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(i))))
             with NullIntegral ->  ()
           done;
         done;

Basis/KinInt.ml

@@ -29,7 +29,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
-    Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom))
+    Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b))
   in
 
   let contracted_class = 
@@ -64,9 +64,9 @@ let contracted_class shell_a shell_b : float Zmap.t =
           sp.(ab).Sp.i, sp.(ab).Sp.j 
         in
         let expo_a = 
-          sp.(ab).Sp.shell_a.expo.(i)
+          (Cs.expo sp.(ab).Sp.shell_a).(i)
         and expo_b = 
-          sp.(ab).Sp.shell_b.expo.(j)
+          (Cs.expo sp.(ab).Sp.shell_b).(j)
         in
 
         let xyz_of_int k =
@@ -141,10 +141,10 @@ let of_basis basis =
       in
 
       Array.iteri (fun j_c powers_j ->
-          let j_c = shell.(j).index + j_c + 1 in
+          let j_c = Cs.index shell.(j) + j_c + 1 in
           let xj = to_powers powers_j in
           Array.iteri (fun i_c powers_i ->
-              let i_c = shell.(i).index + i_c + 1 in
+              let i_c = Cs.index shell.(i) + i_c + 1 in
               let xi = to_powers powers_i in
               let key =
                 Zkey.of_powers_six xi xj
@@ -155,8 +155,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom))
+            ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
-        ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) 
+        ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) 
     done;
   done;
   Mat.detri result;

Basis/NucInt.ml

@@ -75,10 +75,10 @@ let of_basis_nuclei basis nuclei =
 
       (* Write the data in the output file *)
       Array.iteri (fun i_c powers_i ->
-          let i_c = shell.(i).index + i_c + 1 in
+          let i_c = Cs.index shell.(i) + i_c + 1 in
           let xi = to_powers powers_i in
           Array.iteri (fun j_c powers_j ->
-              let j_c = shell.(j).index + j_c + 1 in
+              let j_c = Cs.index shell.(j) + j_c + 1 in
               let xj = to_powers powers_j in
               let key = 
                 Zkey.of_powers_six xi xj
@@ -88,8 +88,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 shell.(j).Cs.totAngMom))
+            ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j))))
-        ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) 
+        ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
     done;
   done;
   Mat.detri eni_array;

Basis/OneElectronRR.ml

@@ -101,12 +101,12 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
   and shell_b = shell_p.Csp.shell_b
   in
   let maxm = 
-    (Am.to_int @@ shell_a.totAngMom) + (Am.to_int @@ shell_b.totAngMom) 
+    Am.to_int (Cs.totAngMom shell_a) + Am.to_int (Cs.totAngMom shell_b) 
   in
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
-    Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom))
+    Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b))
   in
 
   let contracted_class = 
@@ -140,7 +140,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
           shell_p.Csp.shell_pairs.(ab).Sp.center
         in
         let center_pa = 
-          Co.(center_p |- shell_a.center) 
+          Co.(center_p |- Cs.center shell_a) 
         in
 
         for c=0 to Array.length geometry - 1 do
@@ -156,7 +156,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 (shell_a.totAngMom, shell_b.totAngMom) with
+          match Cs.(totAngMom shell_a, totAngMom shell_b) with
           | Am.(S,S) -> 
             let integral = zero_m_array.(0) in
             contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge
@@ -177,7 +177,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
                   hvrr_one_e
                     angMomA  angMomB
                     zero_m_array
-                    shell_b.expo.(b)
+                    (Cs.expo shell_b).(b)
                     shell_p.Csp.expo_inv.(ab)
                     center_ab center_pa center_pc
                     map

Basis/Overlap.ml

@@ -29,7 +29,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
 
   (* Pre-computation of integral class indices *)
   let class_indices = 
-    Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom))
+    Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b))
   in
 
   let contracted_class = 
@@ -113,10 +113,10 @@ let of_basis basis =
       in
 
       Array.iteri (fun j_c powers_j ->
-          let j_c = shell.(j).index + j_c + 1 in
+          let j_c = Cs.index shell.(j) + j_c + 1 in
           let xj = to_powers powers_j in
           Array.iteri (fun i_c powers_i ->
-              let i_c = shell.(i).index + i_c + 1 in
+              let i_c = Cs.index shell.(i) + i_c + 1 in
               let xi = to_powers powers_i in
               let key = 
                 Zkey.of_powers_six xi xj
@@ -127,8 +127,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom))
+            ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i))))
-        ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) 
+        ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j))))
     done;
   done;
   Mat.detri result;

Basis/TwoElectronRR.ml

@@ -286,8 +286,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   (* Pre-computation of integral class indices *)
   let class_indices = 
     Am.zkey_array (Am.Quartet
-         (shell_a.totAngMom, shell_b.totAngMom,
+         Cs.(totAngMom shell_a, totAngMom shell_b,
-          shell_c.totAngMom, shell_d.totAngMom ))
+             totAngMom shell_c, totAngMom shell_d ))
   in
 
   let contracted_class = 
@@ -297,8 +297,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   let monocentric =
       shell_p.Csp.monocentric &&
       shell_q.Csp.monocentric &&
-      shell_p.Csp.shell_a.Cs.center  = 
+      Cs.center shell_p.Csp.shell_a  = 
-      shell_q.Csp.shell_a.Cs.center 
+      Cs.center shell_q.Csp.shell_a 
   in
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
@@ -333,8 +333,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
         in
         begin
-          match (shell_a.totAngMom, shell_b.totAngMom,
+          match Cs.(totAngMom shell_a, totAngMom shell_b,
-                 shell_c.totAngMom, shell_d.totAngMom) with
+                    totAngMom shell_c, totAngMom shell_d) with
           | Am.(S,S,S,S) ->
             let integral =
               zero_m_array.(0)
@@ -402,7 +402,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                     hvrr_two_e
                       angMom_a angMom_b angMom_c angMom_d
                       zero_m_array
-                      shell_b.expo.(b) shell_d.expo.(d)
+                      (Cs.expo shell_b).(b) (Cs.expo shell_d).(d)
                       shell_p.Csp.expo_inv.(ab)
                       shell_q.Csp.expo_inv.(cd)
                       sp.(ab).Sp.center_ab sq.(cd).Sp.center_ab center_pq

Basis/TwoElectronRRVectorized.ml

@@ -4,6 +4,7 @@ open Bigarray
 
 module Am  = AngularMomentum
 module Co  = Coordinate
+module Cs  = ContractedShell
 module Csp = ContractedShellPair
 module Sp  = ShellPair
 module Po  = Powers
@@ -568,8 +569,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   (* Pre-computation of integral class indices *)
   let class_indices = 
     Am.zkey_array (Am.Quartet
-      (shell_a.totAngMom, shell_b.totAngMom,
+      Cs.(totAngMom shell_a, totAngMom shell_b,
-       shell_c.totAngMom, shell_d.totAngMom))
+          totAngMom shell_c, totAngMom shell_d))
   in
 
   let contracted_class = 
@@ -579,8 +580,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   let monocentric =
       shell_p.Csp.monocentric &&
       shell_q.Csp.monocentric &&
-      shell_p.Csp.shell_a.center =
+      Cs.center shell_p.Csp.shell_a =
-      shell_q.Csp.shell_a.center 
+      Cs.center shell_q.Csp.shell_a
   in
 
   (** Screening on the product of coefficients *)
@@ -623,8 +624,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
   begin
-    match (shell_a.totAngMom, shell_b.totAngMom,
+    match Cs.(totAngMom shell_a, totAngMom shell_b,
-           shell_c.totAngMom, shell_d.totAngMom) with
+              totAngMom shell_c, totAngMom shell_d) with
     | Am.(S,S,S,S) ->
       contracted_class.(0) <-
         begin
@@ -684,9 +685,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
       in
 
       let expo_b =
-        Array.map (fun shell_ab -> shell_b.expo.(shell_ab.Sp.j)) sp
+        Array.map (fun shell_ab -> (Cs.expo shell_b).(shell_ab.Sp.j)) sp
       and expo_d =
-        Array.map (fun shell_cd -> shell_d.expo.(shell_cd.Sp.j)) sq
+        Array.map (fun shell_cd -> (Cs.expo shell_d).(shell_cd.Sp.j)) sq
       in
       let norm_coef_scale_p = shell_p.Csp.norm_coef_scale in
 
@@ -719,7 +720,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           Array.init np (fun ab ->
             let shell_ab = sp.(ab) in
                 let cpa = 
-                  Co.(shell_ab.Sp.center |- shell_a.center)
+                  Co.(shell_ab.Sp.center |- Cs.center shell_a)
                 in
                 match xyz with
                 | 0 -> Co.(get X cpa);
@@ -738,7 +739,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           Array.init nq (fun cd ->
             let shell_cd = sq.(cd) in
                 let cqc = 
-                  Co.(shell_cd.Sp.center |- shell_c.center)
+                  Co.(shell_cd.Sp.center |- Cs.center shell_c)
                 in
                 match xyz with
                 | 0 -> Co.(get X cqc);