From b5b6010bbf2b1e8f4ad6db205a42d42e4583e9d5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 27 Mar 2018 19:26:32 +0200 Subject: [PATCH] Introduced ShellPairCouples in ERI --- Basis/ERI.ml | 53 +++++++++++++++++++++--------- Basis/ERI.mli | 4 +-- Basis/TwoElectronRR.ml | 73 ++++++------------------------------------ 3 files changed, 49 insertions(+), 81 deletions(-) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 480402e..675ab55 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -12,6 +12,7 @@ module Asp = AtomicShellPair module Bs = Basis module Cs = ContractedShell module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple module Fis = FourIdxStorage @@ -23,6 +24,7 @@ let set_phys = Fis.set_phys let to_file = Fis.to_file +let cutoff = integrals_cutoff (** (00|00)^m : Fundamental electron repulsion integral @@ -56,9 +58,12 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = -let class_of_contracted_shell_pairs shell_p shell_q = +let class_of_contracted_shell_pair_couple shell_pair_couple = + let shell_p = Cspc.shell_pair_p shell_pair_couple + and shell_q = Cspc.shell_pair_q shell_pair_couple + in if Array.length (Csp.shell_pairs shell_p) + (Array.length (Csp.shell_pairs shell_q)) < 4 then - TwoElectronRR.contracted_class_shell_pairs ~zero_m shell_p shell_q + TwoElectronRR.contracted_class_shell_pair_couple ~zero_m shell_pair_couple else TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m shell_p shell_q @@ -67,6 +72,19 @@ let class_of_contracted_shell_pairs shell_p shell_q = let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = + List.map (fun pair -> + match Cspc.make ~cutoff pair pair with + | Some cspc -> + let cls = class_of_contracted_shell_pair_couple cspc in + (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) + | None -> (pair, -1.) + ) shell_pairs + |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) + |> List.map fst + + +(* TODO +let filter_contracted_shell_pair_couples ?(cutoff=integrals_cutoff) shell_pair_couples = List.map (fun pair -> let cls = class_of_contracted_shell_pairs pair pair @@ -76,10 +94,10 @@ let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) |> List.map fst +*) - -let store_class ?(cutoff=integrals_cutoff) data cls shell_p shell_q = +let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = let to_powers x = let open Zkey in match to_powers x with @@ -87,6 +105,10 @@ let store_class ?(cutoff=integrals_cutoff) data cls shell_p shell_q = | _ -> assert false in + let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple + and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple + in + Array.iteri (fun i_c powers_i -> let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in let xi = to_powers powers_i in @@ -129,7 +151,7 @@ let of_basis basis = let shell_pairs = Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~cutoff:integrals_cutoff + |> filter_contracted_shell_pairs ~cutoff in Printf.printf "%d significant shell pairs computed in %f seconds\n" @@ -156,19 +178,18 @@ let of_basis basis = Cs.index (Csp.shell_a shell_p) then raise Exit in - let sq = - Csp.shell_pairs shell_q + let sq = Csp.shell_pairs shell_q in + let cspc = + if Array.length sp < Array.length sq then + Cspc.make ~cutoff shell_p shell_q + else + Cspc.make ~cutoff shell_q shell_p in - let f p q = - let cls = class_of_contracted_shell_pairs p q in - store_class ~cutoff:integrals_cutoff eri_array cls p q - in - - if Array.length sp > Array.length sq then - f shell_q shell_p - else - f shell_p shell_q + match cspc with + | Some cspc -> let cls = class_of_contracted_shell_pair_couple cspc in + store_class ~cutoff eri_array cspc cls + | None -> () ) shell_pairs with Exit -> () ) shell_pairs ; diff --git a/Basis/ERI.mli b/Basis/ERI.mli index 9969b73..4f13428 100644 --- a/Basis/ERI.mli +++ b/Basis/ERI.mli @@ -10,8 +10,8 @@ val filter_atomic_shell_pairs : ?cutoff:float -> AtomicShellPair.t list -> Atomi val filter_contracted_shell_pairs : ?cutoff:float -> ContractedShellPair.t list -> ContractedShellPair.t list (** Uses Schwartz screening on contracted shell pairs. *) -val class_of_contracted_shell_pairs : ContractedShellPair.t -> ContractedShellPair.t -> float Zmap.t -(** Computes all the ERI of the class built from th combination of two contracted shell pairs. *) +val class_of_contracted_shell_pair_couple : ContractedShellPairCouple.t -> float Zmap.t +(** Computes all the ERI of the class built from a couple of contracted shell pairs. *) val get_chem : t -> int -> int -> int -> int -> float diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 55bcecf..82499bb 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -14,7 +14,6 @@ module Ps = PrimitiveShell let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff -let empty = Zmap.create 0 exception NullQuartet @@ -278,11 +277,7 @@ let rec hvrr_two_e -let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - - match Cspc.make ~cutoff shell_p shell_q with - | None -> empty - | Some shell_pair_couple -> +let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t = let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in @@ -299,6 +294,10 @@ 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 shell_p = Cspc.shell_pair_p shell_pair_couple + and shell_q = Cspc.shell_pair_q shell_pair_couple + in + let center_ab = Csp.a_minus_b shell_p and center_cd = Csp.a_minus_b shell_q in @@ -357,32 +356,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q 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 @@ -415,11 +388,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q -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 contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : float Zmap.t = let maxm = Am.to_int (Aspc.ang_mom atomic_shell_pair_couple) in @@ -434,6 +403,10 @@ let contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p Aspc.monocentric atomic_shell_pair_couple in + let shell_p = Aspc.atomic_shell_pair_p atomic_shell_pair_couple + and shell_q = Aspc.atomic_shell_pair_q 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 @@ -495,32 +468,6 @@ let contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p 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