From 786b680f401bbbf924fc543724b562a2d20fdf7f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 13 Mar 2018 18:56:28 +0100 Subject: [PATCH] ContractedShell has now a hidden type --- Basis/Basis.ml | 6 ++-- Basis/ContractedShell.ml | 32 +++++++++++++++------ Basis/ContractedShell.mli | 48 ++++++++++++++++++++++---------- Basis/ContractedShellPair.ml | 32 ++++++++++----------- Basis/ERI.ml | 20 ++++++------- Basis/KinInt.ml | 14 +++++----- Basis/NucInt.ml | 8 +++--- Basis/OneElectronRR.ml | 10 +++---- Basis/Overlap.ml | 10 +++---- Basis/TwoElectronRR.ml | 14 +++++----- Basis/TwoElectronRRVectorized.ml | 21 +++++++------- 11 files changed, 125 insertions(+), 90 deletions(-) diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 2ad62bd..38c7741 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -27,14 +27,14 @@ let of_nuclei_and_general_basis n b = Array.iteri (fun i x -> if (i > 0) then result.(i) <- Cs.with_index x ( - result.(i-1).index + - (Array.length result.(i-1).norm_coef_scale )) + Cs.index result.(i-1) + + Cs.size_of_shell result.(i-1) ) ) result ; let size = let n = ref 0 in for i=0 to (Array.length result) - 1 do - n := !n + (Array.length (result.(i).norm_coef_scale)) + n := !n + Cs.size_of_shell result.(i) done; !n in { contracted_shells = result ; size } diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml index 168958c..ba16751 100644 --- a/Basis/ContractedShell.ml +++ b/Basis/ContractedShell.ml @@ -3,14 +3,14 @@ open Constants open Coordinate type t = { - expo : float array; - coef : float array; - center : Coordinate.t; - totAngMom : AngularMomentum.t; - size : int; - norm_coef : float array; - norm_coef_scale : float array; - index : int; + expo : float array; (** Array of exponents {% $\alpha_i$ %} *) + coef : float array; (** Array of contraction coefficients {% $d_i$ %} *) + center : Coordinate.t; (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %} *) + totAngMom : AngularMomentum.t; (** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *) + size : int; (** Number of contracted functions, {% $m$ %} in the formula *) + norm_coef : float array; (** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *) + norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *) + index : int; (** Index in the basis set, represented as an array of contracted shells. *) } module Am = AngularMomentum @@ -113,4 +113,20 @@ let compute_norm_coef expo totAngMom = Array.map (fun x -> let f a = x *. factor a in f) expo +let expo x = x.expo +let coef x = x.coef + +let center x = x.center + +let totAngMom x = x.totAngMom + +let size x = Array.length x.coef + +let norm_coef x = x.norm_coef + +let norm_coef_scale x = x.norm_coef_scale + +let index x = x.index + +let size_of_shell x = Array.length x.norm_coef_scale diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli index cacf36c..84a7eb3 100644 --- a/Basis/ContractedShell.mli +++ b/Basis/ContractedShell.mli @@ -28,28 +28,46 @@ f_i = \frac{1}{\mathcal{N}_i} *) -type t = private { - expo : float array; (** Array of exponents {% $\alpha_i$ %} *) - coef : float array; (** Array of contraction coefficients {% $d_i$ %} *) - center : Coordinate.t; (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %} *) - totAngMom : AngularMomentum.t; (** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *) - size : int; (** Number of contracted functions, {% $m$ %} in the formula *) - norm_coef : float array; (** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *) - norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *) - index : int; (** Index in the basis set, represented as an array of contracted shells. *) -} +type t - - -(** Pretty-printing of the contracted shell in a string *) val to_string : t -> string +(** Pretty-printing of the contracted shell in a string *) -(** Creates a contracted shell *) val make : index:int -> expo:float array -> coef:float array -> center:Coordinate.t -> totAngMom:AngularMomentum.t -> t +(** Creates a contracted shell *) -(** Returns a copy of the contracted shell with a modified index *) val with_index : t -> int -> t +(** Returns a copy of the contracted shell with a modified index *) + + +val expo : t -> float array +(** Array of exponents {% $\alpha_i$ %} *) + +val coef : t -> float array +(** Array of contraction coefficients {% $d_i$ %} *) + +val center : t -> Coordinate.t +(** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %} *) + +val totAngMom : t -> AngularMomentum.t +(** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *) + +val size : t -> int +(** Number of contracted functions, {% $m$ %} in the formula *) + +val norm_coef : t -> float array +(** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *) + +val norm_coef_scale : t -> float array +(** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *) + +val index : t -> int +(** Index in the basis set, represented as an array of contracted shells. *) + +val size_of_shell : t -> int +(** Number of contracted functions in the shell *) + diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index f7f43a8..3f70367 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -34,15 +34,15 @@ let create ?cutoff p_a p_b = | Some cutoff -> cutoff, -. (log cutoff) in - let center_ab = Co.( p_a.Cs.center |- p_b.Cs.center ) + let center_ab = Co.( Cs.center p_a |- Cs.center p_b ) in let norm_sq = Co.dot center_ab center_ab in let norm_coef_scale_a = - p_a.norm_coef_scale + Cs. norm_coef_scale p_a and norm_coef_scale_b = - p_b.norm_coef_scale + Cs. norm_coef_scale p_b in let norm_coef_scale = Array.map (fun v1 -> @@ -52,24 +52,24 @@ let create ?cutoff p_a p_b = |> Array.concat in let shell_pairs = - Array.init p_a.size (fun i -> - let p_a_expo_center = Co.(p_a.Cs.expo.(i) |. p_a.Cs.center) in - let norm_coef_a = p_a.norm_coef.(i) in + Array.init (Cs.size p_a) (fun i -> + let p_a_expo_center = Co.( (Cs.expo p_a).(i) |. Cs.center p_a ) in + let norm_coef_a = (Cs.norm_coef p_a).(i) in - Array.init p_b.size (fun j -> + Array.init (Cs.size p_b) (fun j -> try - let norm_coef_b = p_b.norm_coef.(j) in + let norm_coef_b = (Cs.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 = Co.(p_b.expo.(j) |. p_b.center) in - let expo = p_a.expo.(i) +. p_b.expo.(j) in + let p_b_expo_center = Co.( (Cs.expo p_b).(j) |. Cs.center p_b ) in + let expo = (Cs.expo p_a).(i) +. (Cs.expo p_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 = - p_a.Cs.expo.(i) *. p_b.Cs.expo.(j) *. norm_sq *. expo_inv + (Cs.expo p_a).(i) *. (Cs.expo p_b).(j) *. norm_sq *. expo_inv in if (argexpo > log_cutoff) then raise Null_contribution; @@ -77,19 +77,19 @@ let create ?cutoff p_a p_b = (pi *. expo_inv)**(1.5) *. exp (-. argexpo) in let coef = - norm_coef *. p_a.coef.(i) *. p_b.coef.(j) *. g + norm_coef *. (Cs.coef p_a).(i) *. (Cs.coef p_b).(j) *. g in if abs_float coef < cutoff then raise Null_contribution; let center_a = - Co.(center |- p_a.center) + Co.(center |- Cs.center p_a) in let monocentric = - p_a.center = p_b.center + Cs.(center p_a = center p_b) in let totAngMomInt = - Am.to_int p_a.totAngMom + - Am.to_int p_b.totAngMom + Am.to_int (Cs.totAngMom p_a) + + Am.to_int (Cs.totAngMom p_b) in Some { Sp.i ; j ; diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 81b0ade..3ebad23 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -108,7 +108,7 @@ let of_basis basis = let icount = ref 0 in for i=0 to (Array.length shell) - 1 do - print_int shell.(i).index ; print_newline (); + print_int (Cs.index shell.(i)) ; print_newline (); for j=0 to i do let schwartz_p, schwartz_p_max = schwartz.(i).(j) in if (schwartz_p_max >= cutoff) then @@ -132,7 +132,7 @@ let of_basis basis = let inn = ref 0 and out = ref 0 in for i=0 to (Array.length shell) - 1 do - print_int shell.(i).index ; print_newline (); + print_int (Cs.index shell.(i)) ; print_newline (); for j=0 to i do let schwartz_p, schwartz_p_max = schwartz.(i).(j) in try @@ -180,16 +180,16 @@ let of_basis basis = (* Write the data in the output file *) Array.iteri (fun i_c powers_i -> - let i_c = shell.(i).index + i_c + 1 in + let i_c = Cs.index shell.(i) + i_c + 1 in let xi = to_powers powers_i in Array.iteri (fun j_c powers_j -> - let j_c = shell.(j).index + j_c + 1 in + let j_c = Cs.index shell.(j) + j_c + 1 in let xj = to_powers powers_j in Array.iteri (fun k_c powers_k -> - let k_c = shell.(k).index + k_c + 1 in + let k_c = Cs.index shell.(k) + k_c + 1 in let xk = to_powers powers_k in Array.iteri (fun l_c powers_l -> - let l_c = shell.(l).index + l_c + 1 in + let l_c = Cs.index shell.(l) + l_c + 1 in let xl = to_powers powers_l in let key = if swap then @@ -213,10 +213,10 @@ let of_basis basis = ) else out := !out + 1; - ) Am.(zkey_array (Singlet shell.(l).Cs.totAngMom)) - ) Am.(zkey_array (Singlet shell.(k).Cs.totAngMom)) - ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) - ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) + ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(l)))) + ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(k)))) + ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(j)))) + ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(i)))) with NullIntegral -> () done; done; diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index a65be7b..18025d9 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -29,7 +29,7 @@ let contracted_class shell_a shell_b : float Zmap.t = (* Pre-computation of integral class indices *) let class_indices = - Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom)) + Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b)) in let contracted_class = @@ -64,9 +64,9 @@ let contracted_class shell_a shell_b : float Zmap.t = sp.(ab).Sp.i, sp.(ab).Sp.j in let expo_a = - sp.(ab).Sp.shell_a.expo.(i) + (Cs.expo sp.(ab).Sp.shell_a).(i) and expo_b = - sp.(ab).Sp.shell_b.expo.(j) + (Cs.expo sp.(ab).Sp.shell_b).(j) in let xyz_of_int k = @@ -141,10 +141,10 @@ let of_basis basis = in Array.iteri (fun j_c powers_j -> - let j_c = shell.(j).index + j_c + 1 in + let j_c = Cs.index shell.(j) + j_c + 1 in let xj = to_powers powers_j in Array.iteri (fun i_c powers_i -> - let i_c = shell.(i).index + i_c + 1 in + let i_c = Cs.index shell.(i) + i_c + 1 in let xi = to_powers powers_i in let key = Zkey.of_powers_six xi xj @@ -155,8 +155,8 @@ let of_basis basis = in result.{i_c,j_c} <- value; result.{j_c,i_c} <- value; - ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) - ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) + ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i)))) + ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) done; done; Mat.detri result; diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index 24ca506..e86b2f8 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -75,10 +75,10 @@ let of_basis_nuclei basis nuclei = (* Write the data in the output file *) Array.iteri (fun i_c powers_i -> - let i_c = shell.(i).index + i_c + 1 in + let i_c = Cs.index shell.(i) + i_c + 1 in let xi = to_powers powers_i in Array.iteri (fun j_c powers_j -> - let j_c = shell.(j).index + j_c + 1 in + let j_c = Cs.index shell.(j) + j_c + 1 in let xj = to_powers powers_j in let key = Zkey.of_powers_six xi xj @@ -88,8 +88,8 @@ let of_basis_nuclei basis nuclei = in eni_array.{j_c,i_c} <- value; eni_array.{i_c,j_c} <- value; - ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) - ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) + ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) + ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i)))) done; done; Mat.detri eni_array; diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 762242b..1592de1 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -101,12 +101,12 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = and shell_b = shell_p.Csp.shell_b in let maxm = - (Am.to_int @@ shell_a.totAngMom) + (Am.to_int @@ shell_b.totAngMom) + Am.to_int (Cs.totAngMom shell_a) + Am.to_int (Cs.totAngMom shell_b) in (* Pre-computation of integral class indices *) let class_indices = - Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom)) + Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b)) in let contracted_class = @@ -140,7 +140,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = shell_p.Csp.shell_pairs.(ab).Sp.center in let center_pa = - Co.(center_p |- shell_a.center) + Co.(center_p |- Cs.center shell_a) in for c=0 to Array.length geometry - 1 do @@ -156,7 +156,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in - match (shell_a.totAngMom, shell_b.totAngMom) with + match Cs.(totAngMom shell_a, totAngMom shell_b) with | Am.(S,S) -> let integral = zero_m_array.(0) in contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge @@ -177,7 +177,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = hvrr_one_e angMomA angMomB zero_m_array - shell_b.expo.(b) + (Cs.expo shell_b).(b) shell_p.Csp.expo_inv.(ab) center_ab center_pa center_pc map diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 48135ee..c787b38 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -29,7 +29,7 @@ let contracted_class shell_a shell_b : float Zmap.t = (* Pre-computation of integral class indices *) let class_indices = - Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom)) + Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b)) in let contracted_class = @@ -113,10 +113,10 @@ let of_basis basis = in Array.iteri (fun j_c powers_j -> - let j_c = shell.(j).index + j_c + 1 in + let j_c = Cs.index shell.(j) + j_c + 1 in let xj = to_powers powers_j in Array.iteri (fun i_c powers_i -> - let i_c = shell.(i).index + i_c + 1 in + let i_c = Cs.index shell.(i) + i_c + 1 in let xi = to_powers powers_i in let key = Zkey.of_powers_six xi xj @@ -127,8 +127,8 @@ let of_basis basis = in result.{i_c,j_c} <- value; result.{j_c,i_c} <- value; - ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) - ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) + ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i)))) + ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) done; done; Mat.detri result; diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 68e95b4..44bcb8d 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -286,8 +286,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (* Pre-computation of integral class indices *) let class_indices = Am.zkey_array (Am.Quartet - (shell_a.totAngMom, shell_b.totAngMom, - shell_c.totAngMom, shell_d.totAngMom )) + Cs.(totAngMom shell_a, totAngMom shell_b, + totAngMom shell_c, totAngMom shell_d )) in let contracted_class = @@ -297,8 +297,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let monocentric = shell_p.Csp.monocentric && shell_q.Csp.monocentric && - shell_p.Csp.shell_a.Cs.center = - shell_q.Csp.shell_a.Cs.center + Cs.center shell_p.Csp.shell_a = + Cs.center shell_q.Csp.shell_a in (* Compute all integrals in the shell for each pair of significant shell pairs *) @@ -333,8 +333,8 @@ 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 begin - match (shell_a.totAngMom, shell_b.totAngMom, - shell_c.totAngMom, shell_d.totAngMom) with + match Cs.(totAngMom shell_a, totAngMom shell_b, + totAngMom shell_c, totAngMom shell_d) with | Am.(S,S,S,S) -> let integral = zero_m_array.(0) @@ -402,7 +402,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q hvrr_two_e angMom_a angMom_b angMom_c angMom_d zero_m_array - shell_b.expo.(b) shell_d.expo.(d) + (Cs.expo shell_b).(b) (Cs.expo shell_d).(d) shell_p.Csp.expo_inv.(ab) shell_q.Csp.expo_inv.(cd) sp.(ab).Sp.center_ab sq.(cd).Sp.center_ab center_pq diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 7bb64f3..49b65c8 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -4,6 +4,7 @@ open Bigarray module Am = AngularMomentum module Co = Coordinate +module Cs = ContractedShell module Csp = ContractedShellPair module Sp = ShellPair module Po = Powers @@ -568,8 +569,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (* Pre-computation of integral class indices *) let class_indices = Am.zkey_array (Am.Quartet - (shell_a.totAngMom, shell_b.totAngMom, - shell_c.totAngMom, shell_d.totAngMom)) + Cs.(totAngMom shell_a, totAngMom shell_b, + totAngMom shell_c, totAngMom shell_d)) in let contracted_class = @@ -579,8 +580,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let monocentric = shell_p.Csp.monocentric && shell_q.Csp.monocentric && - shell_p.Csp.shell_a.center = - shell_q.Csp.shell_a.center + Cs.center shell_p.Csp.shell_a = + Cs.center shell_q.Csp.shell_a in (** Screening on the product of coefficients *) @@ -623,8 +624,8 @@ 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 *) begin - match (shell_a.totAngMom, shell_b.totAngMom, - shell_c.totAngMom, shell_d.totAngMom) with + match Cs.(totAngMom shell_a, totAngMom shell_b, + totAngMom shell_c, totAngMom shell_d) with | Am.(S,S,S,S) -> contracted_class.(0) <- begin @@ -684,9 +685,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let expo_b = - Array.map (fun shell_ab -> shell_b.expo.(shell_ab.Sp.j)) sp + Array.map (fun shell_ab -> (Cs.expo shell_b).(shell_ab.Sp.j)) sp and expo_d = - Array.map (fun shell_cd -> shell_d.expo.(shell_cd.Sp.j)) sq + Array.map (fun shell_cd -> (Cs.expo shell_d).(shell_cd.Sp.j)) sq in let norm_coef_scale_p = shell_p.Csp.norm_coef_scale in @@ -719,7 +720,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 |- shell_a.center) + Co.(shell_ab.Sp.center |- Cs.center shell_a) in match xyz with | 0 -> Co.(get X cpa); @@ -738,7 +739,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 |- shell_c.center) + Co.(shell_cd.Sp.center |- Cs.center shell_c) in match xyz with | 0 -> Co.(get X cqc);