Hash table for zkey_array

Basis/Basis.ml

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

Basis/ContractedShell.ml

@@ -10,19 +10,12 @@ type t = {
   size      : int;
   norm_coef : float array;
   norm_coef_scale : float array;
-  powers   : Zkey.t array;
   index    : int;
 }
 
 module Am = AngularMomentum
 
-(** 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
@@ -67,7 +60,7 @@ let make ~index ~expo ~coef ~center ~totAngMom =
     ) powers
   in
   { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
-    norm_coef_scale ; powers }
+    norm_coef_scale }
 
 
 let with_index a i =
@@ -79,7 +72,7 @@ let to_string s =
   let open Printf in
   (match s.totAngMom with
    | Am.S -> sprintf "%3d    " (s.index+1)
-   | _                  -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.powers))
+   | _                  -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.norm_coef_scale))
   ) ^
   ( sprintf "%1s %8.3f %8.3f %8.3f " (Am.to_string s.totAngMom)
       (get X coord) (get Y coord) (get Z coord) ) ^

Basis/ContractedShell.mli

@@ -35,8 +35,7 @@ type t = private {
   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 [powers]. *)
+  norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *)
-  powers   : Zkey.t array;       (** Triplets {% $(n_x,n_y,n_z)$ %} encoded in a {!Zkey.t}. *)
   index    : int;                (** Index in the basis set, represented as an array of contracted shells. *)
 }
 

Basis/ERI.ml

@@ -6,6 +6,7 @@ open Bigarray
 
 type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t
 
+module Am  = AngularMomentum
 module Bs  = Basis
 module Cs  = ContractedShell
 module Csp = ContractedShellPair
@@ -212,10 +213,10 @@ let of_basis basis =
                               )
                             else
                               out := !out + 1;
-                          ) shell.(l).powers 
+                          ) Am.(zkey_array (Singlet shell.(l).Cs.totAngMom))
-                      ) shell.(k).powers
+                      ) Am.(zkey_array (Singlet shell.(k).Cs.totAngMom))
-                  ) shell.(j).powers
+                  ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom))
-              ) shell.(i).powers
+              ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom))
             with NullIntegral ->  ()
           done;
         done;

Basis/GamessReader.ml

@@ -1,3 +1,4 @@
+(** Reads a basis set in GAMESS format *)
 let read_basis filename =
   let lexbuf =
     let ic = open_in filename in

Basis/KinInt.ml

@@ -155,8 +155,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) shell.(i).powers
+            ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom))
-        ) shell.(j).powers 
+        ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) 
     done;
   done;
   Mat.detri result;

Basis/NucInt.ml

@@ -6,7 +6,9 @@ open Lacaml.D
 
 type t = Mat.t
 
+module Am = AngularMomentum
 module Bs = Basis
+module Cs = ContractedShell
 
 (** (0|0)^m : Fundamental electron-nucleus repulsion integral
     $ \int \phi_p(r1) 1/r_{C} dr_1 $
@@ -86,8 +88,8 @@ let of_basis_nuclei basis nuclei =
               in
               eni_array.{j_c,i_c} <- value;
               eni_array.{i_c,j_c} <- value;
-            ) shell.(j).powers
+            ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom))
-        ) shell.(i).powers 
+        ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) 
     done;
   done;
   Mat.detri eni_array;

Basis/Overlap.ml

@@ -127,8 +127,8 @@ let of_basis basis =
               in
               result.{i_c,j_c} <- value;
               result.{j_c,i_c} <- value;
-            ) shell.(i).powers
+            ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom))
-        ) shell.(j).powers
+        ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) 
     done;
   done;
   Mat.detri result;

Basis/Overlap_primitives.ml

@@ -8,22 +8,23 @@ let hvrr (angMom_a, angMom_b) (expo_inv_p) (center_ab, center_pa)
   (** Vertical recurrence relations *)
   let rec vrr angMom_a =
 
-    if angMom_a < 0 then 0.
+    if angMom_a > 0 then
-    else if angMom_a = 0 then 1.
+      chop center_pa (fun () -> vrr (angMom_a-1))
-    else
+      +. chop (0.5 *. float_of_int (angMom_a-1) *. expo_inv_p) 
-      (chop center_pa (fun () -> vrr (angMom_a-1)))
-      +. chop (0.5 *. (float_of_int (angMom_a-1)) *. expo_inv_p) 
         (fun () -> vrr (angMom_a-2))
+    else if angMom_a = 0 then 1.
+    else 0.
 
 
   (** Horizontal recurrence relations *)
   and hrr angMom_a angMom_b =
 
-    if angMom_b < 0 then 0.
+    if angMom_b > 0 then
-    else if angMom_b = 0 then vrr angMom_a 
-    else
       hrr (angMom_a+1) (angMom_b-1) +. chop center_ab
       (fun () -> hrr angMom_a (angMom_b-1) )
+    else if angMom_b = 0 then
+      vrr angMom_a 
+    else 0.
 
   in
   hrr angMom_a angMom_b 

HartreeFock/Fock.ml

@@ -1,5 +1,6 @@
 open Lacaml.D
 open Simulation
+open Constants
 
 type t = Mat.t
 
@@ -17,10 +18,11 @@ let make ~density simulation =
     for nu = 1 to nBas do
       for lambda = 1 to nBas do
         let p = m_P.{lambda,sigma} in
-        for mu = 1 to nu do
+        if abs_float p > epsilon then
-          m_F.{mu,nu} <- m_F.{mu,nu} +.  p *. 
+          for mu = 1 to nu do
-                        (m_G.{mu,lambda,nu,sigma} -. 0.5 *. m_G.{mu,lambda,sigma,nu})
+            m_F.{mu,nu} <- m_F.{mu,nu} +.  p *. 
-        done
+                          (m_G.{mu,lambda,nu,sigma} -. 0.5 *. m_G.{mu,lambda,sigma,nu})
+          done
       done
     done
   done;

Utils/AngularMomentum.ml

@@ -2,7 +2,7 @@ open Powers
 
 exception AngularMomentumError of string
 
-type t =
+type t = 
   | S | P | D | F | G | H | I | J | K | L | M | N | O
 
 let of_char = function
@@ -68,8 +68,12 @@ let n_functions a =
   (a*a + 3*a + 2)/2
 
 
+let zkey_array_memo : (kind, Zkey.t array) Hashtbl.t =
+  Hashtbl.create 13 
+
 (** Returns an array of Zkeys corresponding to all possible angular momenta *)
 let zkey_array a = 
+
   let keys_1d l =
     let create_z { x ; y ; _ } =
       Powers.of_int_tuple (x,y,l-(x+y))
@@ -93,41 +97,50 @@ let zkey_array a =
     |> List.rev
   in
 
-  begin
+  try
-    match a with
+    Hashtbl.find zkey_array_memo a
-    | Singlet l1 -> 
-      List.map (fun x -> Zkey.of_powers_three x) (keys_1d @@ to_int l1)
 
-    | Doublet (l1, l2) -> 
+  with Not_found ->
-      List.map (fun a ->
-          List.map (fun b -> Zkey.of_powers_six a b) (keys_1d @@ to_int l2)
-        ) (keys_1d @@ to_int l1)
-      |> List.concat
 
-    | Triplet (l1, l2, l3) ->
+    let result = 
+      begin
+        match a with
+        | Singlet l1 -> 
+          List.map (fun x -> Zkey.of_powers_three x) (keys_1d @@ to_int l1)
 
-      List.map (fun a ->
+        | Doublet (l1, l2) -> 
-          List.map (fun b ->
+          List.map (fun a ->
-              List.map (fun c ->
+              List.map (fun b -> Zkey.of_powers_six a b) (keys_1d @@ to_int l2)
-                  Zkey.of_powers_nine a b c) (keys_1d @@ to_int l3)
+            ) (keys_1d @@ to_int l1)
-            ) (keys_1d @@ to_int l2)
           |> List.concat
-        ) (keys_1d @@ to_int l1)
-      |> List.concat
 
-    | Quartet (l1, l2, l3, l4) ->
+        | Triplet (l1, l2, l3) ->
 
-      List.map (fun a ->
+          List.map (fun a ->
-          List.map (fun b ->
+              List.map (fun b ->
-              List.map (fun c ->
+                  List.map (fun c ->
-                  List.map (fun d ->
+                      Zkey.of_powers_nine a b c) (keys_1d @@ to_int l3)
-                      Zkey.of_powers_twelve a b c d) (keys_1d @@ to_int l4)
+                ) (keys_1d @@ to_int l2)
-                ) (keys_1d @@ to_int l3)
               |> List.concat
-            ) (keys_1d @@ to_int l2)
+            ) (keys_1d @@ to_int l1)
           |> List.concat
-        ) (keys_1d @@ to_int l1)
+
-        |> List.concat
+        | Quartet (l1, l2, l3, l4) ->
-  end
+
-  |> Array.of_list
+          List.map (fun a ->
+              List.map (fun b ->
+                  List.map (fun c ->
+                      List.map (fun d ->
+                          Zkey.of_powers_twelve a b c d) (keys_1d @@ to_int l4)
+                    ) (keys_1d @@ to_int l3)
+                  |> List.concat
+                ) (keys_1d @@ to_int l2)
+              |> List.concat
+            ) (keys_1d @@ to_int l1)
+          |> List.concat
+      end
+      |> Array.of_list
+    in
+    Hashtbl.add zkey_array_memo a result;
+    result
 

Utils/NonNegativeFloat.ml

@@ -1,4 +1,4 @@
-type t = float
+type t = float 
 
 let of_float x =
   if x < 0. then invalid_arg (__FILE__^": of_float");

Utils/Zkey.ml

@@ -2,8 +2,8 @@ open Powers
 
 type t = 
   {
-    mutable left   : int ;
+    mutable left   : int;
-    mutable right  : int ;
+    mutable right  : int;
     kind   : int ;
   }
 

Utils/Zkey.mli

@@ -67,3 +67,4 @@ val equal : t -> t -> bool
 val compare : t -> t -> int
 (** Comparison function, used for sorting. *)
 
+