diff --git a/Basis/Shell_pair.ml b/Basis/ContractedShellPair.ml similarity index 68% rename from Basis/Shell_pair.ml rename to Basis/ContractedShellPair.ml index 37d7ed2..14ae569 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/ContractedShellPair.ml @@ -1,27 +1,15 @@ open Util open Constants -type t = { - expo : float; (* alpha + beta *) - expo_inv : float; (* 1/(alpha + beta) *) - center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *) - center_a : Coordinate.t; (* P - A *) - center_ab: Coordinate.t; (* A - B *) - norm_sq : float; (* |A-B|^2 *) - norm_coef: float; (* norm_coef_a * norm_coef_b *) - coef : float; (* norm_coef * coef_a * coef_b * g, with - g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *) - norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *) - i : int; - j : int; - shell_a : Contracted_shell.t; - shell_b : Contracted_shell.t; - monocentric : bool -} +type t = ShellPair.t array exception Null_contribution -let create_array ?cutoff p_a p_b = +(** Creates an contracted shell pair : an array of pairs of primitive shells. + A contracted shell with N functions combined with a contracted + shell with M functions generates a NxM array of shell pairs. +*) +let create ?cutoff p_a p_b = let cutoff, log_cutoff = match cutoff with | None -> -1., max_float @@ -91,7 +79,7 @@ let create_array ?cutoff p_a p_b = let monocentric = Contracted_shell.center p_a = Contracted_shell.center p_b in - Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale ; monocentric } + Some ShellPair.{ i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale ; monocentric } with | Null_contribution -> None ) @@ -104,25 +92,43 @@ let create_array ?cutoff p_a p_b = |> Array.of_list +(** Returns an integer characteristic of a contracted shell pair *) let hash a = - Hashtbl.hash (a.expo, a.center_a, a.center_ab, a.coef) + 1 (*TODO*) + +(** Comparison function, used for sorting *) let cmp a b = hash a - hash b +(** The array of all shell pairs *) let shell_pairs basis = Array.mapi (fun i shell_a -> Array.map (fun shell_b -> - create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis + create shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis -let unique_shell_pairs basis = - let map = +(** 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 + + +(** A map from a shell pair hash to the list of indices in the array of shell pairs. *) +let indices_of_shell_pairs sp = + let map = Hashtbl.create 129 in - Array.iter (fun spa -> - Array.iter (fun value -> - let key = hash value in - Hashtbl.add map key value) spa) basis; + Array.iteri (fun i s -> + Array.iteri (fun j shell_p -> + let key = + hash shell_p + in + Hashtbl.add map key (i,j); ) s ) sp; map + + + diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 8468fc5..8255358 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -70,11 +70,20 @@ let to_file ~filename basis = let oc = open_out filename in + Printf.printf "%d shells\n" (Array.length basis); + (* Pre-compute all shell pairs *) let shell_pairs = - Shell_pair.shell_pairs basis + ContractedShellPair.shell_pairs basis in - Printf.printf "%d shells\n" (Array.length basis); + +(* + let unique_shell_pairs = + ContractedShellPair.unique_of_shell_pairs shell_pairs + and indices_of_shell_pairs = + ContractedShellPair.indices_of_shell_pairs shell_pairs + in + *) (* Pre-compute diagonal integrals for Schwartz *) let t0 = Unix.gettimeofday () in @@ -124,11 +133,6 @@ let to_file ~filename basis = let shell_p = shell_pairs.(i).(j) in - (* - Printf.printf "%d %d : " i j; - Array.iter (fun i -> Printf.printf "%d " (Shell_pair.hash i)) shell_p; - print_newline (); - *) for k=0 to i do for l=0 to k do diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 31acd68..151d093 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -5,7 +5,7 @@ open Constants let contracted_class shell_a shell_b : float Zmap.t = let shell_p = - Shell_pair.create_array shell_a shell_b + ContractedShellPair.create shell_a shell_b in (* Pre-computation of integral class indices *) @@ -23,30 +23,30 @@ let contracted_class shell_a shell_b : float Zmap.t = for ab=0 to (Array.length shell_p - 1) do let coef_prod = - shell_p.(ab).Shell_pair.coef + shell_p.(ab).ShellPair.coef in (** Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-4*.cutoff then begin let center_ab = - shell_p.(ab).Shell_pair.center_ab + shell_p.(ab).ShellPair.center_ab in let center_pa = - shell_p.(ab).Shell_pair.center_a + shell_p.(ab).ShellPair.center_a in let expo_inv = - shell_p.(ab).Shell_pair.expo_inv + shell_p.(ab).ShellPair.expo_inv in let norm_coef_scale = - shell_p.(ab).Shell_pair.norm_coef_scale + shell_p.(ab).ShellPair.norm_coef_scale in let i, j = - shell_p.(ab).Shell_pair.i, shell_p.(ab).Shell_pair.j + shell_p.(ab).ShellPair.i, shell_p.(ab).ShellPair.j in let expo_a = - Contracted_shell.expo shell_p.(ab).Shell_pair.shell_a i + Contracted_shell.expo shell_p.(ab).ShellPair.shell_a i and expo_b = - Contracted_shell.expo shell_p.(ab).Shell_pair.shell_b j + Contracted_shell.expo shell_p.(ab).ShellPair.shell_b j in Array.iteri (fun i key -> diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index 0e4bda4..d14c934 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -45,7 +45,7 @@ let to_file ~filename basis geometry = (* Pre-compute all shell pairs *) let shell_pairs = Array.mapi (fun i shell_a -> Array.map (fun shell_b -> - Shell_pair.create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis + ContractedShellPair.create shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis in Printf.printf "%d shells\n" (Array.length basis); diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index ea4a0b5..d764e42 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -123,8 +123,8 @@ let hvrr_one_e (** Computes all the one-electron integrals of the contracted shell pair *) let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = - let shell_a = shell_p.(0).Shell_pair.shell_a - and shell_b = shell_p.(0).Shell_pair.shell_b + let shell_a = shell_p.(0).ShellPair.shell_a + and shell_b = shell_p.(0).ShellPair.shell_b in let maxm = let open Angular_momentum in @@ -146,11 +146,11 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = for ab=0 to (Array.length shell_p - 1) do - let b = shell_p.(ab).Shell_pair.j in + let b = shell_p.(ab).ShellPair.j in try begin - let coef_prod = shell_p.(ab).Shell_pair.coef in - let norm_coef_scale_p = shell_p.(ab).Shell_pair.norm_coef_scale in + let coef_prod = shell_p.(ab).ShellPair.coef in + let norm_coef_scale_p = shell_p.(ab).ShellPair.norm_coef_scale in (** Screening on the product of coefficients *) if (abs_float coef_prod) < 1.e-4*.cutoff then @@ -158,14 +158,14 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let expo_pq_inv = - shell_p.(ab).Shell_pair.expo_inv + shell_p.(ab).ShellPair.expo_inv in let center_ab = - shell_p.(ab).Shell_pair.center_ab + shell_p.(ab).ShellPair.center_ab in let center_p = - shell_p.(ab).Shell_pair.center + shell_p.(ab).ShellPair.center in let center_pa = Coordinate.(center_p |- Contracted_shell.center shell_a) @@ -208,7 +208,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b) (maxm, zero_m_array) (Contracted_shell.expo shell_b b) - (shell_p.(ab).Shell_pair.expo_inv) + (shell_p.(ab).ShellPair.expo_inv) (center_ab, center_pa, center_pc) map in diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index f75d013..0e73d50 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -5,7 +5,7 @@ open Constants let contracted_class shell_a shell_b : float Zmap.t = let shell_p = - Shell_pair.create_array shell_a shell_b + ContractedShellPair.create shell_a shell_b in (* Pre-computation of integral class indices *) @@ -23,22 +23,22 @@ let contracted_class shell_a shell_b : float Zmap.t = for ab=0 to (Array.length shell_p - 1) do let coef_prod = - shell_p.(ab).Shell_pair.coef + shell_p.(ab).ShellPair.coef in (** Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-4*.cutoff then begin let center_ab = - shell_p.(ab).Shell_pair.center_ab + shell_p.(ab).ShellPair.center_ab in let center_pa = - shell_p.(ab).Shell_pair.center_a + shell_p.(ab).ShellPair.center_a in let expo_inv = - shell_p.(ab).Shell_pair.expo_inv + shell_p.(ab).ShellPair.expo_inv in let norm_coef_scale = - shell_p.(ab).Shell_pair.norm_coef_scale + shell_p.(ab).ShellPair.norm_coef_scale in Array.iteri (fun i key -> diff --git a/Basis/ShellPair.ml b/Basis/ShellPair.ml new file mode 100644 index 0000000..c81d6c8 --- /dev/null +++ b/Basis/ShellPair.ml @@ -0,0 +1,35 @@ +open Util +open Constants + +type t = { + expo : float; (* alpha + beta *) + expo_inv : float; (* 1/(alpha + beta) *) + center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *) + center_a : Coordinate.t; (* P - A *) + center_ab: Coordinate.t; (* A - B *) + norm_sq : float; (* |A-B|^2 *) + norm_coef: float; (* norm_coef_a * norm_coef_b *) + coef : float; (* norm_coef * coef_a * coef_b * g, with + g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *) + norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *) + i : int; + j : int; + shell_a : Contracted_shell.t; + shell_b : Contracted_shell.t; + monocentric : bool +} + + +(** Returns an integer characteristic of a primitive shell pair *) +let hash a = + Hashtbl.hash (a.expo, a.center_a, a.center_ab, a.coef) + + +(** Comparison function, used for sorting *) +let cmp a b = + hash a - hash b + + + + + diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index af1448d..9bc6986 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -299,10 +299,10 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - let shell_a = shell_p.(0).Shell_pair.shell_a - and shell_b = shell_p.(0).Shell_pair.shell_b - and shell_c = shell_q.(0).Shell_pair.shell_a - and shell_d = shell_q.(0).Shell_pair.shell_b + let shell_a = shell_p.(0).ShellPair.shell_a + and shell_b = shell_p.(0).ShellPair.shell_b + and shell_c = shell_q.(0).ShellPair.shell_a + and shell_d = shell_q.(0).ShellPair.shell_b in let maxm = let open Angular_momentum in @@ -325,12 +325,12 @@ 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 *) for ab=0 to (Array.length shell_p - 1) do - let cab = shell_p.(ab).Shell_pair.coef in - let b = shell_p.(ab).Shell_pair.j in - let norm_coef_scale_p = shell_p.(ab).Shell_pair.norm_coef_scale in + let cab = shell_p.(ab).ShellPair.coef in + let b = shell_p.(ab).ShellPair.j in + let norm_coef_scale_p = shell_p.(ab).ShellPair.norm_coef_scale in for cd=0 to (Array.length shell_q - 1) do let coef_prod = - cab *. shell_q.(cd).Shell_pair.coef + cab *. shell_q.(cd).ShellPair.coef in (** Screening on the product of coefficients *) try @@ -338,10 +338,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q raise NullQuartet; let expo_pq_inv = - shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv + shell_p.(ab).ShellPair.expo_inv +. shell_q.(cd).ShellPair.expo_inv in let center_pq = - Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center) + Coordinate.(shell_p.(ab).ShellPair.center |- shell_q.(cd).ShellPair.center) in let norm_pq_sq = Coordinate.dot center_pq center_pq @@ -359,10 +359,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral | _ -> - let d = shell_q.(cd).Shell_pair.j in + let d = shell_q.(cd).ShellPair.j in let map_1d = Zmap.create (4*maxm) in let map_2d = Zmap.create (Array.length class_indices) in - let norm_coef_scale_q = shell_q.(cd).Shell_pair.norm_coef_scale in + let norm_coef_scale_q = shell_q.(cd).ShellPair.norm_coef_scale in let norm_coef_scale = Array.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q @@ -372,8 +372,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in (* let monocentric = - shell_p.(ab).Shell_pair.monocentric && - shell_q.(cd).Shell_pair.monocentric + shell_p.(ab).ShellPair.monocentric && + shell_q.(cd).ShellPair.monocentric in *) (* Compute the integral class from the primitive shell quartet *) @@ -433,8 +433,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) (maxm, zero_m_array) (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d) - (shell_p.(ab).Shell_pair.expo_inv, shell_q.(cd).Shell_pair.expo_inv) - (shell_p.(ab).Shell_pair.center_ab, shell_q.(cd).Shell_pair.center_ab, center_pq) + (shell_p.(ab).ShellPair.expo_inv, shell_q.(cd).ShellPair.expo_inv) + (shell_p.(ab).ShellPair.center_ab, shell_q.(cd).ShellPair.center_ab, center_pq) map_1d map_2d in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral @@ -456,8 +456,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (** Computes all the two-electron integrals of the contracted shell quartet *) let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = - let shell_p = Shell_pair.create_array ~cutoff shell_a shell_b - and shell_q = Shell_pair.create_array ~cutoff shell_c shell_d + let shell_p = ContractedShellPair.create ~cutoff shell_a shell_b + and shell_q = ContractedShellPair.create ~cutoff shell_c shell_d in contracted_class_shell_pairs ~zero_m shell_p shell_q diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 812c4b6..247b2b0 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -282,10 +282,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - let shell_a = shell_p.(0).Shell_pair.shell_a - and shell_b = shell_p.(0).Shell_pair.shell_b - and shell_c = shell_q.(0).Shell_pair.shell_a - and shell_d = shell_q.(0).Shell_pair.shell_b + let shell_a = shell_p.(0).ShellPair.shell_a + and shell_b = shell_p.(0).ShellPair.shell_b + and shell_c = shell_q.(0).ShellPair.shell_a + and shell_d = shell_q.(0).ShellPair.shell_b in let maxm = let open Angular_momentum in @@ -316,7 +316,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (fun accu shell_ab -> accu +. Array.fold_left (fun accu shell_cd -> let coef_prod = - shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef + shell_ab.ShellPair.coef *. shell_cd.ShellPair.coef in (** Screening on the product of coefficients *) try @@ -324,10 +324,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q raise NullQuartet; let expo_pq_inv = - shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv + shell_ab.ShellPair.expo_inv +. shell_cd.ShellPair.expo_inv in let center_pq = - Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) + Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) in let norm_pq_sq = Coordinate.dot center_pq center_pq @@ -344,18 +344,18 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q | _ -> Array.iter (fun shell_ab -> - let norm_coef_scale_p = shell_ab.Shell_pair.norm_coef_scale in - let b = shell_ab.Shell_pair.j in + let norm_coef_scale_p = shell_ab.ShellPair.norm_coef_scale in + let b = shell_ab.ShellPair.j in let common = Array.map (fun shell_cd -> let coef_prod = - shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef + shell_ab.ShellPair.coef *. shell_cd.ShellPair.coef in let expo_pq_inv = - shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv + shell_ab.ShellPair.expo_inv +. shell_cd.ShellPair.expo_inv in let center_pq = - Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) + Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) in let norm_pq_sq = Coordinate.dot center_pq center_pq @@ -365,10 +365,10 @@ 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 - let d = shell_cd.Shell_pair.j in + let d = shell_cd.ShellPair.j in - (zero_m_array, shell_cd.Shell_pair.expo_inv, - Contracted_shell.expo shell_d d, shell_cd.Shell_pair.center_ab, + (zero_m_array, shell_cd.ShellPair.expo_inv, + Contracted_shell.expo shell_d d, shell_cd.ShellPair.center_ab, center_pq,coef_prod) ) shell_q |> Array.to_list @@ -407,7 +407,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let map_1d = Array.init maxm (fun _ -> Zmap.create (4*maxm)) in let map_2d = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in let norm = - let norm_coef_scale_q = shell_q.(0).Shell_pair.norm_coef_scale in + let norm_coef_scale_q = shell_q.(0).ShellPair.norm_coef_scale in Array.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q ) norm_coef_scale_p @@ -428,8 +428,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) (maxm, zero_m_array) (Contracted_shell.expo shell_b b, d) - (shell_ab.Shell_pair.expo_inv, expo_inv) - (shell_ab.Shell_pair.center_ab, center_cd, center_pq) + (shell_ab.ShellPair.expo_inv, expo_inv) + (shell_ab.ShellPair.center_ab, center_cd, center_pq) coef_prod map_1d map_2d in let x = Array.fold_left (+.) 0. integral in @@ -450,8 +450,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (** Computes all the two-electron integrals of the contracted shell quartet *) let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = - let shell_p = Shell_pair.create_array ~cutoff shell_a shell_b - and shell_q = Shell_pair.create_array ~cutoff shell_c shell_d + let shell_p = ContractedShellPair.create ~cutoff shell_a shell_b + and shell_q = ContractedShellPair.create ~cutoff shell_c shell_d in contracted_class_shell_pairs ~zero_m shell_p shell_q