From 99a78616fdc5ef22df0d09b5f16042bc76910b52 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 28 Sep 2020 00:24:36 +0200 Subject: [PATCH] Added RS operators --- gaussian_integrals/lib/dependencies.ml | 4 - gaussian_integrals/lib/eri_lr.ml | 65 +++++++++++++ gaussian_integrals/lib/screened_eri.ml | 93 +++++++++++++++++++ .../lib/two_electron_integrals.ml | 2 +- .../lib/two_electron_integrals.mli | 4 +- gaussian_integrals/lib/zero_m_parameters.ml | 6 +- operators/lib/f12_operator.ml | 14 +-- operators/lib/gaussian_operator.ml | 6 +- operators/lib/operator.ml | 9 +- operators/lib/operator.mli | 15 ++- operators/lib/rs_operator.ml | 8 ++ 11 files changed, 198 insertions(+), 28 deletions(-) delete mode 100644 gaussian_integrals/lib/dependencies.ml create mode 100644 gaussian_integrals/lib/eri_lr.ml create mode 100644 gaussian_integrals/lib/screened_eri.ml create mode 100644 operators/lib/rs_operator.ml diff --git a/gaussian_integrals/lib/dependencies.ml b/gaussian_integrals/lib/dependencies.ml deleted file mode 100644 index 2b3c528..0000000 --- a/gaussian_integrals/lib/dependencies.ml +++ /dev/null @@ -1,4 +0,0 @@ -include Qcaml_common -include Qcaml_particles -include Qcaml_linear_algebra -include Qcaml_gaussian_basis diff --git a/gaussian_integrals/lib/eri_lr.ml b/gaussian_integrals/lib/eri_lr.ml new file mode 100644 index 0000000..101d556 --- /dev/null +++ b/gaussian_integrals/lib/eri_lr.ml @@ -0,0 +1,65 @@ +(** Long-range electron-electron repulsion integrals. + See Eq(52) in 10.1039/b605188j +*) + +open Qcaml_common +open Qcaml_gaussian_basis +open Qcaml_operators + +module Csp = Contracted_shell_pair +module Cspc = Contracted_shell_pair_couple + +module T = struct + + let name = "Long-range electron repulsion integrals" + + open Zero_m_parameters + + let zero_m z = + let mu_erf = + match z.operator with + | Some (Operator.Range_sep mu) -> mu + | _ -> assert false + in + let m = mu_erf *. mu_erf in + let expo_pq_inv = z.expo_p_inv +. z.expo_q_inv in + let fG_inv = expo_pq_inv +. 1. /. m in + let fG = 1. /. fG_inv in + assert (expo_pq_inv <> 0.); + let t = + if z.norm_pq_sq > Constants.integrals_cutoff then + z.norm_pq_sq *. fG + else 0. + in + let maxm = z.maxm in + let result = Util.boys_function ~maxm t in + let rec aux accu k = function + | 0 -> result.(k) <- result.(k) *. accu + | l -> + begin + result.(k) <- result.(k) *. accu; + let new_accu = -. accu *. fG in + (aux [@tailcall]) new_accu (k+1) (l-1) + end + in + let f = Constants.two_over_sq_pi *. (sqrt fG) in + aux f 0 maxm; + result + + let class_of_contracted_shell_pair_couple ?operator 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 + Two_electron_rr.contracted_class_shell_pair_couple + ?operator ~zero_m shell_pair_couple + else + Two_electron_rr_vectorized.contracted_class_shell_pairs + ?operator ~zero_m shell_p shell_q + +end + +module M = Two_electron_integrals.Make(T) +include M + diff --git a/gaussian_integrals/lib/screened_eri.ml b/gaussian_integrals/lib/screened_eri.ml new file mode 100644 index 0000000..2016e82 --- /dev/null +++ b/gaussian_integrals/lib/screened_eri.ml @@ -0,0 +1,93 @@ +(** Screened Electron-electron repulsion integrals (Yukawa potential). Required for F12/r12. *) + +open Qcaml_common +open Qcaml_gaussian_basis +open Qcaml_operators + +module Csp = Contracted_shell_pair +module Cspc = Contracted_shell_pair_couple + +module T = struct + + let name = "Screened electron repulsion integrals" + + open Zero_m_parameters + + let zero_m z = + let expo_pq_inv = z.expo_p_inv +. z.expo_q_inv in + assert (expo_pq_inv <> 0.); + + let expo_G_inv, coef_G = + let g = + match z.operator with + | Some (Operator.F12 f) -> f.F12_operator.gaussian + | _ -> assert false + in + g.Gaussian_operator.expo_g_inv, + g.Gaussian_operator.coef_g + in + + let expo_pq = 1. /. expo_pq_inv in + let maxm = z.maxm in + + let result = Array.init (maxm+1) (fun _ -> 0.) in + Util.array_range 0 (Array.length coef_G) + |> Array.iter (fun i -> + let fG_inv = expo_pq_inv +. expo_G_inv.(i) in + let fG = 1. /. fG_inv in + let t = + if z.norm_pq_sq > Constants.integrals_cutoff then + z.norm_pq_sq *. (expo_pq -. fG) + else 0. + in + let fm = Util.boys_function ~maxm t in + + let tmp_array = + let result = Array.init (maxm+1) (fun _ -> 1.) in + Util.array_range 1 maxm + |> Array.iter (fun k -> result.(k) <- result.(k-1) *. expo_pq *. expo_G_inv.(i)); + result + in + + let tmp_result = Array.init (maxm+1) (fun _ -> 0.) in + let rec aux accu m = function + | 0 -> tmp_result.(m) <- fm.(m) *. accu + | j -> + begin + let f = + Util.array_range 0 m + |> Array.fold_left (fun v k -> + v +. (Util.binom_float m k) *. tmp_array.(k) *. fm.(k)) 0. + in + tmp_result.(m) <- accu *. f; + let new_accu = -. accu *. expo_pq in + (aux [@tailcall]) new_accu (m+1) (j-1) + end + in + let f = + Constants.two_over_sq_pi *. (sqrt expo_pq) *. fG *. + expo_G_inv.(i) *. exp (-.fG *. z.norm_pq_sq) + in + aux f 0 maxm; + Array.iteri (fun k v -> + result.(k) <- result.(k) +. coef_G.(i) *. v + ) tmp_result) ; + result + + let class_of_contracted_shell_pair_couple ?operator 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 + Two_electron_rr.contracted_class_shell_pair_couple + ?operator ~zero_m shell_pair_couple + else + Two_electron_rr_vectorized.contracted_class_shell_pairs + ?operator ~zero_m shell_p shell_q + +end + +module M = Two_electron_integrals.Make(T) +include M + diff --git a/gaussian_integrals/lib/two_electron_integrals.ml b/gaussian_integrals/lib/two_electron_integrals.ml index 0ad9d44..2859a6a 100644 --- a/gaussian_integrals/lib/two_electron_integrals.ml +++ b/gaussian_integrals/lib/two_electron_integrals.ml @@ -19,7 +19,7 @@ module type Two_ei_structure = sig val name : string val class_of_contracted_shell_pair_couple : - ?operator:'a Operator.t -> Cspc.t -> float Zmap.t + ?operator:Operator.t -> Cspc.t -> float Zmap.t end diff --git a/gaussian_integrals/lib/two_electron_integrals.mli b/gaussian_integrals/lib/two_electron_integrals.mli index 60f2709..83ce1f9 100644 --- a/gaussian_integrals/lib/two_electron_integrals.mli +++ b/gaussian_integrals/lib/two_electron_integrals.mli @@ -16,7 +16,7 @@ val name : string (** Name of the kind of integrals, for printing purposes. *) val class_of_contracted_shell_pair_couple : - ?operator:'a Operator.t -> Contracted_shell_pair_couple.t -> float Zmap.t + ?operator:Operator.t -> Contracted_shell_pair_couple.t -> float Zmap.t (** Returns an integral class from a couple of contracted shells. The results is stored in a Zmap. *) @@ -28,7 +28,7 @@ module Make : functor (T : Two_ei_structure) -> sig include module type of Four_idx_storage - val of_basis : ?operator:'a Operator.t -> Basis.t -> t + val of_basis : ?operator:Operator.t -> Basis.t -> t (** Compute all ERI's for a given {!Basis.t}. *) end diff --git a/gaussian_integrals/lib/zero_m_parameters.ml b/gaussian_integrals/lib/zero_m_parameters.ml index f44f45b..996771d 100644 --- a/gaussian_integrals/lib/zero_m_parameters.ml +++ b/gaussian_integrals/lib/zero_m_parameters.ml @@ -1,7 +1,7 @@ open Qcaml_common open Qcaml_operators -type 'a t = +type t = { expo_p_inv : float ; expo_q_inv : float ; @@ -10,8 +10,8 @@ type 'a t = center_pq : Coordinate.t ; center_pa : Coordinate.t ; center_qc : Coordinate.t ; - zero_m_func : 'a t -> float array ; - operator : 'a Operator.t option; + zero_m_func : t -> float array ; + operator : Operator.t option; } let zero ?operator zero_m_func = diff --git a/operators/lib/f12_operator.ml b/operators/lib/f12_operator.ml index 90b1267..7f98c22 100644 --- a/operators/lib/f12_operator.ml +++ b/operators/lib/f12_operator.ml @@ -7,9 +7,9 @@ type t = } -let make_gaussian_corr_factor expo_s coef_g expo_sg = - let expo_sg = - Array.map (fun x -> x *. expo_s *. expo_s) expo_sg +let make_gaussian_corr_factor expo_s coef_g expo_sg = + let expo_sg = + Array.map (fun x -> x *. expo_s *. expo_s) expo_sg in let gaussian = Gaussian_operator.make coef_g expo_sg in { expo_s ; gaussian } @@ -60,10 +60,10 @@ let simple_gaussian_geminal' expo_s = |] and expo_sg = [| - 0.63172472556807146570889699432882480323314666748046875; - 26.3759196683467962429858744144439697265625; - 0.63172102793029016876147352377302013337612152099609375; - 7.08429025944207335641067402320913970470428466796875; + 0.63172472556807146570889699432882480323314666748046875; + 26.3759196683467962429858744144439697265625; + 0.63172102793029016876147352377302013337612152099609375; + 7.08429025944207335641067402320913970470428466796875; 42.4442841447001910637482069432735443115234375; 391.44036073596890901171718724071979522705078125 ; |] diff --git a/operators/lib/gaussian_operator.ml b/operators/lib/gaussian_operator.ml index 0f47cd5..64bfe08 100644 --- a/operators/lib/gaussian_operator.ml +++ b/operators/lib/gaussian_operator.ml @@ -8,14 +8,14 @@ type t = } let make coef_g expo_g = - let expo_g_inv = - Array.map (fun x -> 1. /. x ) expo_g + let expo_g_inv = + Array.map (fun x -> 1. /. x ) expo_g in { coef_g ; expo_g ; expo_g_inv } let one_over_r = - + let coef_g = [| 841.88478132 ; 70.590185207 ; 18.3616020768 ; 7.2608642093 ; 3.57483416444 ; 2.01376031082 ; 1.24216542801 ; 0.81754348620 ; diff --git a/operators/lib/operator.ml b/operators/lib/operator.ml index 86f04e3..5c6c520 100644 --- a/operators/lib/operator.ml +++ b/operators/lib/operator.ml @@ -1,6 +1,9 @@ -type 'a t = -| F12 of F12_operator.t -| Gaussian of Gaussian_operator.t +type t = + | F12 of F12_operator.t + | Gaussian of Gaussian_operator.t + | Range_sep of Rs_operator.t let of_f12 f = F12 f let of_gaussian g = Gaussian g +let of_range_separation mu = Range_sep mu + diff --git a/operators/lib/operator.mli b/operators/lib/operator.mli index 82af419..c70f950 100644 --- a/operators/lib/operator.mli +++ b/operators/lib/operator.mli @@ -1,9 +1,14 @@ -type 'a t = -| F12 of F12_operator.t -| Gaussian of Gaussian_operator.t +type t = + | F12 of F12_operator.t + | Gaussian of Gaussian_operator.t + | Range_sep of Rs_operator.t -val of_f12 : F12_operator.t -> F12_operator.t t +val of_f12 : F12_operator.t -> t (** Creates a F12 operator *) -val of_gaussian: Gaussian_operator.t -> Gaussian_operator.t t +val of_gaussian: Gaussian_operator.t -> t (** Creates a Gaussian operator *) + +val of_range_separation: Rs_operator.t -> t +(** Creates a range-separated operator *) + diff --git a/operators/lib/rs_operator.ml b/operators/lib/rs_operator.ml new file mode 100644 index 0000000..507312f --- /dev/null +++ b/operators/lib/rs_operator.ml @@ -0,0 +1,8 @@ +(** Type for range-separation parameter *) + +type t = float + + + + +