diff --git a/Basis/AtomicShellPair.ml b/Basis/AtomicShellPair.ml index f133eaf..370b19c 100644 --- a/Basis/AtomicShellPair.ml +++ b/Basis/AtomicShellPair.ml @@ -15,12 +15,11 @@ module Am = AngularMomentum module As = AtomicShell module Co = Coordinate module Cs = ContractedShell -module Ps = PrimitiveShell -module Psp = PrimitiveShellPair +module Csp = ContractedShellPair (** Creates an atomic shell pair : an array of pairs of contracted shells. *) -let make ?(cutoff=1.e-32) atomic_shell_a atomic_shell_b = +let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b = let l_a = Array.to_list (As.contracted_shells atomic_shell_a) and l_b = Array.to_list (As.contracted_shells atomic_shell_b) @@ -32,27 +31,33 @@ let make ?(cutoff=1.e-32) atomic_shell_a atomic_shell_b = Csp.make ~cutoff s_a s_b ) l_b ) l_a + |> List.concat + |> List.filter (function None -> false | _ -> true) + |> List.map (function None -> assert false | Some x -> x) in + match contracted_shell_pairs with + | [] -> None + | _ -> Some { atomic_shell_a ; atomic_shell_b ; contracted_shell_pairs } let atomic_shell_a x = x.atomic_shell_a let atomic_shell_b x = x.atomic_shell_b -let contracted_shell_pairs = x.contracted_shell_pairs +let contracted_shell_pairs x = x.contracted_shell_pairs let monocentric x = Csp.monocentric @@ List.hd x.contracted_shell_pairs -let center_ab x = Csp.center_ab @@ List.hd x.contracted_shell_pairs +let a_minus_b x = Csp.a_minus_b @@ List.hd x.contracted_shell_pairs -let totAngMon x = Csp.totAngMon @@ List.hd x.contracted_shell_pairs +let a_minus_b_sq x = Csp.a_minus_b_sq @@ List.hd x.contracted_shell_pairs + +let ang_mom x = Csp.ang_mom @@ List.hd x.contracted_shell_pairs let norm_scales x = Csp.norm_scales @@ List.hd x.contracted_shell_pairs -let norm_sq x = Csp.norm_sq @@ List.hd x.contracted_shell_pairs - (** The array of all shell pairs with their correspondance in the list of contracted shells. *) -let of_atomic_shell_array basis = +let of_atomic_shell_array ?(cutoff=Constants.epsilon) basis = Array.mapi (fun i shell_a -> Array.mapi (fun j shell_b -> make ~cutoff shell_a shell_b) diff --git a/Basis/AtomicShellPair.mli b/Basis/AtomicShellPair.mli index 1c6cb88..674cca3 100644 --- a/Basis/AtomicShellPair.mli +++ b/Basis/AtomicShellPair.mli @@ -14,7 +14,7 @@ val make : ?cutoff:float -> AtomicShell.t -> AtomicShell.t -> t option [None]. *) -val of_atomic_shell_array : AtomicShell.t array -> t option list +val of_atomic_shell_array : ?cutoff:float -> AtomicShell.t array -> t option array array (** Creates all possible atomic shell pairs from an array of atomic shells. If an atomic shell pair is not significant, sets the value to [None]. *) @@ -34,10 +34,10 @@ val contracted_shell_pairs : t -> ContractedShellPair.t list contracted functions used to build the atomic shell pair. *) -val center_ab : t -> Coordinate.t +val a_minus_b : t -> Coordinate.t (* A-B *) -val norm_sq : t -> float +val a_minus_b_sq : t -> float (* |A-B|^2 *) val norm_scales : t -> float array diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 3fde3ef..b92732b 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -21,7 +21,7 @@ module Psp = PrimitiveShellPair A contracted shell with N functions combined with a contracted shell with M functions generates a NxM array of shell pairs. *) -let make ?(cutoff=1.e-32) s_a s_b = +let make ?(cutoff=Constants.epsilon) s_a s_b = let make = Psp.create_make_of (Cs.primitives s_a).(0) (Cs.primitives s_b).(0) in @@ -119,10 +119,10 @@ let cmp a b = (** The array of all shell pairs with their correspondance in the list of contracted shells. *) -let of_contracted_shell_array basis = +let of_contracted_shell_array ?(cutoff=Constants.epsilon) basis = Array.mapi (fun i shell_a -> Array.mapi (fun j shell_b -> - make shell_a shell_b) + make ~cutoff shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli index faccce9..f2c926d 100644 --- a/Basis/ContractedShellPair.mli +++ b/Basis/ContractedShellPair.mli @@ -27,7 +27,7 @@ val make : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option *) -val of_contracted_shell_array : ContractedShell.t array -> t option array array +val of_contracted_shell_array : ?cutoff:float -> ContractedShell.t array -> t option array array (** Creates all possible contracted shell pairs from a list of contracted shells. If a shell pair is not significant, sets the value to [None]: diff --git a/Basis/ERI.ml b/Basis/ERI.ml index bacbb65..e34d81e 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -42,20 +42,12 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = (** Compute all the integrals of a contracted class *) -(* -let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = - TwoElectronRRVectorized.contracted_class ~zero_m shell_a shell_b shell_c shell_d - *) -let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = - TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d - -(** Compute all the integrals of a contracted class *) -let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q - let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = TwoElectronRR.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q +let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = + TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q + let cutoff2 = cutoff *. cutoff (* diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 4d37f05..55bcecf 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -1,6 +1,8 @@ open Util module Am = AngularMomentum +module Asp = AtomicShellPair +module Aspc = AtomicShellPairCouple module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair @@ -408,15 +410,140 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q -(** Computes all the two-electron integrals of the contracted shell quartet *) -let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = - - let shell_p = Csp.make ~cutoff shell_a shell_b - and shell_q = Csp.make ~cutoff shell_c shell_d - in - match shell_p, shell_q with - | Some shell_p, Some shell_q -> - contracted_class_shell_pairs ~zero_m shell_p shell_q - | _ -> empty + + + +let contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = + + match Aspc.make shell_p shell_q with + | None -> empty + | Some atomic_shell_pair_couple -> + + let maxm = Am.to_int (Aspc.ang_mom atomic_shell_pair_couple) in + + (* Pre-computation of integral class indices *) + let class_indices = Aspc.zkey_array atomic_shell_pair_couple in + + let contracted_class = + Array.make (Array.length class_indices) 0.; + in + + let monocentric = + Aspc.monocentric atomic_shell_pair_couple + in + + (* Compute all integrals in the shell for each pair of significant shell pairs *) + + let center_ab = Asp.a_minus_b shell_p + and center_cd = Asp.a_minus_b shell_q + in + + let norm_scales = Aspc.norm_scales atomic_shell_pair_couple in + + + List.iter (fun cspc -> + List.iter (fun (coef_prod, spc) -> + let sp_ab = Pspc.shell_pair_p spc + and sp_cd = Pspc.shell_pair_q spc + in + + let expo_inv_p = Psp.exponent_inv sp_ab + 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.exponent_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 Aspc.ang_mom atomic_shell_pair_couple with + | Am.S -> + let integral = zero_m_array.(0) in + contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral + | _ -> + let expo_b = Ps.exponent (Psp.shell_b sp_ab) + and expo_d = Ps.exponent (Psp.shell_b sp_cd) + and center_pa = Psp.center_minus_a sp_ab + in + let map_1d = Zmap.create (4*maxm) + and map_2d = Zmap.create (Array.length class_indices) + in + let center_qc = Psp.center_minus_a sp_cd + in + + (* Compute the integral class from the primitive shell quartet *) + class_indices + |> Array.iteri (fun i key -> + let (angMom_a,angMom_b,angMom_c,angMom_d) = + match Zkey.to_powers key with + | Zkey.Twelve x -> x + | _ -> assert false + in + try + if monocentric then + begin + if ( ((1 land angMom_a.Po.x + angMom_b.Po.x + angMom_c.Po.x + angMom_d.Po.x)=1) || + ((1 land angMom_a.Po.y + angMom_b.Po.y + angMom_c.Po.y + angMom_d.Po.y)=1) || + ((1 land angMom_a.Po.z + angMom_b.Po.z + angMom_c.Po.z + angMom_d.Po.z)=1) + ) then + raise NullQuartet + end; + + (* Schwartz screening *) + (* + if (maxm > 8) then + ( + let schwartz_p = + let key = + Zkey.of_powers_twelve angMom_a angMom_b angMom_a angMom_b + in + match schwartz_p with + | None -> 1. + | Some schwartz_p -> Zmap.find schwartz_p key + in + if schwartz_p < cutoff then raise NullQuartet; + let schwartz_q = + let key = + Zkey.of_powers_twelve angMom_c angMom_d angMom_c angMom_d + in + match schwartz_q with + | None -> 1. + | Some schwartz_q -> Zmap.find schwartz_q key + in + if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; + ); + *) + + + let norm = norm_scales.(i) in + let coef_prod = coef_prod *. norm in + + let integral = + hvrr_two_e + angMom_a angMom_b angMom_c angMom_d + zero_m_array + expo_b expo_d + expo_inv_p expo_inv_q + center_ab center_cd center_pq + center_pa center_qc + map_1d map_2d + in + contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral + with NullQuartet -> () + ) + end + ) (Cspc.coefs_and_shell_pair_couples cspc) + ) (Aspc.contracted_shell_pair_couples atomic_shell_pair_couple); + + 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/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 7a33819..d8abcf1 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -822,15 +822,3 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q with NullQuartet -> empty - -(** Computes all the two-electron integrals of the contracted shell quartet *) -let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = - - let shell_p = Csp.make ~cutoff shell_a shell_b - and shell_q = Csp.make ~cutoff shell_c shell_d - in - match shell_p, shell_q with - | Some shell_p, Some shell_q -> - contracted_class_shell_pairs ~zero_m shell_p shell_q - | _ -> empty -