From 03c5d1d64ae607a3b6ac8a29059780f7483e8b57 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 15 Mar 2018 16:03:43 +0100 Subject: [PATCH] Introduced PrimitiveShellPair --- Basis/ContractedShellPair.ml | 134 +++++++------------------------ Basis/ContractedShellPair.mli | 6 +- Basis/KinInt.ml | 12 ++- Basis/OneElectronRR.ml | 19 +++-- Basis/Overlap.ml | 4 +- Basis/PrimitiveShellPair.ml | 5 ++ Basis/PrimitiveShellPair.mli | 6 ++ Basis/ShellPair.ml | 40 --------- Basis/TwoElectronRR.ml | 54 ++++++------- Basis/TwoElectronRRVectorized.ml | 23 +++--- 10 files changed, 100 insertions(+), 203 deletions(-) delete mode 100644 Basis/ShellPair.ml diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index be17671..68f4fca 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -7,7 +7,7 @@ type t = { shell_a : ContractedShell.t; shell_b : ContractedShell.t; - shell_pairs : ShellPair.t array; + shell_pairs : PrimitiveShellPair.t array; coef : float array; expo_inv : float array; center_ab : Coordinate.t; (* A-B *) @@ -22,7 +22,6 @@ module Co = Coordinate module Cs = ContractedShell module Ps = PrimitiveShell module Psp = PrimitiveShellPair -module Sp = ShellPair (** Creates an contracted shell pair : an array of pairs of primitive shells. A contracted shell with N functions combined with a contracted @@ -37,115 +36,36 @@ Format.printf "@[<2>shell_b :@ %a@]@;" Cs.pp s_b; let make = Psp.create_make_of (Cs.prim s_a).(0) (Cs.prim s_b).(0) in -let center_ab = Co.( Cs.center s_a |- Cs.center s_b ) in - -(* -Format.printf "@[center_ab :@ %a@]@;" Coordinate.pp_angstrom center_ab; -Format.printf "@[a_minus_b :@ %a@]@." Coordinate.pp_angstrom (Psp.a_minus_b ( - match make 0 (Cs.prim s_a).(0) 0 (Cs.prim s_b).(0) 0. - with Some x -> x | _ -> assert false)); - *) - -let norm_sq = Co.dot center_ab center_ab in - -let norm_coef_scale_a = - Cs.norm_coef_scale s_a -and norm_coef_scale_b = - Cs.norm_coef_scale s_b -in -let norm_coef_scale = - Array.map (fun v1 -> - Array.map (fun v2 -> v1 *. v2) norm_coef_scale_b - ) norm_coef_scale_a - |> Array.to_list - |> Array.concat -in - -assert (norm_coef_scale = Psp.norm_coef_scale ( - match make 0 (Cs.prim s_a).(0) 0 (Cs.prim s_b).(0) 0. - with Some x -> x | _ -> assert false)); - let shell_pairs = - Array.init (Cs.size s_a) (fun i -> -let p_a = (Cs.prim s_a).(i) in -let p_a_expo_center = Co.( (Cs.expo s_a).(i) |. Cs.center s_a ) in -let norm_coef_a = (Cs.norm_coef s_a).(i) in -assert (norm_coef_a = Ps.norm_coef p_a); - - let make = make 0 p_a in - - Array.init (Cs.size s_b) (fun j -> -let p_b = (Cs.prim s_b).(j) in - - try - let sp = make 0 p_b cutoff in - -let sp_ = match sp with Some x -> x | None -> raise Null_contribution in - - let norm_coef_b = (Cs.norm_coef s_b).(j) in -assert (norm_coef_b = Ps.norm_coef p_b); - - let norm_coef = norm_coef_a *. norm_coef_b in - - if norm_coef < cutoff then - raise Null_contribution; - let p_b_expo_center = Co.( (Cs.expo s_b).(j) |. Cs.center s_b ) in - let expo = (Cs.expo s_a).(i) +. (Cs.expo s_b).(j) in - let expo_inv = 1. /. expo in - let center = Co.(expo_inv |. (p_a_expo_center |+ p_b_expo_center ) ) - in - let argexpo = - (Cs.expo s_a).(i) *. (Cs.expo s_b).(j) *. norm_sq *. expo_inv - in - let g = - (pi *. expo_inv)**(1.5) *. exp (-. argexpo) - in - let coef = - norm_coef *. g - in - if abs_float coef < cutoff then - raise Null_contribution; - let center_a = - Co.(center |- Cs.center s_a) - in - let monocentric = - Cs.(center s_a = center s_b) - in - let totAngMomInt = - Am.(Cs.totAngMom s_a + Cs.totAngMom s_b) - |> Am.to_int - in -assert (expo= Psp.expo sp_ ); -assert (expo_inv= Psp.expo_inv sp_ ); -assert (center= Psp.center sp_ ); - Some ( (Cs.coef s_a).(i) *. (Cs.coef s_b).(j), { - Sp.i ; j ; - shell_a=s_a ; shell_b=s_b ; - coef ; - expo ; expo_inv ; - center ; center_a ; center_ab ; - norm_sq ; monocentric ; totAngMomInt - }) - with - | Null_contribution -> None - ) - ) + Array.mapi (fun i p_a -> + let c_a = (Cs.coef s_a).(i) in + let make = make i p_a in + Array.mapi (fun j p_b -> + let c_b = (Cs.coef s_b).(j) in + let coef = c_a *. c_b in + assert (coef <> 0.); + let cutoff = cutoff /. abs_float coef in + coef, make j p_b cutoff) (Cs.prim s_b)) (Cs.prim s_a) |> Array.to_list |> Array.concat |> Array.to_list - |> List.filter (function Some _ -> true | None -> false) - |> List.map (function Some x -> x | None -> assert false) + |> List.filter (function (_, Some _) -> true | _ -> false) + |> List.map (function (c, Some x) -> (c,x) | _ -> assert false) |> Array.of_list in - let coef = Array.map (fun (c,y) -> c *. y.Sp.coef) shell_pairs - and expo_inv = Array.map (fun (_,y) -> y.Sp.expo_inv) shell_pairs + + + let coef = Array.map (fun (c,y) -> c *. Psp.norm_coef y) shell_pairs + and expo_inv = Array.map (fun (_,y) -> Psp.expo_inv y) shell_pairs in let shell_pairs = Array.map snd shell_pairs in + let root = shell_pairs.(0) in { - shell_a = s_a ; shell_b = s_b ; coef ; expo_inv ; - shell_pairs ; center_ab=shell_pairs.(0).center_ab; - norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq; - totAngMomInt = shell_pairs.(0).Sp.totAngMomInt; + shell_a = s_a ; shell_b = s_b ; coef ; expo_inv ; shell_pairs ; + center_ab = Psp.a_minus_b root; + norm_coef_scale = Psp.norm_coef_scale root; + norm_sq=Psp.a_minus_b_sq root; + totAngMomInt = Psp.totAngMom root |> Am.to_int; } @@ -159,7 +79,7 @@ let norm_sq x = x.norm_sq let totAngMomInt x = x.totAngMomInt let norm_coef_scale x = x.norm_coef_scale -let monocentric x = x.shell_pairs.(0).Sp.monocentric +let monocentric x = Psp.monocentric x.shell_pairs.(0) (** Returns an integer characteristic of a contracted shell pair *) let hash a = @@ -199,8 +119,12 @@ let of_basis basis = let equivalent x y = (Array.length x = Array.length y) && - (Array.init (Array.length x) (fun k -> Sp.equivalent x.(k) y.(k)) - |> Array.fold_left (fun accu x -> x && accu) true) + let rec eqv = function + | 0 -> true + | k -> if Psp.equivalent x.(k) y.(k) then + eqv (k-1) + else false + in eqv (Array.length x - 1) (** A list of unique shell pairs *) diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli index fbf4cf9..15faa95 100644 --- a/Basis/ContractedShellPair.mli +++ b/Basis/ContractedShellPair.mli @@ -37,9 +37,9 @@ val shell_b : t -> ContractedShell.t build the contracted shell pair. *) -val shell_pairs : t -> ShellPair.t array -(** Returns an array of {!ShellPair.t}, containing all the pairs of primitive functions. - build the contracted shell pair. +val shell_pairs : t -> PrimitiveShellPair.t array +(** Returns an array of {!PrimitiveShellPair.t}, containing all the pairs of + primitive functions used to build the contracted shell pair. *) val coef : t -> float array diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 4dd344c..02fdcc2 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -8,7 +8,8 @@ module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair module Po = Powers -module Sp = ShellPair +module Ps = PrimitiveShell +module Psp = PrimitiveShellPair type t = Mat.t @@ -55,18 +56,15 @@ let contracted_class shell_a shell_b : float Zmap.t = if (abs_float coef_prod) > 1.e-4*.cutoff then begin let center_pa = - sp.(ab).Sp.center_a + Psp.center_minus_a sp.(ab) in let expo_inv = (Csp.expo_inv shell_p).(ab) in - let i, j = - sp.(ab).Sp.i, sp.(ab).Sp.j - in let expo_a = - (Cs.expo sp.(ab).Sp.shell_a).(i) + Ps.expo (Psp.shell_a sp.(ab)) and expo_b = - (Cs.expo sp.(ab).Sp.shell_b).(j) + Ps.expo (Psp.shell_b sp.(ab)) in let xyz_of_int k = diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 640f4a3..5c60cd4 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -8,7 +8,8 @@ module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair module Po = Powers -module Sp = ShellPair +module Ps = PrimitiveShell +module Psp = PrimitiveShellPair @@ -118,9 +119,9 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let norm_coef_scale_p = Csp.norm_coef_scale shell_p in - for ab=0 to (Array.length (Csp.shell_pairs shell_p) - 1) + let sp = Csp.shell_pairs shell_p in + for ab=0 to Array.length sp - 1 do - let b = (Csp.shell_pairs shell_p).(ab).Sp.j in try begin let coef_prod = (Csp.coef shell_p).(ab) in @@ -134,14 +135,18 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = (Csp.expo_inv shell_p).(ab) in + let expo_b = + Ps.expo (Psp.shell_b sp.(ab)) + in + let center_ab = Csp.center_ab shell_p in let center_p = - (Csp.shell_pairs shell_p).(ab).Sp.center + Psp.center sp.(ab) in let center_pa = - Co.(center_p |- Cs.center shell_a) + Psp.center_minus_a sp.(ab) in for c=0 to Array.length geometry - 1 do @@ -178,8 +183,8 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = hvrr_one_e angMomA angMomB zero_m_array - (Cs.expo shell_b).(b) - (Csp.expo_inv shell_p).(ab) + expo_b + expo_pq_inv center_ab center_pa center_pc map in diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 9f1a4b7..2f5fa64 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -10,7 +10,7 @@ module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair module Po = Powers -module Sp = ShellPair +module Psp = PrimitiveShellPair let cutoff = integrals_cutoff @@ -66,7 +66,7 @@ let contracted_class shell_a shell_b : float Zmap.t = (Csp.expo_inv shell_p).(ab) in let center_pa = - sp.(ab).Sp.center_a + Psp.center_minus_a sp.(ab) in Array.iteri (fun i key -> diff --git a/Basis/PrimitiveShellPair.ml b/Basis/PrimitiveShellPair.ml index 4099a43..3c037bd 100644 --- a/Basis/PrimitiveShellPair.ml +++ b/Basis/PrimitiveShellPair.ml @@ -153,3 +153,8 @@ let norm_coef x = x.norm_coef let expo x = x.expo let center x = x.center + +let shell_a x = x.shell_a + +let shell_b x = x.shell_b + diff --git a/Basis/PrimitiveShellPair.mli b/Basis/PrimitiveShellPair.mli index 2df93ee..d3bf483 100644 --- a/Basis/PrimitiveShellPair.mli +++ b/Basis/PrimitiveShellPair.mli @@ -111,3 +111,9 @@ val center_minus_a : t -> Coordinate.t val norm_coef : t -> float (** Normalization coefficient of the shell pair. *) +val shell_a : t -> PrimitiveShell.t +(** Returns the first primitive shell that was used to build the shell pair. *) + +val shell_b : t -> PrimitiveShell.t +(** Returns the second primitive shell that was used to build the shell pair. *) + diff --git a/Basis/ShellPair.ml b/Basis/ShellPair.ml deleted file mode 100644 index 80b44d3..0000000 --- a/Basis/ShellPair.ml +++ /dev/null @@ -1,40 +0,0 @@ -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 *) - coef : float; (* norm_coef * coef_a * coef_b * g, with - g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *) - totAngMomInt : int ; - i : int; - j : int; - shell_a : ContractedShell.t; - shell_b : ContractedShell.t; - monocentric : bool -} - - -(** Returns an integer characteristic of a primitive shell pair *) -let hash a = - Hashtbl.hash a - -let equivalent a b = - a = b -(* - Hashtbl.hash (a.expo, a.center_a, a.center_ab, a.coef, ContractedShell.totAngMom a.shell_a, ContractedShell.totAngMom a.shell_b) -*) - - -(** Comparison function, used for sorting *) -let cmp a b = - hash a - hash b - - - - - diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index fd6710b..fcdf0b0 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -4,8 +4,9 @@ module Am = AngularMomentum module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair -module Sp = ShellPair module Po = Powers +module Psp = PrimitiveShellPair +module Ps = PrimitiveShell let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff @@ -279,7 +280,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q and shell_c = Csp.shell_a shell_q and shell_d = Csp.shell_b shell_q and sp = Csp.shell_pairs shell_p - and sq = Csp.shell_pairs shell_q in let maxm = Csp.totAngMomInt shell_p + Csp.totAngMomInt shell_q in @@ -301,38 +301,38 @@ 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 norm_coef_scale_p = Csp.norm_coef_scale shell_p in + let norm_coef_scale_p_list = Array.to_list (Csp.norm_coef_scale shell_p) in let norm_coef_scale_q = Csp.norm_coef_scale shell_q in + for ab=0 to (Array.length sp - 1) do - let cab = (Csp.coef shell_p).(ab) in - let b = (Csp.shell_pairs shell_p).(ab).Sp.j in - let expo_b = (Cs.expo shell_b).(b) in - let expo_inv_p = (Csp.expo_inv shell_p).(ab) in - let center_ab = sp.(ab).Sp.center_ab in - let center_pa = sp.(ab).Sp.center_a in + + let sp_ab = (Csp.shell_pairs shell_p).(ab) in + let c_ab = (Csp.coef shell_p).(ab) in + let expo_b = Ps.expo (Psp.shell_b sp_ab) in + let expo_inv_p = Psp.expo_inv sp_ab in + let center_ab = Psp.a_minus_b sp_ab in + let center_pa = Psp.center_minus_a sp_ab in + for cd=0 to (Array.length (Csp.shell_pairs shell_q) - 1) do - let coef_prod = - cab *. (Csp.coef shell_q).(cd) - in + + let sp_cd = (Csp.shell_pairs shell_q).(cd) in + let c_cd = (Csp.coef shell_q).(cd) in + let coef_prod = c_ab *. c_cd in + (** Screening on the product of coefficients *) try if (abs_float coef_prod) < 1.e-3 *. cutoff then raise NullQuartet; - let center_pq = - Co.(sp.(ab).Sp.center |- sq.(cd).Sp.center) - in - let norm_pq_sq = - Co.dot center_pq center_pq - in - - let expo_inv_q = (Csp.expo_inv shell_q).(cd) in - + let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in + let norm_pq_sq = Co.dot center_pq center_pq in + let expo_inv_q = Psp.expo_inv sp_cd in let expo_pq_inv = expo_inv_p +. expo_inv_q in let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in + begin match Cs.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d) with @@ -342,16 +342,14 @@ 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 = (Csp.shell_pairs shell_q).(cd).Sp.j in - let expo_d = (Cs.expo shell_d).(d) in + let expo_d = Ps.expo (Psp.shell_b sp_cd) in let map_1d = Zmap.create (4*maxm) in let map_2d = Zmap.create (Array.length class_indices) in - let center_cd = sq.(cd).Sp.center_ab in - let center_qc = sq.(cd).Sp.center_a in + let center_cd = Psp.a_minus_b sp_cd in + let center_qc = Psp.center_minus_a sp_cd in let norm_coef_scale = - Array.to_list norm_coef_scale_p - |> List.map (fun v1 -> - Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) + List.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) + norm_coef_scale_p_list |> Array.concat in diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 3a7450e..745e0ae 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -6,8 +6,9 @@ module Am = AngularMomentum module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair -module Sp = ShellPair module Po = Powers +module Psp = PrimitiveShellPair +module Ps = PrimitiveShell exception NullQuartet exception Found @@ -629,9 +630,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q begin try let expo_inv_p = - Vec.init np (fun ab -> sp.(ab-1).Sp.expo_inv) + Vec.init np (fun ab -> Psp.expo_inv sp.(ab-1)) and expo_inv_q = - Vec.init nq (fun cd -> sq.(cd-1).Sp.expo_inv) + Vec.init nq (fun cd -> Psp.expo_inv sq.(cd-1)) in let coef = @@ -652,7 +653,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let center_pq = - Co.(sp.(i-1).Sp.center |- sq.(j-1).Sp.center) + Co.(Psp.center sp.(i-1) |- Psp.center sq.(j-1)) in let norm_pq_sq = Co.dot center_pq center_pq @@ -677,15 +678,15 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let expo_inv_p = - Array.map (fun shell_ab -> shell_ab.Sp.expo_inv) sp + Array.map (fun shell_ab -> Psp.expo_inv shell_ab) sp and expo_inv_q = - Array.map (fun shell_cd -> shell_cd.Sp.expo_inv) sq + Array.map (fun shell_cd -> Psp.expo_inv shell_cd) sq in let expo_b = - Array.map (fun shell_ab -> (Cs.expo shell_b).(shell_ab.Sp.j)) sp + Array.map (fun shell_ab -> Ps.expo (Psp.shell_b shell_ab) ) sp and expo_d = - Array.map (fun shell_cd -> (Cs.expo shell_d).(shell_cd.Sp.j)) sq + Array.map (fun shell_cd -> Ps.expo (Psp.shell_b shell_cd) ) sq in let norm_coef_scale_p = Csp.norm_coef_scale shell_p in @@ -698,7 +699,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let shell_cd = sq.(cd) in let cpq = - Co.(shell_ab.Sp.center |- shell_cd.Sp.center) + Co.(Psp.center shell_ab |- Psp.center shell_cd) in match xyz with | 0 -> Co.get X cpq; @@ -718,7 +719,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.init np (fun ab -> let shell_ab = sp.(ab) in let cpa = - Co.(shell_ab.Sp.center |- Cs.center shell_a) + Co.(Psp.center shell_ab |- Cs.center shell_a) in match xyz with | 0 -> Co.(get X cpa); @@ -737,7 +738,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.init nq (fun cd -> let shell_cd = sq.(cd) in let cqc = - Co.(shell_cd.Sp.center |- Cs.center shell_c) + Co.(Psp.center shell_cd |- Cs.center shell_c) in match xyz with | 0 -> Co.(get X cqc);