From 365735d111210d4bff777ae33e6d11b37ba95d11 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 9 Feb 2018 19:41:22 +0100 Subject: [PATCH] Working on contraction --- Basis/ContractedShellPair.ml | 132 +++++----- Basis/ERI.ml | 13 +- Basis/KinInt.ml | 28 ++- Basis/OneElectronRR.ml | 21 +- Basis/Overlap.ml | 26 +- Basis/ShellPair.ml | 1 - Basis/TwoElectronRR.ml | 36 +-- Basis/TwoElectronRRVectorized.ml | 405 +++++++++++++++++-------------- 8 files changed, 366 insertions(+), 296 deletions(-) diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 20a1be1..632def2 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -6,7 +6,17 @@ 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 +type t = +{ + shell_a : Contracted_shell.t; + shell_b : Contracted_shell.t; + shell_pairs : ShellPair.t array; + coef : float array; + expo_inv : float array; + center_ab : Coordinate.t; (* A-B *) + norm_sq : float; (* |A-B|^2 *) + norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *) +} (** Creates an contracted shell pair : an array of pairs of primitive shells. @@ -38,62 +48,72 @@ let create ?cutoff p_a p_b = |> Array.to_list |> Array.concat in - Array.init (Contracted_shell.size p_a) (fun i -> - let p_a_expo_center = Coordinate.( - Contracted_shell.expo p_a i |. Contracted_shell.center p_a ) - in - let norm_coef_a = - Contracted_shell.norm_coef p_a i - in + let shell_pairs = + Array.init (Contracted_shell.size p_a) (fun i -> + let p_a_expo_center = Coordinate.( + Contracted_shell.expo p_a i |. Contracted_shell.center p_a ) + in + let norm_coef_a = + Contracted_shell.norm_coef p_a i + in - Array.init (Contracted_shell.size p_b) (fun j -> - try - let norm_coef_b = - Contracted_shell.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 = Coordinate.( - Contracted_shell.expo p_b j |. Contracted_shell.center p_b ) - in - let expo = Contracted_shell.(expo p_a i +. expo p_b j) in - let expo_inv = 1. /. expo in - let center = Coordinate.( - expo_inv |. (p_a_expo_center |+ p_b_expo_center ) ) - in - let argexpo = - Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv - in - if (argexpo > log_cutoff) then - raise Null_contribution; - let g = - (pi *. expo_inv)**(1.5) *. exp(-. argexpo) - in - let coef = - norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g - in - if (abs_float coef < cutoff) then - raise Null_contribution; - let center_a = - Coordinate.(center |- Contracted_shell.center p_a) - in - let monocentric = - Contracted_shell.center p_a = Contracted_shell.center p_b - in - 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 - ) - ) - |> Array.to_list - |> Array.concat - |> Array.to_list - |> List.filter (function Some _ -> true | None -> false) - |> List.map (function Some x -> x | None -> assert false) - |> Array.of_list + Array.init (Contracted_shell.size p_b) (fun j -> + try + let norm_coef_b = + Contracted_shell.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 = Coordinate.( + Contracted_shell.expo p_b j |. Contracted_shell.center p_b ) + in + let expo = Contracted_shell.(expo p_a i +. expo p_b j) in + let expo_inv = 1. /. expo in + let center = Coordinate.( + expo_inv |. (p_a_expo_center |+ p_b_expo_center ) ) + in + let argexpo = + Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv + in + if (argexpo > log_cutoff) then + raise Null_contribution; + let g = + (pi *. expo_inv)**(1.5) *. exp(-. argexpo) + in + let coef = + norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g + in + if (abs_float coef < cutoff) then + raise Null_contribution; + let center_a = + Coordinate.(center |- Contracted_shell.center p_a) + in + let monocentric = + Contracted_shell.center p_a = Contracted_shell.center p_b + in + 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 ; monocentric } + with + | Null_contribution -> None + ) + ) + |> Array.to_list + |> Array.concat + |> Array.to_list + |> List.filter (function Some _ -> true | None -> false) + |> List.map (function Some x -> x | None -> assert false) + |> Array.of_list + in + let coef = Array.map (fun x -> (fun y -> y.ShellPair.coef) x) shell_pairs + and expo_inv = Array.map (fun x -> (fun y -> y.ShellPair.expo_inv) x) shell_pairs + in + { + shell_a = p_a ; shell_b = p_b ; coef ; expo_inv ; + shell_pairs ; center_ab=shell_pairs.(0).center_ab; + norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq + } (** Returns an integer characteristic of a contracted shell pair *) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 7e13d3f..d01c185 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -134,6 +134,10 @@ let of_basis basis = shell_p = shell_pairs.(i).(j) in + let sp = + shell_p.ContractedShellPair.shell_pairs + in + for k=0 to i do for l=0 to k do let schwartz_q, schwartz_q_max = schwartz.(k).(l) in @@ -143,20 +147,23 @@ let of_basis basis = let shell_q = shell_pairs.(k).(l) in + let sq = + shell_q.ContractedShellPair.shell_pairs + in let swap = - Array.length shell_q < Array.length shell_p + Array.length sp > Array.length sq in (* Compute all the integrals of the class *) let cls = if swap then - if Array.length shell_p < 2 then + if Array.length sp < 2 then contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p else contracted_class_shell_pairs_vec ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p else - if Array.length shell_q < 2 then + if Array.length sq < 2 then contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q else contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index d192137..57106bc 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -24,33 +24,35 @@ let contracted_class shell_a shell_b : float Zmap.t = (* Compute all integrals in the shell for each pair of significant shell pairs *) - for ab=0 to (Array.length shell_p - 1) + let sp = shell_p.ContractedShellPair.shell_pairs in + let center_ab = + shell_p.ContractedShellPair.center_ab + in + let norm_coef_scale = + shell_p.ContractedShellPair.norm_coef_scale + in + + for ab=0 to (Array.length sp - 1) do let coef_prod = - shell_p.(ab).ShellPair.coef + shell_p.ContractedShellPair.coef.(ab) in (** Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-4*.cutoff then begin - let center_ab = - shell_p.(ab).ShellPair.center_ab - in let center_pa = - shell_p.(ab).ShellPair.center_a + sp.(ab).ShellPair.center_a in let expo_inv = - shell_p.(ab).ShellPair.expo_inv - in - let norm_coef_scale = - shell_p.(ab).ShellPair.norm_coef_scale + shell_p.ContractedShellPair.expo_inv.(ab) in let i, j = - shell_p.(ab).ShellPair.i, shell_p.(ab).ShellPair.j + sp.(ab).ShellPair.i, sp.(ab).ShellPair.j in let expo_a = - Contracted_shell.expo shell_p.(ab).ShellPair.shell_a i + Contracted_shell.expo sp.(ab).ShellPair.shell_a i and expo_b = - Contracted_shell.expo shell_p.(ab).ShellPair.shell_b j + Contracted_shell.expo sp.(ab).ShellPair.shell_b j in Array.iteri (fun i key -> diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index d764e42..f517ba2 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).ShellPair.shell_a - and shell_b = shell_p.(0).ShellPair.shell_b + let shell_a = shell_p.ContractedShellPair.shell_a + and shell_b = shell_p.ContractedShellPair.shell_b in let maxm = let open Angular_momentum in @@ -144,13 +144,14 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = (* Compute all integrals in the shell for each pair of significant shell pairs *) - for ab=0 to (Array.length shell_p - 1) + let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale + in + for ab=0 to (Array.length shell_p.ContractedShellPair.shell_pairs - 1) do - let b = shell_p.(ab).ShellPair.j in + let b = shell_p.ContractedShellPair.shell_pairs.(ab).ShellPair.j in try begin - let coef_prod = shell_p.(ab).ShellPair.coef in - let norm_coef_scale_p = shell_p.(ab).ShellPair.norm_coef_scale in + let coef_prod = shell_p.ContractedShellPair.coef.(ab) in (** Screening on the product of coefficients *) if (abs_float coef_prod) < 1.e-4*.cutoff then @@ -158,14 +159,14 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let expo_pq_inv = - shell_p.(ab).ShellPair.expo_inv + shell_p.ContractedShellPair.expo_inv.(ab) in let center_ab = - shell_p.(ab).ShellPair.center_ab + shell_p.ContractedShellPair.center_ab in let center_p = - shell_p.(ab).ShellPair.center + shell_p.ContractedShellPair.shell_pairs.(ab).ShellPair.center in let center_pa = Coordinate.(center_p |- Contracted_shell.center shell_a) @@ -208,7 +209,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).ShellPair.expo_inv) + (shell_p.ContractedShellPair.expo_inv.(ab)) (center_ab, center_pa, center_pc) map in diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 740d40c..d65e22a 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -22,27 +22,31 @@ let contracted_class shell_a shell_b : float Zmap.t = Array.make (Array.length class_indices) 0. in + let sp = + shell_p.ContractedShellPair.shell_pairs + in + let center_ab = + shell_p.ContractedShellPair.center_ab + in + let norm_coef_scale = + shell_p.ContractedShellPair.norm_coef_scale + in + (* Compute all integrals in the shell for each pair of significant shell pairs *) - for ab=0 to (Array.length shell_p - 1) + for ab=0 to (Array.length sp - 1) do let coef_prod = - shell_p.(ab).ShellPair.coef + shell_p.ContractedShellPair.coef.(ab) in (** Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-4*.cutoff then begin - let center_ab = - shell_p.(ab).ShellPair.center_ab + let expo_inv = + shell_p.ContractedShellPair.expo_inv.(ab) in let center_pa = - shell_p.(ab).ShellPair.center_a - in - let expo_inv = - shell_p.(ab).ShellPair.expo_inv - in - let norm_coef_scale = - shell_p.(ab).ShellPair.norm_coef_scale + sp.(ab).ShellPair.center_a in Array.iteri (fun i key -> diff --git a/Basis/ShellPair.ml b/Basis/ShellPair.ml index 60cbe10..cb9736e 100644 --- a/Basis/ShellPair.ml +++ b/Basis/ShellPair.ml @@ -11,7 +11,6 @@ type t = { 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; diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 2a24c48..64fbdc7 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -310,10 +310,12 @@ let rec 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).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 + let shell_a = shell_p.ContractedShellPair.shell_a + and shell_b = shell_p.ContractedShellPair.shell_b + and shell_c = shell_q.ContractedShellPair.shell_a + and shell_d = shell_q.ContractedShellPair.shell_b + and sp = shell_p.ContractedShellPair.shell_pairs + and sq = shell_q.ContractedShellPair.shell_pairs in let maxm = let open Angular_momentum in @@ -335,13 +337,14 @@ 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).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 norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in + let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in + for ab=0 to (Array.length sp - 1) do + let cab = shell_p.ContractedShellPair.coef.(ab) in + let b = sp.(ab).ShellPair.j in + for cd=0 to (Array.length shell_q.ContractedShellPair.shell_pairs - 1) do let coef_prod = - cab *. shell_q.(cd).ShellPair.coef + cab *. shell_q.ContractedShellPair.coef.(cd) in (** Screening on the product of coefficients *) try @@ -349,10 +352,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q raise NullQuartet; let expo_pq_inv = - shell_p.(ab).ShellPair.expo_inv +. shell_q.(cd).ShellPair.expo_inv + shell_p.ContractedShellPair.expo_inv.(ab) +. + shell_q.ContractedShellPair.expo_inv.(cd) in let center_pq = - Coordinate.(shell_p.(ab).ShellPair.center |- shell_q.(cd).ShellPair.center) + Coordinate.(sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center) in let norm_pq_sq = Coordinate.dot center_pq center_pq @@ -370,12 +374,11 @@ 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).ShellPair.j in + let d = shell_q.ContractedShellPair.shell_pairs.(cd).ShellPair.j in let map_1d = Zmap.create (4*maxm) in let map_2d = Zmap.create (Array.length class_indices) 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).ShellPair.norm_coef_scale in let norm_coef_scale = Array.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q @@ -446,8 +449,9 @@ 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).ShellPair.expo_inv, shell_q.(cd).ShellPair.expo_inv) - (shell_p.(ab).ShellPair.center_ab, shell_q.(cd).ShellPair.center_ab, center_pq) + (shell_p.ContractedShellPair.expo_inv.(ab), + shell_q.ContractedShellPair.expo_inv.(cd) ) + (sp.(ab).ShellPair.center_ab, sq.(cd).ShellPair.center_ab, center_pq) map_1d map_2d map_1d' map_2d' in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 79a9f29..33021b3 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -1,4 +1,6 @@ open Util +open Lacaml.D +open Bigarray let cutoff = Constants.cutoff let cutoff2 = cutoff *. cutoff @@ -8,9 +10,13 @@ exception Found let at_least_one_valid arr = try - Array.fold_left (fun _ x -> if (abs_float x > cutoff) then raise Found else false ) false arr + Vec.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false with Found -> true +(*TODO : REMOVE *) +let sum integral = + Array.fold_left (+.) 0. integral + (** Horizontal and Vertical Recurrence Relations (HVRR) *) let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (totAngMom_a, totAngMom_b, totAngMom_c, totAngMom_d) @@ -18,12 +24,14 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) - coef_prod map_1d map_2d + coef_prod map_1d map_2d = - let ncoef = (Array.length coef_prod) in + let nq = Mat.dim1 coef_prod in + let np = Mat.dim2 coef_prod in + let empty = - Array.make ncoef 0. + Array.make nq 0. in let totAngMom_a = Angular_momentum.to_int totAngMom_a @@ -33,7 +41,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in (** Vertical recurrence relations *) - let rec vrr0_v m angMom_a = function + let rec vrr0_v l m angMom_a = function | 1 -> let xyz = match angMom_a with @@ -41,16 +49,16 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | (_,1,_) -> 1 | _ -> 2 in - let f = expo_b *. (Coordinate.coord center_ab xyz) in - Array.init ncoef (fun k -> coef_prod.(k) *. expo_inv_p *. - ( (Coordinate.coord center_pq.(k) xyz) *. zero_m_array.(m+1).(k) - -. f *. zero_m_array.(m).(k) ) ) - | 0 -> Array.map2 ( *. ) zero_m_array.(m) coef_prod + let f = expo_b.{l} *. (Coordinate.coord center_ab xyz) in + Array.init nq (fun k -> coef_prod.{k+1,l} *. expo_inv_p.{l} *. + (center_pq.{xyz+1,k+1,l} *. zero_m_array.(m+1).{k+1,l} + -. f *. zero_m_array.(m).{k+1,l} ) ) + | 0 -> Array.init nq (fun k -> zero_m_array.(m).{k+1,l} *. coef_prod.{k+1,l}) | totAngMom_a -> let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in - try Zmap.find map_1d.(m) key with + try Zmap.find map_1d.(m).(l-1) key with | Not_found -> let result = let am, amm, amxyz, xyz = @@ -62,43 +70,43 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) if amxyz < 0 then empty else let v1 = let f = - -. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz) + -. expo_b.{l} *. expo_inv_p.{l} *. (Coordinate.coord center_ab xyz) in if (abs_float f < cutoff) then empty else - Array.map (fun v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) ) + Array.map (fun v1k -> f *. v1k) (vrr0_v l m am (totAngMom_a-1) ) in let p1 = - Array.mapi (fun k v2k -> v1.(k) +. expo_inv_p *. (Coordinate.coord center_pq.(k) xyz) *. v2k) (vrr0_v (m+1) am (totAngMom_a-1)) + Array.mapi (fun k v2k -> v1.(k) +. expo_inv_p.{l} *. center_pq.{xyz+1,k+1,l} *. v2k) (vrr0_v l (m+1) am (totAngMom_a-1)) in if amxyz < 1 then p1 else - let f = (float_of_int amxyz) *. expo_inv_p *. 0.5 + let f = (float_of_int amxyz) *. expo_inv_p.{l} *. 0.5 in if (abs_float f < cutoff) then empty else - let v1 = vrr0_v m amm (totAngMom_a-2) + let v1 = vrr0_v l m amm (totAngMom_a-2) in let v2 = - if (abs_float (f *. expo_inv_p)) < cutoff then empty else - vrr0_v (m+1) amm (totAngMom_a-2) + if (abs_float (f *. expo_inv_p.{l})) < cutoff then empty else + vrr0_v l (m+1) amm (totAngMom_a-2) in - Array.init ncoef (fun k -> p1.(k) +. - f *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) - in Zmap.add map_1d.(m) key result; + Array.init nq (fun k -> p1.(k) +. + f *. (v1.(k) +. v2.(k) *. expo_inv_p.{l} ) ) + in Zmap.add map_1d.(m).(l-1) key result; result - and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c = + and vrr_v l m angMom_a angMom_c totAngMom_a totAngMom_c = match (totAngMom_a, totAngMom_c) with | (i,0) -> if (i>0) then - vrr0_v m angMom_a totAngMom_a + vrr0_v l m angMom_a totAngMom_a else - Array.map2 ( *. ) zero_m_array.(m) coef_prod + Array.init nq (fun k -> zero_m_array.(m).{k+1,l} *. coef_prod.{k+1,l}) | (_,_) -> let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c)) in - try Zmap.find map_2d.(m) key with + try Zmap.find map_2d.(m).(l-1) key with | Not_found -> let result = begin @@ -126,24 +134,26 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) *) let f1 = let f = (Coordinate.coord center_cd xyz) in - Array.init ncoef (fun k -> - expo_d.(k) *. expo_inv_q.(k) *. f) + Array.init nq (fun k -> + expo_d.{k+1} *. expo_inv_q.{k+1} *. f) + |> Vec.of_array in let f2 = - Array.init ncoef (fun k -> - expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz) ) + Array.init nq (fun k -> + expo_inv_q.{k+1} *. center_pq.{xyz+1,k+1,l} ) + |> Vec.of_array in let v1 = if (at_least_one_valid f1) then - vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) + vrr_v l m angMom_a cm totAngMom_a (totAngMom_c-1) else empty and v2 = if (at_least_one_valid f2) then - vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) + vrr_v l (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) else empty in let p1 = - Array.init ncoef (fun k -> -. v1.(k) *. f1.(k) -. v2.(k) *. f2.(k)) + Array.init nq (fun k -> -. v1.(k) *. f1.{k+1} -. v2.(k) *. f2.{k+1}) in let p2 = if cxyz < 2 then p1 else @@ -151,54 +161,56 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (float_of_int (cxyz-1)) *. 0.5 in let f1 = - Array.map (fun e -> fcm *. e) expo_inv_q + Vec.map (fun e -> fcm *. e) expo_inv_q in let f2 = - Array.map2 ( *. ) f1 expo_inv_q + Vec.mul f1 expo_inv_q in let v1 = if (at_least_one_valid f1) then - vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) + vrr_v l m angMom_a cmm totAngMom_a (totAngMom_c-2) else empty in let v2 = if (at_least_one_valid f2) then - vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) + vrr_v l (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) else empty in - Array.init ncoef (fun k -> p1.(k) +. f1.(k) *. v1.(k) +. f2.(k) *. v2.(k)) + Array.init nq (fun k -> p1.(k) +. f1.{k+1} *. v1.(k) +. f2.{k+1} *. v2.(k)) in if (axyz < 1) || (cxyz < 1) then p2 else let fa = - (float_of_int axyz) *. expo_inv_p *. 0.5 + (float_of_int axyz) *. expo_inv_p.{l} *. 0.5 in let f1 = - Array.map (fun e -> fa *. e ) expo_inv_q + Vec.map (fun e -> fa *. e ) expo_inv_q in if (at_least_one_valid f1) then let v = - vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) + vrr_v l (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) in - Array.init ncoef (fun k -> p2.(k) -. f1.(k) *. v.(k)) + Array.init nq (fun k -> p2.(k) -. f1.{k+1} *. v.(k)) else p2 end - in Zmap.add map_2d.(m) key result; + in Zmap.add map_2d.(m).(l-1) key result; result (** Horizontal recurrence relations *) - and hrr0_v angMom_a angMom_b angMom_c + and hrr0_v l angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c = match totAngMom_b with | 0 -> begin match (totAngMom_a, totAngMom_c) with - | (0,0) -> Array.map2 ( *. ) zero_m_array.(0) coef_prod - | (_,0) -> vrr0_v 0 angMom_a totAngMom_a - | (_,_) -> vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c + | (0,0) -> + Array.init nq (fun k -> zero_m_array.(0).{k+1,l} *. coef_prod.{k+1,l}) + |> sum + | (_,0) -> vrr0_v l 0 angMom_a totAngMom_a |> sum + | (_,_) -> vrr_v l 0 angMom_a angMom_c totAngMom_a totAngMom_c |> sum end | 1 -> let (aax, aay, aaz) = angMom_a in @@ -210,13 +222,13 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in let f = Coordinate.coord center_ab xyz in let v1 = - vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c + vrr_v l 0 ap angMom_c (totAngMom_a+1) totAngMom_c in - if (abs_float f < cutoff) then v1 else + if (abs_float f < cutoff) then sum v1 else let v2 = - vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c + vrr_v l 0 angMom_a angMom_c totAngMom_a totAngMom_c in - Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2 + Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2 |> sum | _ -> let (aax, aay, aaz) = angMom_a and (abx, aby, abz) = angMom_b in @@ -226,7 +238,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | (0,_,_) -> aby, 1 | _ -> abx, 0 in - if (bxyz < 1) then empty else + if (bxyz < 1) then 0. else let ap, bm = match xyz with | 0 -> (aax+1,aay,aaz),(abx-1,aby,abz) @@ -235,22 +247,24 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in let h1 = - hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c + hrr0_v l ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in let f = (Coordinate.coord center_ab xyz) in if (abs_float f < cutoff) then h1 else let h2 = - hrr0_v angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c - in Array.map2 (fun h1 h2 -> h1 +. h2 *. f) h1 h2 + hrr0_v l angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c + in + h1 +. h2 *. f - and hrr_v angMom_a angMom_b angMom_c angMom_d + and hrr_v l angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d = match (totAngMom_b, totAngMom_d) with | (_,0) -> if (totAngMom_b = 0) then - vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c + vrr_v l 0 angMom_a angMom_c totAngMom_a totAngMom_c + |> sum else - hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c + hrr0_v l angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> let (acx, acy, acz) = angMom_c and (adx, ady, adz) = angMom_d in @@ -261,19 +275,21 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | _ -> (acx, acy, acz+1), (adx, ady, adz-1), 2 in let h1 = - hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) + hrr_v l angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) and h2 = - hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) + hrr_v l angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in let f = (Coordinate.coord center_cd xyz) in - Array.init ncoef (fun k -> h1.(k) +. h2.(k) *. f) + h1 +. f *. h2 in - hrr_v - (angMom_a.(0),angMom_a.(1),angMom_a.(2)) - (angMom_b.(0),angMom_b.(1),angMom_b.(2)) - (angMom_c.(0),angMom_c.(1),angMom_c.(2)) - (angMom_d.(0),angMom_d.(1),angMom_d.(2)) - totAngMom_a totAngMom_b totAngMom_c totAngMom_d + Array.init np (fun ab -> + hrr_v (ab+1) + (angMom_a.(0),angMom_a.(1),angMom_a.(2)) + (angMom_b.(0),angMom_b.(1),angMom_b.(2)) + (angMom_c.(0),angMom_c.(1),angMom_c.(2)) + (angMom_d.(0),angMom_d.(1),angMom_d.(2)) + totAngMom_a totAngMom_b totAngMom_c totAngMom_d + ) |> sum @@ -283,10 +299,12 @@ 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).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 + let shell_a = shell_p.ContractedShellPair.shell_a + and shell_b = shell_p.ContractedShellPair.shell_b + and shell_c = shell_q.ContractedShellPair.shell_a + and shell_d = shell_q.ContractedShellPair.shell_b + and sp = shell_p.ContractedShellPair.shell_pairs + and sq = shell_q.ContractedShellPair.shell_pairs in let maxm = let open Angular_momentum in @@ -308,134 +326,149 @@ 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 *) + let expo_inv_p = + Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp + |> Vec.of_array + and expo_inv_q = + Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq + |> Vec.of_array + in + let np, nq = + Vec.dim expo_inv_p, Vec.dim expo_inv_q + in + let coef = + let result = Mat.make0 nq np in + Lacaml.D.ger + (Vec.of_array shell_q.ContractedShellPair.coef) + (Vec.of_array shell_p.ContractedShellPair.coef) + result; + result + in + begin match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d) with | Angular_momentum.(S,S,S,S) -> - contracted_class.(0) <- - Array.fold_left - (fun accu shell_ab -> accu +. - Array.fold_left (fun accu shell_cd -> - let coef_prod = - shell_ab.ShellPair.coef *. shell_cd.ShellPair.coef - in - (** Screening on the product of coefficients *) - try - if (abs_float coef_prod) < 1.e-3*.cutoff then - raise NullQuartet; + contracted_class.(0) <- + let zm_array = Mat.init_rows np nq (fun i j -> + (** Screening on the product of coefficients *) + try + if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then + raise NullQuartet; - let expo_pq_inv = - shell_ab.ShellPair.expo_inv +. shell_cd.ShellPair.expo_inv - in - let center_pq = - Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in - - let zero_m_array = - zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq - in - - accu +. coef_prod *. zero_m_array.(0) - with NullQuartet -> accu - ) 0. shell_q - ) 0. shell_p + let expo_pq_inv = + expo_inv_p.{i} +. expo_inv_q.{j} + in + let center_pq = + Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center) + in + let norm_pq_sq = + Coordinate.dot center_pq center_pq + in + let zero_m_array = + zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq + in + zero_m_array.(0) + with NullQuartet -> 0. + ) in + Mat.gemm_trace zm_array coef | _ -> - Array.iter (fun shell_ab -> - 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.ShellPair.coef *. shell_cd.ShellPair.coef - in - let expo_pq_inv = - shell_ab.ShellPair.expo_inv +. shell_cd.ShellPair.expo_inv - in - let center_pq = - Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in + let expo_b = + Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp + |> Vec.of_array + and expo_d = + Array.map (fun shell_cd -> Contracted_shell.expo shell_d shell_cd.ShellPair.j) sq + |> Vec.of_array + in + let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in + + + let center_pq = + let result = + Array3.create Float64 fortran_layout 3 nq np + in + Array.iteri (fun ab shell_ab -> + Array.iteri (fun cd shell_cd -> + let cpq = + Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) + in + result.{1,cd+1,ab+1} <- Coordinate.x cpq; + result.{2,cd+1,ab+1} <- Coordinate.y cpq; + result.{3,cd+1,ab+1} <- Coordinate.z cpq; + ) sq + ) sp; + result + in + let zero_m_array = + let result = + Array.init (maxm+1) (fun _ -> Mat.make0 nq np) + in + Array.iteri (fun ab shell_ab -> + let zero_m_array_tmp = + Array.mapi (fun cd shell_cd -> + let expo_pq_inv = + expo_inv_p.{ab+1} +. expo_inv_q.{cd+1} + in + let norm_pq_sq = + center_pq.{1,cd+1,ab+1} *. center_pq.{1,cd+1,ab+1} +. + center_pq.{2,cd+1,ab+1} *. center_pq.{2,cd+1,ab+1} +. + center_pq.{3,cd+1,ab+1} *. center_pq.{3,cd+1,ab+1} + in - let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - in - - let d = shell_cd.ShellPair.j in - - (zero_m_array, shell_cd.ShellPair.expo_inv, - Contracted_shell.expo shell_d d, - center_pq,coef_prod) - ) shell_q - |> Array.to_list - |> List.filter (fun (zero_m_array, expo_inv, d, - center_pq,coef_prod) -> abs_float coef_prod >= 1.e-4 *. cutoff) - |> Array.of_list - in - let zero_m_array = Array.map (fun (zero_m_array, expo_inv, d, - center_pq,coef_prod) -> zero_m_array) common - and expo_inv = Array.map (fun (zero_m_array, expo_inv, d, - center_pq,coef_prod) -> expo_inv ) common - and d = Array.map (fun (zero_m_array, expo_inv, d, - center_pq,coef_prod) -> d) common - and center_pq = Array.map (fun (zero_m_array, expo_inv, d, - center_pq,coef_prod) -> center_pq) common - and coef_prod = Array.map (fun (zero_m_array, expo_inv, d, - center_pq,coef_prod) -> coef_prod) common - in - (* Transpose zero_m_array - *) - let zero_m_array = - let result = Array.init (maxm+1) (fun _ -> - Array.make (Array.length coef_prod) 0.) + ) sq + (* + |> Array.to_list + |> List.filter (fun (zero_m_array, d, + center_pq,coef_prod) -> abs_float coef_prod >= 1.e-4 *. cutoff) + |> Array.of_list + *) in - for m=0 to maxm do - for k=0 to (Array.length coef_prod-1) do - result.(m).(k) <- zero_m_array.(k).(m) - done; - done; - result - in + (* Transpose result *) + for m=0 to maxm do + for cd=1 to nq do + result.(m).{cd,ab+1} <- zero_m_array_tmp.(cd-1).(m) + done + done + ) sp; + result + in + let norm = + let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in + Array.map (fun v1 -> + Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q + ) norm_coef_scale_p + |> Array.to_list + |> Array.concat + in + + let map_1d = Array.init maxm (fun _ -> Array.init np (fun _ -> Zmap.create (4*maxm))) + and map_2d = Array.init maxm (fun _ -> Array.init np (fun _ -> Zmap.create (Array.length class_indices))) + in (* Compute the integral class from the primitive shell quartet *) - 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).ShellPair.norm_coef_scale in - Array.map (fun v1 -> - Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q - ) norm_coef_scale_p - |> Array.to_list - |> Array.concat - in - Array.iteri (fun i key -> - let a = Zkey.to_int_array Zkey.Kind_12 key in - let (angMomA,angMomB,angMomC,angMomD) = - ( [| a.(0) ; a.(1) ; a.(2) |], - [| a.(3) ; a.(4) ; a.(5) |], - [| a.(6) ; a.(7) ; a.(8) |], - [| a.(9) ; a.(10) ; a.(11) |] ) - in - let integral = - hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD) - (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, - Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) - (maxm, zero_m_array) - (Contracted_shell.expo shell_b b, d) - (shell_ab.ShellPair.expo_inv, expo_inv) - (shell_p.(0).ShellPair.center_ab, - shell_q.(0).ShellPair.center_ab, center_pq) - coef_prod map_1d map_2d - in - let x = Array.fold_left (+.) 0. integral in - contracted_class.(i) <- contracted_class.(i) +. x *. norm.(i) - ) class_indices - ) shell_p + Array.iteri (fun i key -> + let a = Zkey.to_int_array Zkey.Kind_12 key in + let (angMomA,angMomB,angMomC,angMomD) = + ( [| a.(0) ; a.(1) ; a.(2) |], + [| a.(3) ; a.(4) ; a.(5) |], + [| a.(6) ; a.(7) ; a.(8) |], + [| a.(9) ; a.(10) ; a.(11) |] ) + in + let integral = + hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD) + (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, + Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) + (maxm, zero_m_array) + (expo_b, expo_d) + (expo_inv_p, expo_inv_q) + (shell_p.ContractedShellPair.center_ab, + shell_q.ContractedShellPair.center_ab, center_pq) + coef map_1d map_2d + in + contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i) + ) class_indices end;