diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 14ae569..20a1be1 100644 --- a/Basis/ContractedShellPair.ml +++ b/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 -> - create shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis + Array.mapi (fun i shell_a -> + 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 = - Array.to_list sp - |> Array.concat - |> Array.to_list - |> List.sort_uniq cmp +let unique sp = + let sp = + Array.to_list sp + |> Array.concat + |> Array.to_list + 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. *) -let indices_of_shell_pairs sp = +(** A map from a shell pair hash to the list of indices in the array + 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 diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 8255358..a3623cd 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -76,14 +76,12 @@ let to_file ~filename basis = let shell_pairs = ContractedShellPair.shell_pairs basis in - -(* - let unique_shell_pairs = - ContractedShellPair.unique_of_shell_pairs shell_pairs - and indices_of_shell_pairs = - ContractedShellPair.indices_of_shell_pairs shell_pairs + let indices_of_shell_pairs = + ContractedShellPair.indices shell_pairs + in + let unique_shell_pairs = + ContractedShellPair.unique 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.{(j_c-1),(k_c-1),(i_c-1),(l_c-1)} <- value; - eri_array.{(i_c-1),(l_c-1),(j_c-1),(k_c-1)} <- value; - eri_array.{(j_c-1),(l_c-1),(i_c-1),(k_c-1)} <- value; - eri_array.{(k_c-1),(i_c-1),(l_c-1),(j_c-1)} <- value; - eri_array.{(k_c-1),(j_c-1),(l_c-1),(i_c-1)} <- value; - eri_array.{(l_c-1),(i_c-1),(k_c-1),(j_c-1)} <- value; - eri_array.{(l_c-1),(j_c-1),(k_c-1),(i_c-1)} <- value; + eri_array.{i_c,k_c,j_c,l_c} <- value; + eri_array.{j_c,k_c,i_c,l_c} <- value; + eri_array.{i_c,l_c,j_c,k_c} <- value; + eri_array.{j_c,l_c,i_c,k_c} <- value; + eri_array.{k_c,i_c,l_c,j_c} <- value; + eri_array.{k_c,j_c,l_c,i_c} <- value; + eri_array.{l_c,i_c,k_c,j_c} <- 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 -> () + ); + 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); diff --git a/Basis/ShellPair.ml b/Basis/ShellPair.ml index c81d6c8..60cbe10 100644 --- a/Basis/ShellPair.ml +++ b/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 *)