Working on contracted shell pairs

Basis/ContractedShellPair.ml

@@ -1,9 +1,13 @@
 open Util
 open Constants
 
+exception Null_contribution
+
+(** Array of shell pairs obtained from combining two contracted shells. The
+    two integers correspond to the indices of the combined shells.
+*)
 type t = ShellPair.t array
 
-exception Null_contribution
 
 (** Creates an contracted shell pair : an array of pairs of primitive shells.
     A contracted shell with N functions combined with a contracted
@@ -94,30 +98,66 @@ let create ?cutoff p_a p_b =
 
 (** Returns an integer characteristic of a contracted shell pair *)
 let hash a = 
-  1 (*TODO*)
+  Array.map Hashtbl.hash a
 
 
 (** Comparison function, used for sorting *)
 let cmp a b =
-  hash a - hash b
+  if a = b then 0
+  else if (Array.length a < Array.length b) then -1
+  else if (Array.length a > Array.length b) then  1
+  else 
+    let out = ref 0 in
+    begin
+     try
+        for k=0 to (Array.length a - 1) do
+          if a.(k) < b.(k) then
+            (out := (-1); raise Not_found)
+          else if a.(k) > b.(k) then
+            (out :=  1; raise Not_found);
+        done
+     with Not_found -> ();
+    end;
+    !out
 
 
-(** The array of all shell pairs *)
+(** The array of all shell pairs with their correspondance in the list
+    of contracted shells.
+*)
 let shell_pairs basis =
-    Array.mapi (fun i shell_a -> Array.map (fun shell_b -> 
+  Array.mapi (fun i shell_a ->
-        create shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis
+    Array.mapi (fun j shell_b -> 
+      create shell_a shell_b)
+    (Array.sub basis 0 (i+1)) 
+  ) basis
+
+
+let equivalent x y =
+  (Array.length x = Array.length y) &&
+  (Array.init (Array.length x) (fun k -> ShellPair.equivalent x.(k) y.(k))
+  |> Array.fold_left (fun accu x -> x && accu) true)
 
 
 (** A list of unique shell pairs *)
-let unique_of_shell_pairs sp =
+let unique sp =
+  let sp = 
     Array.to_list sp
     |> Array.concat
     |> Array.to_list
-  |> List.sort_uniq cmp
+  in
+  let rec aux accu = function
+  | [] -> accu
+  | x::rest -> 
+      try ignore @@ List.find (fun y -> equivalent x y) accu; aux accu rest
+      with Not_found -> aux (x::accu) rest
+  in
+  aux [] sp
 
 
-(** A map from a shell pair hash to the list of indices in the array of shell pairs. *)
+(** A map from a shell pair hash to the list of indices in the array
-let indices_of_shell_pairs sp =
+    of shell pairs.
+*)
+let indices sp =
   let map =
     Hashtbl.create 129
   in
@@ -126,7 +166,8 @@ let indices_of_shell_pairs sp =
       let key =
         hash shell_p
       in
-      Hashtbl.add map key (i,j); ) s ) sp;
+      Hashtbl.add map key (i,j); ) s
+  ) sp;
   map
       
 

Basis/ERI.ml

@@ -76,14 +76,12 @@ let to_file ~filename basis =
   let shell_pairs =
     ContractedShellPair.shell_pairs basis
   in
-
+  let indices_of_shell_pairs =
-(*
+    ContractedShellPair.indices shell_pairs
-  let unique_shell_pairs =
+  in
-    ContractedShellPair.unique_of_shell_pairs shell_pairs
+  let unique_shell_pairs =
-  and indices_of_shell_pairs =
+    ContractedShellPair.unique shell_pairs
-    ContractedShellPair.indices_of_shell_pairs shell_pairs
   in
-  *)
 
   (* Pre-compute diagonal integrals for Schwartz *)
   let t0 = Unix.gettimeofday () in
@@ -108,6 +106,7 @@ let to_file ~filename basis =
   done;
   Printf.printf "%d shell pairs computed in %f seconds\n" !icount (Unix.gettimeofday () -. t0);
 
+  (* 4D data initialization *)
   let eri_array = 
 
     let n = ref 0 in
@@ -124,15 +123,19 @@ let to_file ~filename basis =
 
   let inn = ref 0 and out = ref 0 in
 
+  (*
   for i=0 to (Array.length basis) - 1 do
     print_int (Contracted_shell.index basis.(i)) ; print_newline ();
     for j=0 to i do
+    *)
+  List.iter (fun shell_p ->  
+    let i,j = 
+      Hashtbl.find indices_of_shell_pairs (ContractedShellPair.hash shell_p)
+    in
+    assert (ContractedShellPair.equivalent (ContractedShellPair.create basis.(i) basis.(j)) shell_p);
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
       try
         if (schwartz_p_max < cutoff) then raise NullIntegral;
-        let
-          shell_p = shell_pairs.(i).(j)
-        in
 
         for k=0 to i do
           for l=0 to k do
@@ -162,18 +165,21 @@ let to_file ~filename basis =
                   contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q 
             in
 
+      Hashtbl.find_all indices_of_shell_pairs (ContractedShellPair.hash shell_p)
+      |> List.iter (fun (i,j) ->
+            Printf.printf "%d %d\n" i j;
             (* Write the data in the output file *)
             Array.iteri (fun i_c powers_i ->
-                let i_c = (Contracted_shell.index basis.(i)) + i_c + 1 in
+                let i_c = (Contracted_shell.index basis.(i)) + i_c  in
                 let xi = to_int_tuple powers_i in
                 Array.iteri (fun j_c powers_j ->
-                    let j_c = (Contracted_shell.index basis.(j)) + j_c + 1 in
+                    let j_c = (Contracted_shell.index basis.(j)) + j_c  in
                     let xj = to_int_tuple powers_j in
                     Array.iteri (fun k_c powers_k ->
-                        let k_c = (Contracted_shell.index basis.(k)) + k_c + 1 in 
+                        let k_c = (Contracted_shell.index basis.(k)) + k_c  in 
                         let xk = to_int_tuple powers_k in
                         Array.iteri (fun l_c powers_l ->
-                            let l_c = (Contracted_shell.index basis.(l)) + l_c + 1 in
+                            let l_c = (Contracted_shell.index basis.(l)) + l_c  in
                             let xl = to_int_tuple powers_l in
                             let key = 
                               if swap then
@@ -184,14 +190,14 @@ let to_file ~filename basis =
                             let value = 
                               Zmap.find cls key 
                             in
-                              eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} <- value;
+                              eri_array.{i_c,k_c,j_c,l_c} <- value;
-                              eri_array.{(j_c-1),(k_c-1),(i_c-1),(l_c-1)} <- value;
+                              eri_array.{j_c,k_c,i_c,l_c} <- value;
-                              eri_array.{(i_c-1),(l_c-1),(j_c-1),(k_c-1)} <- value;
+                              eri_array.{i_c,l_c,j_c,k_c} <- value;
-                              eri_array.{(j_c-1),(l_c-1),(i_c-1),(k_c-1)} <- value;
+                              eri_array.{j_c,l_c,i_c,k_c} <- value;
-                              eri_array.{(k_c-1),(i_c-1),(l_c-1),(j_c-1)} <- value;
+                              eri_array.{k_c,i_c,l_c,j_c} <- value;
-                              eri_array.{(k_c-1),(j_c-1),(l_c-1),(i_c-1)} <- value;
+                              eri_array.{k_c,j_c,l_c,i_c} <- value;
-                              eri_array.{(l_c-1),(i_c-1),(k_c-1),(j_c-1)} <- value;
+                              eri_array.{l_c,i_c,k_c,j_c} <- value;
-                              eri_array.{(l_c-1),(j_c-1),(k_c-1),(i_c-1)} <- value;
+                              eri_array.{l_c,j_c,k_c,i_c} <- value;
                             if (abs_float value > cutoff) then
                               (inn := !inn + 1;
                               )
@@ -201,12 +207,16 @@ let to_file ~filename basis =
                       ) (Contracted_shell.powers basis.(k))
                   ) (Contracted_shell.powers basis.(j))
               ) (Contracted_shell.powers basis.(i));
+            );
             with NullIntegral ->  ()
           done;
         done;
       with NullIntegral -> ()
+      ) unique_shell_pairs;
+      (*
     done;
   done;
+  *)
 
   Printf.printf "Computed %d non-zero ERIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0);
 

Basis/ShellPair.ml

@@ -22,7 +22,13 @@ type t = {
 
 (** Returns an integer characteristic of a primitive shell pair *)
 let hash a = 
-  Hashtbl.hash (a.expo, a.center_a, a.center_ab, a.coef) 
+  Hashtbl.hash a
+
+let equivalent a b =
+  a = b
+(*
+  Hashtbl.hash (a.expo, a.center_a, a.center_ab, a.coef, Contracted_shell.totAngMom a.shell_a, Contracted_shell.totAngMom a.shell_b) 
+*)
 
 
 (** Comparison function, used for sorting *)