From bb0b52f4ea1d715f734521dce15ae35d966a1934 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 13 Mar 2019 22:02:08 +0100 Subject: [PATCH] F12 integrals --- Basis/AOBasis.ml | 8 +- Basis/AOBasis.mli | 3 + Basis/ERI.ml | 2 +- Basis/F12.ml | 4 + Basis/F12RR.ml | 290 +++++++++++++++++ ...ml => TwoElectronIntegralsNonSeparable.ml} | 57 ++-- ...i => TwoElectronIntegralsNonSeparable.mli} | 5 - Basis/TwoElectronIntegralsSeparable.ml | 307 ++++++++++++++++++ Basis/TwoElectronRR.ml | 27 +- Basis/TwoElectronRRVectorized.ml | 48 ++- Basis/Zero_m_parameters.ml | 10 +- run_integrals.ml | 5 +- 12 files changed, 700 insertions(+), 66 deletions(-) create mode 100644 Basis/F12.ml create mode 100644 Basis/F12RR.ml rename Basis/{TwoElectronIntegrals.ml => TwoElectronIntegralsNonSeparable.ml} (89%) rename Basis/{TwoElectronIntegrals.mli => TwoElectronIntegralsNonSeparable.mli} (86%) create mode 100644 Basis/TwoElectronIntegralsSeparable.ml diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index 5de5b24..4d0fb89 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -9,6 +9,7 @@ type t = eN_ints : NucInt.t lazy_t; kin_ints : KinInt.t lazy_t; ee_ints : ERI.t lazy_t; + f12_ints : F12.t lazy_t; cartesian : bool; } @@ -18,6 +19,7 @@ let ortho t = Lazy.force t.ortho let eN_ints t = Lazy.force t.eN_ints let kin_ints t = Lazy.force t.kin_ints let ee_ints t = Lazy.force t.ee_ints +let f12_ints t = Lazy.force t.f12_ints let cartesian t = t.cartesian @@ -44,7 +46,11 @@ let make ~cartesian ~basis nuclei = ERI.of_basis basis ) in - { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; + let f12_ints = lazy ( + F12.of_basis basis + ) in + + { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ; cartesian ; } diff --git a/Basis/AOBasis.mli b/Basis/AOBasis.mli index ea9d66e..8694136 100644 --- a/Basis/AOBasis.mli +++ b/Basis/AOBasis.mli @@ -19,6 +19,9 @@ val eN_ints : t -> NucInt.t val ee_ints : t -> ERI.t (** Electron-electron potential integrals *) +val f12_ints : t -> F12.t +(** Electron-electron potential integrals *) + val kin_ints : t -> KinInt.t (** Kinetic energy integrals *) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 1237a38..697abb9 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -37,6 +37,6 @@ module Zm = struct end -module M = TwoElectronIntegrals.Make(Zm) +module M = TwoElectronIntegralsNonSeparable.Make(Zm) include M diff --git a/Basis/F12.ml b/Basis/F12.ml new file mode 100644 index 0000000..29d832f --- /dev/null +++ b/Basis/F12.ml @@ -0,0 +1,4 @@ +(** Two-electron integrals over Slater geminals via a fit using Gaussian geminals. +*) + +include TwoElectronIntegralsSeparable diff --git a/Basis/F12RR.ml b/Basis/F12RR.ml new file mode 100644 index 0000000..3f61c06 --- /dev/null +++ b/Basis/F12RR.ml @@ -0,0 +1,290 @@ +open Util +open Constants + +module Am = AngularMomentum +module Asp = AtomicShellPair +module Aspc = AtomicShellPairCouple +module Co = Coordinate +module Cs = ContractedShell +module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple +module Po = Powers +module Psp = PrimitiveShellPair +module Pspc = PrimitiveShellPairCouple +module Ps = PrimitiveShell +module Zp = Zero_m_parameters + +let cutoff = Constants.integrals_cutoff +let cutoff2 = cutoff *. cutoff + +exception NullQuartet + +type four_idx_intermediates = + { + zero : float array array; + gp : float array; + gg : float array; + gq : float array; + coef_g : float array ; + center_ra : Co.t array ; + center_rc : Co.t array ; + center_ab : Co.t ; + center_cd : Co.t ; + } + +(** Horizontal and Vertical Recurrence Relations (HVRR) *) +let rec hvrr angMom_a angMom_b angMom_c angMom_d + abcd map_x map_y map_z = + + let gp = abcd.gp in + let gg = abcd.gg in + let gq = abcd.gq in + + let coef_g = abcd.coef_g in + + let run xyz map = + + let zero = + match xyz with + | Co.X -> abcd.zero.(0) + | Co.Y -> abcd.zero.(1) + | Co.Z -> abcd.zero.(2) + in + let angMom_a = Po.get xyz angMom_a in + let angMom_b = Po.get xyz angMom_b in + let angMom_c = Po.get xyz angMom_c in + let angMom_d = Po.get xyz angMom_d in + let center_ab = Co.get xyz abcd.center_ab in + let center_cd = Co.get xyz abcd.center_cd in + let center_ra = Array.map (fun x -> Co.get xyz x) abcd.center_ra in + let center_rc = Array.map (fun x -> Co.get xyz x) abcd.center_rc in + + let rec vrr angMom_a angMom_c = + match angMom_a, angMom_c with + | 0, -1 + | -1, 0 -> assert false + | 0, 0 -> zero + | 1, 0 -> + let v1 = zero in + Array.mapi (fun i v1i -> center_ra.(i) *. v1i ) v1 + | 0, 1 -> + let v1 = zero in + Array.mapi (fun i v1i -> center_rc.(i) *. v1i ) v1 + | 1, 1 -> + let v1 = vrr 1 0 in + let v2 = zero in + Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. + gg.(i) *. v2.(i) ) v1 + | 2, 0 -> + let v1 = vrr 1 0 in + let v2 = zero in + Array.mapi (fun i v1i -> center_ra.(i) *. v1i +. gp.(i) *. v2.(i)) v1 + | _, 0 -> + let v1 = + vrr (angMom_a-1) 0 + in + let a = float_of_int (angMom_a-1) in + let v2 = + vrr (angMom_a-2) 0 + in + Array.mapi (fun i v1i -> center_ra.(i) *. v1i +. a *. gp.(i) *. v2.(i)) v1 + | _, 1 -> + let v1 = + vrr angMom_a 0 + in + let a = float_of_int angMom_a in + let v2 = + vrr (angMom_a-1) 0 + in + Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. + a *. gg.(i) *. v2.(i) ) v1 + | 0, _ -> + let v1 = + vrr 0 (angMom_c-1) + in + let b = float_of_int (angMom_c-1) in + let v3 = + vrr 0 (angMom_c-2) + in + Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. + b *. gq.(i) *. v3.(i)) v1 + | _ -> + let v1 = + vrr angMom_a (angMom_c-1) + in + let a = float_of_int angMom_a in + let b = float_of_int (angMom_c-1) in + let v2 = + vrr (angMom_a-1) (angMom_c-1) + in + let v3 = + vrr angMom_a (angMom_c-2) + in + Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. + a *. gg.(i) *. v2.(i) +. b *. gq.(i) *. v3.(i)) v1 + in + + let rec hrr angMom_a angMom_b angMom_c angMom_d = + let key = Zkey.of_int_four angMom_a angMom_b angMom_c angMom_d in + + try Zmap.find map key with + | Not_found -> + let result = + match angMom_b, angMom_d with + | -1, 0 + | 0, -1 -> assert false + | 0, 0 -> + vrr angMom_a angMom_c + | _, 0 -> + let h1 = + hrr (angMom_a+1) (angMom_b-1) angMom_c angMom_d + in + if center_ab = 0. then + h1 + else + let h2 = + hrr angMom_a (angMom_b-1) angMom_c angMom_d + in + Array.mapi (fun i h1i -> h1i +. center_ab *. h2.(i)) h1 + | _ -> + let h1 = + hrr angMom_a angMom_b (angMom_c+1) (angMom_d-1) + in + if center_cd = 0. then + h1 + else + let h2 = + hrr angMom_a angMom_b angMom_c (angMom_d-1) + in + Array.mapi (fun i h1i -> h1i +. center_cd *. h2.(i)) h1 + in (Zmap.add map key result; result) + + in + hrr angMom_a angMom_b angMom_c angMom_d + in + let x, y, z = + (run Co.X map_x), (run Co.Y map_y), (run Co.Z map_z) + in + let rec aux accu = function + | 0 -> accu +. coef_g.(0) *. x.(0) *. y.(0) *. z.(0) + | i -> aux (accu +. coef_g.(i) *. x.(i) *. y.(i) *. z.(i)) (i-1) + in + aux 0. (Array.length x - 1) + + +let contracted_class_shell_pair_couple expo_g_inv coef_g shell_pair_couple : float Zmap.t = + + (* Pre-computation of integral class indices *) + let class_indices = Cspc.zkey_array shell_pair_couple in + + let contracted_class = + Array.make (Array.length class_indices) 0.; + in + + (* 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 + + let norm_scales = Cspc.norm_scales shell_pair_couple in + + 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_p_inv = Psp.exponent_inv sp_ab in + let expo_q_inv = Psp.exponent_inv sp_cd in + let expo_pgq = Array.map (fun e -> + 1. /. (expo_p_inv +. expo_q_inv +. e) + ) expo_g_inv + in + + let fp = Array.map (fun e -> expo_p_inv *. e) expo_pgq in + let fq = Array.map (fun e -> expo_q_inv *. e) expo_pgq in + + let gp = + let x = 0.5 *. expo_p_inv in + Array.map (fun f -> (1. -. f) *. x) fp + in + let gq = + let x = 0.5 *. expo_q_inv in + Array.map (fun f -> (1. -. f) *. x) fq + in + let gg = + let x = 0.5 *. expo_q_inv in + Array.map (fun f -> f *. x) fp + in + + let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in + let center_qc = Psp.center_minus_a sp_cd in + let center_pa = Psp.center_minus_a sp_ab in + + let center_ra = Array.map (fun f -> Co.(center_pa |- (f |. center_pq))) fp in + let center_rc = Array.map (fun f -> Co.(center_qc |+ (f |. center_pq))) fq in + + (* zero_m is defined here *) + let zero = Array.map (fun xyz -> + let x = Co.get xyz center_pq in + Array.mapi (fun i ei -> + let fg = expo_g_inv.(i) *. ei in + (sqrt fg) *. exp (-. x *. x *. ei )) expo_pgq + ) Co.[| X ; Y ; Z |] + in + begin + match Cspc.ang_mom shell_pair_couple with + (* + | Am.S -> + let integral = + zero.(0) *. zero.(1) *. zero.(2) + in + contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral + *) + | _ -> + let map_x, map_y, map_z = + Zmap.create (Array.length class_indices), + Zmap.create (Array.length class_indices), + Zmap.create (Array.length class_indices) + 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 + let norm = norm_scales.(i) in + let coef_prod = coef_prod *. norm in + let abcd = { + zero ; gp ; gg ; gq ; + center_ab ; center_cd ; + center_ra ; center_rc ; + coef_g ; + } in + let integral = + hvrr + angMom_a angMom_b angMom_c angMom_d + abcd map_x map_y map_z + in + contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral + ) + end + ) (Cspc.coefs_and_shell_pair_couples 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/TwoElectronIntegrals.ml b/Basis/TwoElectronIntegralsNonSeparable.ml similarity index 89% rename from Basis/TwoElectronIntegrals.ml rename to Basis/TwoElectronIntegralsNonSeparable.ml index 786e1b8..0b58a80 100644 --- a/Basis/TwoElectronIntegrals.ml +++ b/Basis/TwoElectronIntegralsNonSeparable.ml @@ -1,19 +1,16 @@ -(** Two electron integral functor, parameterized by the zero_m function. +(** Two electron integral functor for operators that are not separable among %{ $(x,y,z)$ %}. + It is parameterized by the [zero_m] function. *) open Constants let cutoff = integrals_cutoff -module Am = AngularMomentum -module As = AtomicShell -module Asp = AtomicShellPair module Bs = Basis module Cs = ContractedShell module Csp = ContractedShellPair module Cspc = ContractedShellPairCouple module Fis = FourIdxStorage - module type Zero_mType = sig val name : string @@ -39,34 +36,32 @@ module Make(Zero_m : Zero_mType) = struct TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m 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. ) - (* TODO \sum_k |coef_k * integral_k| *) - | None -> (pair, -1.) + 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. ) + (* TODO \sum_k |coef_k * integral_k| *) + | 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 - in - (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - ) 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 + in + (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) + ) shell_pairs + |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) + |> List.map fst - *) +*) let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = @@ -114,9 +109,9 @@ module Make(Zero_m : Zero_mType) = struct let eri_array = Fis.create ~size:n `Dense -(* + (* Fis.create ~size:n `Sparse -*) + *) in let t0 = Unix.gettimeofday () in @@ -308,7 +303,3 @@ end - - - - diff --git a/Basis/TwoElectronIntegrals.mli b/Basis/TwoElectronIntegralsNonSeparable.mli similarity index 86% rename from Basis/TwoElectronIntegrals.mli rename to Basis/TwoElectronIntegralsNonSeparable.mli index c95761b..43ecb24 100644 --- a/Basis/TwoElectronIntegrals.mli +++ b/Basis/TwoElectronIntegralsNonSeparable.mli @@ -30,11 +30,6 @@ module Make : functor (Zero_m : Zero_mType) -> sig include module type of FourIdxStorage - 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 filter_contracted_shell_pairs : ?cutoff:float -> ContractedShellPair.t list -> ContractedShellPair.t list diff --git a/Basis/TwoElectronIntegralsSeparable.ml b/Basis/TwoElectronIntegralsSeparable.ml new file mode 100644 index 0000000..ebc54cb --- /dev/null +++ b/Basis/TwoElectronIntegralsSeparable.ml @@ -0,0 +1,307 @@ +(** Two electron integral functor for operators that are separable among %{ $(x,y,z)$ %}. + It is parameterized by the [zero_m] function. +*) + +open Constants +let cutoff = integrals_cutoff + +module Bs = Basis +module Cs = ContractedShell +module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple +module Fis = FourIdxStorage + + +include FourIdxStorage + +(** Exponent of the geminal *) +let expo_s = 1.0 + +(** Coefficients and exponents of the Gaussian fit *) +let coef_g = + [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] + +let expo_sg_inv = + Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) + [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] + + +let class_of_contracted_shell_pair_couple shell_pair_couple = + F12RR.contracted_class_shell_pair_couple + expo_sg_inv coef_g shell_pair_couple + + + + +module Zero_m = struct + let name = "F12" +end + + + +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. ) + (* TODO \sum_k |coef_k * integral_k| *) + | 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 + in + (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) + ) shell_pairs + |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) + |> List.map fst + +*) + + +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 + | Three x -> x + | _ -> 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 + Array.iteri (fun j_c powers_j -> + let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in + let xj = to_powers powers_j in + Array.iteri (fun k_c powers_k -> + let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in + let xk = to_powers powers_k in + Array.iteri (fun l_c powers_l -> + let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in + let xl = to_powers powers_l in + let key = Zkey.of_powers_twelve xi xj xk xl in + let value = Zmap.find cls key in + set_chem data i_c j_c k_c l_c value + ) (Cs.zkey_array (Csp.shell_b shell_q)) + ) (Cs.zkey_array (Csp.shell_a shell_q)) + ) (Cs.zkey_array (Csp.shell_b shell_p)) + ) (Cs.zkey_array (Csp.shell_a shell_p)) + + + + + + +let of_basis_serial basis = + + let n = Bs.size basis + and shell = Bs.contracted_shells basis + in + + let eri_array = + Fis.create ~size:n `Dense +(* + Fis.create ~size:n `Sparse +*) + in + + let t0 = Unix.gettimeofday () in + + let shell_pairs = + Csp.of_contracted_shell_array shell + |> filter_contracted_shell_pairs ~cutoff + in + + Printf.printf "%d significant shell pairs computed in %f seconds\n" + (List.length shell_pairs) (Unix.gettimeofday () -. t0); + + + let t0 = Unix.gettimeofday () in + let ishell = ref 0 in + + List.iter (fun shell_p -> + let () = + if (Cs.index (Csp.shell_a shell_p) > !ishell) then + (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) + in + + let sp = + Csp.shell_pairs shell_p + in + + try + List.iter (fun shell_q -> + let () = + if Cs.index (Csp.shell_a shell_q) > + Cs.index (Csp.shell_a shell_p) then + raise Exit + in + 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 + + 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 ; + Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); + eri_array + + + + + + + +(* Parallel functions *) + + + +let of_basis_parallel basis = + + let n = Bs.size basis + and shell = Bs.contracted_shells basis + in + + let store_class_parallel + ?(cutoff=integrals_cutoff) contracted_shell_pair_couple cls = + let to_powers x = + let open Zkey in + match to_powers x with + | Three x -> x + | _ -> 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 + + let result = ref [] 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 + Array.iteri (fun j_c powers_j -> + let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in + let xj = to_powers powers_j in + Array.iteri (fun k_c powers_k -> + let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in + let xk = to_powers powers_k in + Array.iteri (fun l_c powers_l -> + let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in + let xl = to_powers powers_l in + let key = Zkey.of_powers_twelve xi xj xk xl in + let value = Zmap.find cls key in + result := (i_c, j_c, k_c, l_c, value) :: !result + ) (Cs.zkey_array (Csp.shell_b shell_q)) + ) (Cs.zkey_array (Csp.shell_a shell_q)) + ) (Cs.zkey_array (Csp.shell_b shell_p)) + ) (Cs.zkey_array (Csp.shell_a shell_p)); + !result + in + + + + let t0 = Unix.gettimeofday () in + + let shell_pairs = + Csp.of_contracted_shell_array shell + |> filter_contracted_shell_pairs ~cutoff + in + + if Parallel.master then + Printf.printf "%d significant shell pairs computed in %f seconds\n" + (List.length shell_pairs) (Unix.gettimeofday () -. t0); + + + let t0 = Unix.gettimeofday () in + let ishell = ref max_int in + + let input_stream = Stream.of_list (List.rev shell_pairs) in + + let f shell_p = + let () = + if Parallel.rank < 2 && Cs.index (Csp.shell_a shell_p) < !ishell then + (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) + in + + let sp = + Csp.shell_pairs shell_p + in + + let result = ref [] in + try + List.iter (fun shell_q -> + let () = + if Cs.index (Csp.shell_a shell_q) > + Cs.index (Csp.shell_a shell_p) then + raise Exit + in + 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 + + match cspc with + | Some cspc -> + let cls = + class_of_contracted_shell_pair_couple cspc + in + result := (store_class_parallel ~cutoff cspc cls) :: !result; + | None -> () + ) shell_pairs; + raise Exit + with Exit -> List.concat !result |> Array.of_list + in + + let eri_array = + if Parallel.master then + Fis.create ~size:n `Dense + else + Fis.create ~size:0 `Dense + in + Farm.run ~ordered:false ~f input_stream + |> Stream.iter (fun l -> + Array.iter (fun (i_c,j_c,k_c,l_c,value) -> + set_chem eri_array i_c j_c k_c l_c value) l); + + if Parallel.master then + Printf.printf + "Computed %s Integrals in parallel in %f seconds\n%!" Zero_m.name (Unix.gettimeofday () -. t0); + Parallel.broadcast (lazy eri_array) + + + +let of_basis = + match Parallel.size with + | 1 -> of_basis_serial + | _ -> of_basis_parallel + + + + + diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 387e041..2a51148 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -337,13 +337,16 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t in let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in + let center_pa = Psp.center_minus_a sp_ab in + let center_qc = Psp.center_minus_a sp_cd in let norm_pq_sq = Co.dot center_pq center_pq in let expo_q_inv = Psp.exponent_inv sp_cd in - + let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in let zero_m_array = - Zp.{ - maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq - } |> zero_m + zero_m Zp.{ + maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; + center_pq ; center_pa ; center_qc ; normalization ; + } in begin @@ -354,13 +357,10 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t | _ -> 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 @@ -450,13 +450,17 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : in let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in + let center_qc = Psp.center_minus_a sp_cd in + let center_pa = Psp.center_minus_a sp_ab in let norm_pq_sq = Co.dot center_pq center_pq in let expo_q_inv = Psp.exponent_inv sp_cd in + let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in let zero_m_array = - Zp.{ - maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq - } |> zero_m + zero_m Zp.{ + maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; + center_pq ; center_pa ; center_qc ; normalization ; + } in begin @@ -467,13 +471,10 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : | _ -> 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 diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index a5856c5..d911356 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -643,15 +643,23 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let center_pq = Co.(Psp.center sp.(i-1) |- Psp.center sq.(j-1)) + and center_pa = + Co.(Psp.center sp.(i-1) |- Cs.center shell_a) + and center_qc = + Co.(Psp.center sq.(i-1) |- Cs.center shell_c) in let norm_pq_sq = Co.dot center_pq center_pq in + let normalization = Psp.normalization sp.(i-1) *. + Psp.normalization sq.(i-1) + in let zero_m_array = - Zp.{ - maxm=0 ; expo_p_inv ; expo_q_inv ; norm_pq_sq - } |> zero_m + zero_m Zp.{ + maxm=0 ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; + center_pq ; center_pa ; center_qc ; normalization ; + } in zero_m_array.(0) with NullQuartet -> 0. @@ -747,7 +755,20 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.init np (fun _ -> Array.make nq 0. ) ) in let empty = Array.make (maxm+1) 0. in + let center_qc_tmp = Array.init nq (fun cd -> + Bohr.make { Point. + x = (center_qc Co.X).(cd) ; + y = (center_qc Co.Y).(cd) ; + z = (center_qc Co.Z).(cd) ; + }) + in Array.iteri (fun ab shell_ab -> + let center_pa = Bohr.make { Point. + x = (center_pa Co.X).(ab) ; + y = (center_pa Co.Y).(ab) ; + z = (center_pa Co.Z).(ab) ; + } + in let zero_m_array_tmp = Array.mapi (fun cd shell_cd -> if (abs_float coef.(ab).(cd) < cutoff) then @@ -756,17 +777,22 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let expo_p_inv, expo_q_inv = expo_p_inv.(ab), expo_q_inv.(cd) in + let x = (center_pq X).(ab).(cd) + and y = (center_pq Y).(ab).(cd) + and z = (center_pq Z).(ab).(cd) + in let norm_pq_sq = - let x = (center_pq X).(ab).(cd) - and y = (center_pq Y).(ab).(cd) - and z = (center_pq Z).(ab).(cd) - in x *. x +. y *. y +. z *. z in - - Zp.{ - maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq - } |> zero_m + let normalization = Psp.normalization shell_ab *. + Psp.normalization shell_cd + in + zero_m Zp.{ + maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; + center_pq = Bohr.make Point.{x ; y ; z} ; + center_pa ; center_qc = center_qc_tmp.(cd) ; + normalization ; + } ) sq in (* Transpose result *) diff --git a/Basis/Zero_m_parameters.ml b/Basis/Zero_m_parameters.ml index efb6131..0d8ebd8 100644 --- a/Basis/Zero_m_parameters.ml +++ b/Basis/Zero_m_parameters.ml @@ -3,7 +3,11 @@ type t = expo_p_inv : float ; expo_q_inv : float ; norm_pq_sq : float ; - maxm : int; + normalization : float ; + maxm : int ; + center_pq : Coordinate.t ; + center_pa : Coordinate.t ; + center_qc : Coordinate.t ; } let zero = @@ -12,4 +16,8 @@ let zero = expo_p_inv = 0.; expo_q_inv = 0.; norm_pq_sq = 0.; + normalization = 0.; + center_pq = Coordinate.zero ; + center_pa = Coordinate.zero ; + center_qc = Coordinate.zero ; } diff --git a/run_integrals.ml b/run_integrals.ml index 738e19f..cf8f827 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -41,10 +41,13 @@ let run ~out = let eN_ints = AOBasis.eN_ints ao_basis in let kin_ints = AOBasis.kin_ints ao_basis in let ee_ints = AOBasis.ee_ints ao_basis in + let f12_ints = AOBasis.f12_ints ao_basis in Overlap.to_file ~filename:(out_file^".overlap") overlap; NucInt.to_file ~filename:(out_file^".nuc") eN_ints; KinInt.to_file ~filename:(out_file^".kin") kin_ints; - ERI.to_file ~filename:(out_file^".eri") ee_ints + ERI.to_file ~filename:(out_file^".eri") ee_ints; + F12.to_file ~filename:(out_file^".f12") f12_ints; + () let () =