From a50ffd83f80dbdfb8e3437bb2bb25e8024747634 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 20 Mar 2018 19:02:58 +0100 Subject: [PATCH] List of shell pairs --- Basis/ContractedShellPair.ml | 79 +++++++++++++++++------------- Basis/ContractedShellPair.mli | 6 +++ Basis/OneElectronRR.ml | 53 ++++++++------------- Basis/Overlap.ml | 90 ++++++++++++++++------------------- Basis/TwoElectronRR.ml | 22 ++++----- 5 files changed, 119 insertions(+), 131 deletions(-) diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 586b9a0..ece2116 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -7,13 +7,7 @@ type t = { shell_a : ContractedShell.t; shell_b : ContractedShell.t; - shell_pairs : PrimitiveShellPair.t list; - coefficients : float array; - exponents_inv : float array; - center_ab : Coordinate.t; (* A-B *) - norm_sq : float; (* |A-B|^2 *) - norm_scales : float array; (* norm_coef.(i) / norm_coef.(0) *) - totAngMom : AngularMomentum.t; (* Total angular Momentum *) + coefs_and_shell_pairs : (float * PrimitiveShellPair.t) list; } @@ -31,7 +25,7 @@ let make ?(cutoff=1.e-32) s_a s_b = let make = Psp.create_make_of (Cs.primitives s_a).(0) (Cs.primitives s_b).(0) in - let shell_pairs = + let coefs_and_shell_pairs = Array.mapi (fun i p_a -> let c_a = (Cs.coefficients s_a).(i) in let make = make p_a in @@ -45,38 +39,55 @@ let make ?(cutoff=1.e-32) s_a s_b = |> Array.concat |> Array.to_list |> List.filter (function (_, Some _) -> true | _ -> false) - |> List.map (function (c, Some x) -> (c,x) | _ -> assert false) + |> List.map (function (c, Some x) -> (c *. Psp.normalization x, x) | _ -> assert false) in - - match shell_pairs with + match coefs_and_shell_pairs with | [] -> None - | head :: _ -> - let coefficients = List.map (fun (c,y) -> c *. Psp.normalization y) shell_pairs |> Array.of_list - and exponents_inv = List.map (fun (_,y) -> Psp.exponent_inv y) shell_pairs |> Array.of_list - in - let shell_pairs = List.map snd shell_pairs in - let root = snd head in - Some { - shell_a = s_a ; shell_b = s_b ; coefficients ; exponents_inv ; shell_pairs ; - center_ab = Psp.a_minus_b root; - norm_scales = Psp.norm_scales root; - norm_sq=Psp.a_minus_b_sq root; - totAngMom = Psp.totAngMom root; - } + | coefs_and_shell_pairs -> Some { shell_a = s_a ; shell_b = s_b ; coefs_and_shell_pairs } -let shell_a x = x.shell_a -let shell_b x = x.shell_b -let shell_pairs x = Array.of_list x.shell_pairs -let coefficients x = x.coefficients -let exponents_inv x = x.exponents_inv -let center_ab x = x.center_ab -let norm_sq x = x.norm_sq -let totAngMom x = x.totAngMom -let norm_scales x = x.norm_scales +let shell_a x = x.shell_a +let shell_b x = x.shell_b +let coefs_and_shell_pairs x = x.coefs_and_shell_pairs + +let shell_pairs x = + List.map snd x.coefs_and_shell_pairs + |> Array.of_list + +let coefficients x = + List.map fst x.coefs_and_shell_pairs + |> Array.of_list + +let exponents_inv x = + List.map (fun (_,sp) -> Psp.exponent_inv sp) x.coefs_and_shell_pairs + |> Array.of_list + +let center_ab x = + match x.coefs_and_shell_pairs with + | [] -> assert false + | (_,sp)::_ -> Psp.a_minus_b sp + +let norm_sq x = + match x.coefs_and_shell_pairs with + | [] -> assert false + | (_,sp)::_ -> Psp.a_minus_b_sq sp + +let totAngMom x = + match x.coefs_and_shell_pairs with + | [] -> assert false + | (_,sp)::_ -> Psp.totAngMom sp + +let norm_scales x = + match x.coefs_and_shell_pairs with + | [] -> assert false + | (_,sp)::_ -> Psp.norm_scales sp + +let monocentric x = + match x.coefs_and_shell_pairs with + | [] -> assert false + | (_,sp)::_ -> Psp.monocentric sp -let monocentric x = Psp.monocentric (List.hd x.shell_pairs) (** Returns an integer characteristic of a contracted shell pair *) let hash a = diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli index 7afd25f..087352a 100644 --- a/Basis/ContractedShellPair.mli +++ b/Basis/ContractedShellPair.mli @@ -45,6 +45,11 @@ val shell_b : t -> ContractedShell.t build the contracted shell pair. *) +val coefs_and_shell_pairs : t -> (float * PrimitiveShellPair.t) list +(** Returns an list of coefficients and of {!PrimitiveShellPair.t}, containing all + the pairs of primitive functions used to 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. @@ -54,6 +59,7 @@ val coefficients : t -> float array val exponents_inv : t -> float array + val center_ab : t -> Coordinate.t (* A-B *) diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index b7fb00c..e1103c1 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -117,40 +117,24 @@ 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 *) - let norm_scales_p = Csp.norm_scales shell_p - in - let sp = Csp.shell_pairs shell_p in - for ab=0 to Array.length sp - 1 - do + let norm_scales_p = Csp.norm_scales shell_p in + + let center_ab = Csp.center_ab shell_p in + + List.iter (fun (coef_prod, psp) -> try begin - let coef_prod = (Csp.coefficients shell_p).(ab) in - (** Screening on the product of coefficients *) if abs_float coef_prod < 1.e-3 *. integrals_cutoff then raise NullPair; - - let expo_pq_inv = - (Csp.exponents_inv shell_p).(ab) + let expo_pq_inv = Psp.exponent_inv psp + and expo_b = Ps.exponent (Psp.shell_b psp) + and center_p = Psp.center psp + and center_pa = Psp.center_minus_a psp in - - let expo_b = - Ps.exponent (Psp.shell_b sp.(ab)) - in - - let center_ab = - Csp.center_ab shell_p - in - let center_p = - Psp.center sp.(ab) - in - let center_pa = - Psp.center_minus_a sp.(ab) - in - - for c=0 to Array.length geometry - 1 do - let element, nucl_coord = geometry.(c) in + + Array.iter (fun (element, nucl_coord) -> let charge = Element.to_charge element |> Charge.to_float in let center_pc = Co.(center_p |- nucl_coord ) @@ -190,15 +174,16 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = in contracted_class.(i) <- contracted_class.(i) -. coef_prod *. integral *. charge ) - done + ) geometry end with NullPair -> () -done; -let result = - Zmap.create (Array.length contracted_class) -in -Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; -result + ) (Csp.coefs_and_shell_pairs shell_p); + + let result = + Zmap.create (Array.length contracted_class) + in + Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; + result diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index a0a6225..5447a1b 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -36,61 +36,51 @@ let contracted_class shell_a shell_b : float Zmap.t = Array.make (Array.length class_indices) 0. in - let sp = - Csp.shell_pairs shell_p + let center_ab = + Csp.center_ab shell_p + in + let norm_coef_scales = + Csp.norm_scales shell_p in - let center_ab = - Csp.center_ab shell_p - in - let norm_coef_scales = - Csp.norm_scales shell_p - in + (* Compute all integrals in the shell for each pair of significant shell pairs *) - (* Compute all integrals in the shell for each pair of significant shell pairs *) + let xyz_of_int k = + match k with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in - let xyz_of_int k = - match k with - | 0 -> Co.X - | 1 -> Co.Y - | _ -> Co.Z - in - for ab=0 to (Array.length sp - 1) - do - let coef_prod = - (Csp.coefficients shell_p).(ab) - in - (** Screening on thr product of coefficients *) - if (abs_float coef_prod) > 1.e-3*.cutoff then - begin - let expo_inv = - (Csp.exponents_inv shell_p).(ab) - in - let center_pa = - Psp.center_minus_a sp.(ab) - in + List.iter (fun (coef_prod, psp) -> + (** Screening on thr product of coefficients *) + if (abs_float coef_prod) > 1.e-3*.cutoff then + begin + let expo_inv = Psp.exponent_inv psp + and center_pa = Psp.center_minus_a psp + in - Array.iteri (fun i key -> - let (angMomA,angMomB) = to_powers key in - let f k = - let xyz = xyz_of_int k in - Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) - expo_inv - (Co.get xyz center_ab, - Co.get xyz center_pa) - in - let norm = norm_coef_scales.(i) in - let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in - contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral - ) class_indices - end - done; - let result = - Zmap.create (Array.length contracted_class) - in - Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; - result - end + Array.iteri (fun i key -> + let (angMomA,angMomB) = to_powers key in + let f k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) + expo_inv + (Co.get xyz center_ab, + Co.get xyz center_pa) + in + let norm = norm_coef_scales.(i) in + let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in + contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral + ) class_indices + end + ) (Csp.coefs_and_shell_pairs shell_p); + let result = + Zmap.create (Array.length contracted_class) + in + Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; + result + end | None -> Zmap.create 0 diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index b736d3e..fdcbbb0 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -279,7 +279,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q and shell_b = Csp.shell_b shell_p and shell_c = Csp.shell_a shell_q and shell_d = Csp.shell_b shell_q - and sp = Csp.shell_pairs shell_p in let maxm = Am.(Csp.totAngMom shell_p + Csp.totAngMom shell_q |> to_int) in @@ -304,19 +303,16 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm_coef_scale_p_list = Array.to_list (Csp.norm_scales shell_p) in let norm_coef_scale_q = Csp.norm_scales shell_q in - for ab=0 to (Array.length sp - 1) do + let center_ab = Csp.center_ab shell_p in + List.iter (fun (c_ab, sp_ab) -> - let sp_ab = (Csp.shell_pairs shell_p).(ab) in - let c_ab = (Csp.coefficients shell_p).(ab) in - let expo_b = Ps.exponent (Psp.shell_b sp_ab) in - let expo_inv_p = Psp.exponent_inv sp_ab in - let center_ab = Psp.a_minus_b sp_ab in - let center_pa = Psp.center_minus_a sp_ab in + let expo_b = Ps.exponent (Psp.shell_b sp_ab) + and expo_inv_p = Psp.exponent_inv sp_ab + and center_pa = Psp.center_minus_a sp_ab + in - for cd=0 to (Array.length (Csp.shell_pairs shell_q) - 1) do + List.iter (fun (c_cd, sp_cd) -> - let sp_cd = (Csp.shell_pairs shell_q).(cd) in - let c_cd = (Csp.coefficients shell_q).(cd) in let coef_prod = c_ab *. c_cd in (** Screening on the product of coefficients *) @@ -415,8 +411,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q ) end with NullQuartet -> () - done - done; + ) (Csp.coefs_and_shell_pairs shell_q) + ) (Csp.coefs_and_shell_pairs shell_p); let result = Zmap.create (Array.length contracted_class)