diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli index 84a7eb3..5ebe8dd 100644 --- a/Basis/ContractedShell.mli +++ b/Basis/ContractedShell.mli @@ -1,7 +1,7 @@ (** Set of contracted Gaussians with a given {!AngularMomentum.t} {% \\[ -(x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} \sum_{i=1}^{m} \mathcal{N}_i f_i d_i \exp \left( -\alpha_i |r-R_A|^2 \right) +\chi(r) = (x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} \sum_{i=1}^{m} \mathcal{N}_i f_i d_i \exp \left( -\alpha_i |r-R_A|^2 \right) \\] %} where: diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 3f70367..c34b105 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -14,7 +14,6 @@ type t = norm_sq : float; (* |A-B|^2 *) norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *) totAngMomInt : int; (* Total angular Momentum *) - monocentric : bool; } @@ -118,10 +117,21 @@ let create ?cutoff p_a p_b = 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; - monocentric = shell_pairs.(0).Sp.monocentric; } +let shell_a x = x.shell_a +let shell_b x = x.shell_b +let shell_pairs x = x.shell_pairs +let coef x = x.coef +let expo_inv x = x.expo_inv +let center_ab x = x.center_ab +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 + (** Returns an integer characteristic of a contracted shell pair *) let hash a = Array.map Hashtbl.hash a @@ -150,7 +160,7 @@ let cmp a b = (** The array of all shell pairs with their correspondance in the list of contracted shells. *) -let shell_pairs basis = +let of_basis basis = Array.mapi (fun i shell_a -> Array.mapi (fun j shell_b -> create shell_a shell_b) diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli new file mode 100644 index 0000000..4ec7760 --- /dev/null +++ b/Basis/ContractedShellPair.mli @@ -0,0 +1,68 @@ +(** A datastructure to represent pairs of contracted shells. + +A contracted shell pair is a product of two {!ContractedShell.t}: + +{% \\[ +\varphi_{ab}(r) & = \chi_a(r) \times \chi_b(r) \\ + & = P_a \sum_{i=1}^m_a \mathcal{N}_i^a f_i[P_a] d_i^a \exp \left( -\alpha_i^a |r-R_a|^2 \right) \times \\ + & P_b \sum_{j=1}^m_b \mathcal{N}_j^b f_j[P_b] d_j^b \exp \left( -\alpha_j^b |r-R_b|^2 \right) + & = P_a P_b \sum_{i=1}^m_a \sum_{j=1}^m_b \mathcal{N}_i^a \mathcal{N}_j^b f_i[P_a] f_j[P_b] d_i^a d_j^b \exp \left( -\alpha_i^a |r-R_a|^2 \right -\alpha_j^b |r-R_b|^2 \right) + \\] + +*) + +type t + + +val create : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t +(** Creates an contracted shell pair {% $\varphi_{ab}$ %} from a contracted + shell {% $\chi_a$ %} (first argument) and a contracted shell {% $\chi_b$ %} + (second argument). + + The contracted shell pair contains the pairs of primitives for which + the norm is greater than [cutoff], and for which ... TODO .... + All other pairs are discarded. +*) + +val of_basis : ContractedShell.t array -> t array array +(** Creates all possible contracted shell pairs from the basis set. *) + +val shell_a : t -> ContractedShell.t +(** Returns the first {%ContractedShell.t} {% $\chi_a$ %} which was used to + build the contracted shell pair. +*) + +val shell_b : t -> ContractedShell.t +(** Returns the second {%ContractedShell.t} {% $\chi_b$ %} which was used to + 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 coef : t -> float array + +val expo_inv : t -> float array + +val center_ab : t -> Coordinate.t + (* A-B *) + +val norm_sq : t -> float + (* |A-B|^2 *) + +val norm_coef_scale : t -> float array + (* norm_coef.(i) / norm_coef.(0) *) + +val totAngMomInt : t -> int + (* Total angular Momentum *) + +val monocentric : t -> bool +(** If true, the two contracted shells have the same center. *) + +val hash : 'a array -> int array +val cmp : 'a array -> 'a array -> int +val equivalent : 'a array -> 'a array -> bool +val unique : 'a array array array -> 'a array list +val indices : 'a array array array -> (int array, int * int) Hashtbl.t diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 3ebad23..0dc0cdd 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -91,7 +91,7 @@ let of_basis basis = (* Pre-compute all shell pairs *) let shell_pairs = - Csp.shell_pairs shell + Csp.of_basis shell in (* Pre-compute diagonal integrals for Schwartz *) @@ -143,7 +143,7 @@ let of_basis basis = in let sp = - shell_p.Csp.shell_pairs + Csp.shell_pairs shell_p in for k=0 to i do @@ -156,7 +156,7 @@ let of_basis basis = shell_q = shell_pairs.(k).(l) in let sq = - shell_q.Csp.shell_pairs + Csp.shell_pairs shell_q in let swap = diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 18025d9..4dd344c 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -38,18 +38,18 @@ let contracted_class shell_a shell_b : float Zmap.t = (* Compute all integrals in the shell for each pair of significant shell pairs *) - let sp = shell_p.Csp.shell_pairs in + let sp = Csp.shell_pairs shell_p in let center_ab = - shell_p.Csp.center_ab + Csp.center_ab shell_p in let norm_coef_scale = - shell_p.Csp.norm_coef_scale + Csp.norm_coef_scale shell_p in for ab=0 to (Array.length sp - 1) do let coef_prod = - shell_p.Csp.coef.(ab) + (Csp.coef shell_p).(ab) in (** Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-4*.cutoff then @@ -58,7 +58,7 @@ let contracted_class shell_a shell_b : float Zmap.t = sp.(ab).Sp.center_a in let expo_inv = - shell_p.Csp.expo_inv.(ab) + (Csp.expo_inv shell_p).(ab) in let i, j = sp.(ab).Sp.i, sp.(ab).Sp.j diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 1592de1..a7d800c 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -97,8 +97,8 @@ let hvrr_one_e angMom_a angMom_b (** 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.Csp.shell_a - and shell_b = shell_p.Csp.shell_b + let shell_a = Csp.shell_a shell_p + and shell_b = Csp.shell_b shell_p in let maxm = Am.to_int (Cs.totAngMom shell_a) + Am.to_int (Cs.totAngMom shell_b) @@ -115,14 +115,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 *) - let norm_coef_scale_p = shell_p.Csp.norm_coef_scale + let norm_coef_scale_p = Csp.norm_coef_scale shell_p in - for ab=0 to (Array.length shell_p.Csp.shell_pairs - 1) + for ab=0 to (Array.length (Csp.shell_pairs shell_p) - 1) do - let b = shell_p.Csp.shell_pairs.(ab).Sp.j in + let b = (Csp.shell_pairs shell_p).(ab).Sp.j in try begin - let coef_prod = shell_p.Csp.coef.(ab) in + let coef_prod = (Csp.coef shell_p).(ab) in (** Screening on the product of coefficients *) if abs_float coef_prod < 1.e-3 *. integrals_cutoff then @@ -130,14 +130,14 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let expo_pq_inv = - shell_p.Csp.expo_inv.(ab) + (Csp.expo_inv shell_p).(ab) in let center_ab = - shell_p.Csp.center_ab + Csp.center_ab shell_p in let center_p = - shell_p.Csp.shell_pairs.(ab).Sp.center + (Csp.shell_pairs shell_p).(ab).Sp.center in let center_pa = Co.(center_p |- Cs.center shell_a) @@ -178,7 +178,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = angMomA angMomB zero_m_array (Cs.expo shell_b).(b) - shell_p.Csp.expo_inv.(ab) + (Csp.expo_inv shell_p).(ab) center_ab center_pa center_pc map in diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index c787b38..9f1a4b7 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -37,13 +37,13 @@ let contracted_class shell_a shell_b : float Zmap.t = in let sp = - shell_p.Csp.shell_pairs + Csp.shell_pairs shell_p in let center_ab = - shell_p.Csp.center_ab + Csp.center_ab shell_p in let norm_coef_scale = - shell_p.Csp.norm_coef_scale + Csp.norm_coef_scale shell_p in (* Compute all integrals in the shell for each pair of significant shell pairs *) @@ -57,13 +57,13 @@ let contracted_class shell_a shell_b : float Zmap.t = for ab=0 to (Array.length sp - 1) do let coef_prod = - shell_p.Csp.coef.(ab) + (Csp.coef shell_p).(ab) in (** Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-3*.cutoff then begin let expo_inv = - shell_p.Csp.expo_inv.(ab) + (Csp.expo_inv shell_p).(ab) in let center_pa = sp.(ab).Sp.center_a diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 44bcb8d..6645c24 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -274,14 +274,14 @@ let rec hvrr_two_e let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - let shell_a = shell_p.Csp.shell_a - and shell_b = shell_p.Csp.shell_b - and shell_c = shell_q.Csp.shell_a - and shell_d = shell_q.Csp.shell_b - and sp = shell_p.Csp.shell_pairs - and sq = shell_q.Csp.shell_pairs + let shell_a = Csp.shell_a shell_p + 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 + and sq = Csp.shell_pairs shell_q in - let maxm = shell_p.Csp.totAngMomInt + shell_q.Csp.totAngMomInt in + let maxm = Csp.totAngMomInt shell_p + Csp.totAngMomInt shell_q in (* Pre-computation of integral class indices *) let class_indices = @@ -295,22 +295,20 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let monocentric = - shell_p.Csp.monocentric && - shell_q.Csp.monocentric && - Cs.center shell_p.Csp.shell_a = - Cs.center shell_q.Csp.shell_a + Csp.monocentric shell_p && Csp.monocentric shell_q && + Cs.center (Csp.shell_a shell_p) = Cs.center (Csp.shell_a shell_q) in (* Compute all integrals in the shell for each pair of significant shell pairs *) - let norm_coef_scale_p = shell_p.Csp.norm_coef_scale in - let norm_coef_scale_q = shell_q.Csp.norm_coef_scale in + let norm_coef_scale_p = 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 = shell_p.Csp.coef.(ab) in + let cab = (Csp.coef shell_p).(ab) in let b = sp.(ab).Sp.j in - for cd=0 to (Array.length shell_q.Csp.shell_pairs - 1) do + for cd=0 to (Array.length (Csp.shell_pairs shell_q) - 1) do let coef_prod = - cab *. shell_q.Csp.coef.(cd) + cab *. (Csp.coef shell_q).(cd) in (** Screening on the product of coefficients *) try @@ -325,8 +323,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let expo_pq_inv = - shell_p.Csp.expo_inv.(ab) +. - shell_q.Csp.expo_inv.(cd) + (Csp.expo_inv shell_p).(ab) +. + (Csp.expo_inv shell_q).(cd) in let zero_m_array = @@ -341,7 +339,7 @@ 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.Csp.shell_pairs.(cd).Sp.j in + let d = (Csp.shell_pairs shell_q).(cd).Sp.j in let map_1d = Zmap.create (4*maxm) in let map_2d = Zmap.create (Array.length class_indices) in let norm_coef_scale = @@ -403,8 +401,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q angMom_a angMom_b angMom_c angMom_d zero_m_array (Cs.expo shell_b).(b) (Cs.expo shell_d).(d) - shell_p.Csp.expo_inv.(ab) - shell_q.Csp.expo_inv.(cd) + (Csp.expo_inv shell_p).(ab) + (Csp.expo_inv shell_q).(cd) sp.(ab).Sp.center_ab sq.(cd).Sp.center_ab center_pq sp.(ab).Sp.center_a sq.(cd).Sp.center_a map_1d map_2d diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 49b65c8..816f3e4 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -554,16 +554,15 @@ 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.Csp.shell_a - and shell_b = shell_p.Csp.shell_b - and shell_c = shell_q.Csp.shell_a - and shell_d = shell_q.Csp.shell_b - and sp = shell_p.Csp.shell_pairs - and sq = shell_q.Csp.shell_pairs + let shell_a = Csp.shell_a shell_p + 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 + and sq = Csp.shell_pairs shell_q in let maxm = - shell_p.Csp.totAngMomInt + - shell_q.Csp.totAngMomInt + Csp.totAngMomInt shell_p + Csp.totAngMomInt shell_q in (* Pre-computation of integral class indices *) @@ -578,21 +577,20 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let monocentric = - shell_p.Csp.monocentric && - shell_q.Csp.monocentric && - Cs.center shell_p.Csp.shell_a = - Cs.center shell_q.Csp.shell_a + Csp.monocentric shell_p && Csp.monocentric shell_q && + Cs.center (Csp.shell_a shell_p) = + Cs.center (Csp.shell_a shell_q) in (** Screening on the product of coefficients *) let coef_max_p = Array.fold_left (fun accu x -> if (abs_float x) > accu then (abs_float x) else accu) - 0. shell_p.Csp.coef + 0. (Csp.coef shell_p) and coef_max_q = Array.fold_left (fun accu x -> if (abs_float x) > accu then (abs_float x) else accu) - 0. shell_q.Csp.coef + 0. (Csp.coef shell_q) in let rec build_list cutoff vec accu = function @@ -602,10 +600,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q else accu ) (k-1) in let p_list = - let vec = shell_p.Csp.coef in + let vec = (Csp.coef shell_p) in build_list (cutoff /. coef_max_q) vec [] (Array.length vec - 1) and q_list = - let vec = shell_q.Csp.coef in + let vec = (Csp.coef shell_q) in build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1) in @@ -639,8 +637,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let coef = let result = Mat.make0 nq np in Lacaml.D.ger - (Vec.of_array @@ filter_q shell_q.Csp.coef) - (Vec.of_array @@ filter_p shell_p.Csp.coef) + (Vec.of_array @@ filter_q (Csp.coef shell_q)) + (Vec.of_array @@ filter_p (Csp.coef shell_p)) result; result in @@ -672,8 +670,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q | _ -> let coef = - let cp = filter_p shell_p.Csp.coef - and cq = filter_q shell_q.Csp.coef + let cp = filter_p (Csp.coef shell_p) + and cq = filter_q (Csp.coef shell_q) in Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) in @@ -689,7 +687,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q and expo_d = 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 + let norm_coef_scale_p = Csp.norm_coef_scale shell_p in let center_pq = let result = @@ -790,7 +788,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm = let norm_coef_scale_q = - shell_q.Csp.norm_coef_scale + Csp.norm_coef_scale shell_q in Array.to_list norm_coef_scale_p |> List.map (fun v1 -> @@ -846,8 +844,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q zero_m_array (expo_b, expo_d) (expo_inv_p, expo_inv_q) - (shell_p.Csp.center_ab, - shell_q.Csp.center_ab, center_pq) + (Csp.center_ab shell_p, + Csp.center_ab shell_q, center_pq) (center_pa, center_qc) map_1d map_2d np nq in