From 8cb1d58a279a4191cf5f26ff5530a7a540445054 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 27 Sep 2020 23:55:42 +0200 Subject: [PATCH] Added operators --- gaussian_integrals/lib/dune | 1 + gaussian_integrals/lib/eri.ml | 56 ++++ gaussian_integrals/lib/f12.ml | 32 ++ gaussian_integrals/lib/f12_rr.ml | 287 ++++++++++++++++++ gaussian_integrals/lib/kinetic.ml | 56 ++-- gaussian_integrals/lib/overlap.ml | 9 +- .../lib/two_electron_integrals.ml | 36 +-- .../lib/two_electron_integrals.mli | 10 +- gaussian_integrals/lib/two_electron_rr.ml | 8 +- .../lib/two_electron_rr_vectorized.ml | 6 +- gaussian_integrals/lib/zero_m_parameters.ml | 12 +- operators/lib/dune | 10 + operators/lib/f12_operator.ml | 76 +++++ operators/lib/gaussian_operator.ml | 37 +++ operators/lib/operator.ml | 6 + operators/lib/operator.mli | 9 + 16 files changed, 567 insertions(+), 84 deletions(-) create mode 100644 gaussian_integrals/lib/eri.ml create mode 100644 gaussian_integrals/lib/f12.ml create mode 100644 gaussian_integrals/lib/f12_rr.ml create mode 100644 operators/lib/dune create mode 100644 operators/lib/f12_operator.ml create mode 100644 operators/lib/gaussian_operator.ml create mode 100644 operators/lib/operator.ml create mode 100644 operators/lib/operator.mli diff --git a/gaussian_integrals/lib/dune b/gaussian_integrals/lib/dune index a467a86..d5c4472 100644 --- a/gaussian_integrals/lib/dune +++ b/gaussian_integrals/lib/dune @@ -8,6 +8,7 @@ qcaml.common qcaml.linear_algebra qcaml.gaussian_basis + qcaml.operators ) (modules_without_implementation matrix_on_basis) (synopsis "Integrals on the Gaussian basis sets.")) diff --git a/gaussian_integrals/lib/eri.ml b/gaussian_integrals/lib/eri.ml new file mode 100644 index 0000000..b8c7640 --- /dev/null +++ b/gaussian_integrals/lib/eri.ml @@ -0,0 +1,56 @@ +(** Electron-electron repulsion integrals *) + +open Qcaml_common +open Qcaml_gaussian_basis + +module Csp = Contracted_shell_pair +module Cspc = Contracted_shell_pair_couple + +module T = struct + + let name = "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_pq = 1. /. expo_pq_inv in + let t = + if z.norm_pq_sq > Constants.integrals_cutoff then + z.norm_pq_sq *. expo_pq + 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 *. expo_pq in + (aux [@tailcall]) new_accu (k+1) (l-1) + end + in + let f = Constants.two_over_sq_pi *. (sqrt expo_pq) in + aux f 0 maxm; + result + + let class_of_contracted_shell_pair_couple ?operator shell_pair_couple = + assert (operator = None); + 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 + ~zero_m shell_pair_couple + else + Two_electron_rr_vectorized.contracted_class_shell_pairs + ~zero_m shell_p shell_q + +end + +module M = Two_electron_integrals.Make(T) +include M + diff --git a/gaussian_integrals/lib/f12.ml b/gaussian_integrals/lib/f12.ml new file mode 100644 index 0000000..1f9b69e --- /dev/null +++ b/gaussian_integrals/lib/f12.ml @@ -0,0 +1,32 @@ +(** Two electron integral functor for operators that are separable among %{ $(x,y,z)$ %}. + It is parameterized by the [zero_m] function. +*) + +open Qcaml_common +open Qcaml_operators + +open Constants +let cutoff = integrals_cutoff + +module T = struct + + let name = "F12" + + let class_of_contracted_shell_pair_couple ?operator shell_pair_couple = + let g = + match operator with + | Some (Operator.F12 f) -> f.gaussian + | _ -> failwith "Expected F12 operator" + in + F12_rr.contracted_class_shell_pair_couple + g.Gaussian_operator.expo_g_inv + g.Gaussian_operator.coef_g + shell_pair_couple + +end + +module M = Two_electron_integrals.Make(T) +include M + + + diff --git a/gaussian_integrals/lib/f12_rr.ml b/gaussian_integrals/lib/f12_rr.ml new file mode 100644 index 0000000..5af2337 --- /dev/null +++ b/gaussian_integrals/lib/f12_rr.ml @@ -0,0 +1,287 @@ +open Qcaml_common +open Qcaml_gaussian_basis + +module Am = Angular_momentum +module Asp = Atomic_shell_pair +module Aspc = Atomic_shell_pair_couple +module Co = Coordinate +module Cs = Contracted_shell +module Csp = Contracted_shell_pair +module Cspc = Contracted_shell_pair_couple +module Po = Powers +module Psp = Primitive_shell_pair +module Pspc = Primitive_shell_pair_couple +module Ps = Primitive_shell + +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 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 = Util.float_of_int_fast (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 = Util.float_of_int_fast 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 = Util.float_of_int_fast (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 = Util.float_of_int_fast angMom_a in + let b = Util.float_of_int_fast (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 [@tailcall]) (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/gaussian_integrals/lib/kinetic.ml b/gaussian_integrals/lib/kinetic.ml index 29c7ea9..5e28e5a 100644 --- a/gaussian_integrals/lib/kinetic.ml +++ b/gaussian_integrals/lib/kinetic.ml @@ -10,8 +10,8 @@ module Co = Coordinate module Cs = Contracted_shell module Csp = Contracted_shell_pair module Po = Powers -module Ps = Primitive_shell module Psp = Primitive_shell_pair +module Ps = Primitive_shell type t = Matrix.t external matrix : t -> Matrix.t = "%identity" @@ -19,52 +19,54 @@ external of_matrix : Matrix.t -> t = "%identity" let cutoff = integrals_cutoff -let to_powers x = +let to_powers x = let open Zkey in match to_powers x with | Six x -> x | _ -> assert false (** Computes all the kinetic integrals of the contracted shell pair *) -let contracted_class shell_a shell_b : float Zmap.t = +let contracted_class shell_a shell_b : float Zmap.t = match Csp.make shell_a shell_b with + | None -> Zmap.create 0 | Some shell_p -> begin - (* Pre-computation of integral class indices *) - let class_indices = Csp.zkey_array shell_p in - let contracted_class = + (* Pre-computation of integral class indices *) + let class_indices = Csp.zkey_array shell_p 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 sp = Csp.shell_pairs shell_p in - let a_minus_b = + let a_minus_b = Csp.a_minus_b shell_p in let norm_coef_scales = Csp.norm_scales shell_p in + (* Compute all integrals in the shell for each pair of significant shell pairs *) + + let sp = Csp.shell_pairs shell_p in for ab=0 to (Array.length sp - 1) do let coef_prod = - (Csp.coefficients shell_p).(ab) + (Csp.coefficients shell_p).(ab) in (* Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-4*.cutoff then begin - let center_pa = + let center_pa = Psp.center_minus_a sp.(ab) in - let expo_inv = + let expo_inv = (Csp.exponents_inv shell_p).(ab) in - let expo_a = + let expo_a = Ps.exponent (Psp.shell_a sp.(ab)) - and expo_b = + and expo_b = Ps.exponent (Psp.shell_b sp.(ab)) in @@ -76,26 +78,26 @@ let contracted_class shell_a shell_b : float Zmap.t = in Array.iteri (fun i key -> let (angMomA,angMomB) = to_powers key in - let ov a b k = + let ov a b k = let xyz = xyz_of_int k in Overlap_primitives.hvrr (a, b) expo_inv (Co.get xyz a_minus_b, - Co.get xyz center_pa) + Co.get xyz center_pa) in - let f k = - let xyz = xyz_of_int k in + let f k = + let xyz = xyz_of_int k in ov (Po.get xyz angMomA) (Po.get xyz angMomB) k and g k = - let xyz = xyz_of_int k in + let xyz = xyz_of_int k in let s1 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB - 1) k and s2 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB - 1) k and s3 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB + 1) k and s4 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB + 1) k and a = float_of_int_fast (Po.get xyz angMomA) - and b = float_of_int_fast (Po.get xyz angMomB) + and b = float_of_int_fast (Po.get xyz angMomB) in - 0.5 *. a *. b *. s1 -. expo_a *. b *. s2 -. expo_b *. a *. s3 +. + 0.5 *. a *. b *. s1 -. expo_a *. b *. s2 -. expo_b *. a *. s3 +. 2.0 *. expo_a *. expo_b *. s4 in let s = Array.init 3 f @@ -105,19 +107,19 @@ let contracted_class shell_a shell_b : float Zmap.t = let integral = chop norm (fun () -> k.(0)*.s.(1)*.s.(2) +. s.(0)*.k.(1)*.s.(2) +. - s.(0)*.s.(1)*.k.(2) - ) in + s.(0)*.s.(1)*.k.(2) + ) in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral ) class_indices end done; - let result = + + let result = Zmap.create (Array.length contracted_class) in Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; result end - | None -> Zmap.create 0 (** Create kinetic energy matrix *) @@ -158,7 +160,7 @@ let of_basis basis = result_x.{i_c,j_c} <- value; result_x.{j_c,i_c} <- value; ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) - ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) done; done; Matrix.detri_inplace result; diff --git a/gaussian_integrals/lib/overlap.ml b/gaussian_integrals/lib/overlap.ml index 181bc93..16a0803 100644 --- a/gaussian_integrals/lib/overlap.ml +++ b/gaussian_integrals/lib/overlap.ml @@ -4,11 +4,6 @@ open Qcaml_gaussian_basis open Util open Constants -type t = Matrix.t -external matrix : t -> Matrix.t = "%identity" -external of_matrix : Matrix.t -> t = "%identity" - - module Am = Angular_momentum module Bs = Basis module Co = Coordinate @@ -17,6 +12,9 @@ module Csp = Contracted_shell_pair module Po = Powers module Psp = Primitive_shell_pair +type t = Matrix.t +external matrix : t -> Matrix.t = "%identity" +external of_matrix : Matrix.t -> t = "%identity" let cutoff = integrals_cutoff @@ -80,6 +78,7 @@ let contracted_class shell_a shell_b : float Zmap.t = ) class_indices end ) (Csp.coefs_and_shell_pairs shell_p); + let result = Zmap.create (Array.length contracted_class) in diff --git a/gaussian_integrals/lib/two_electron_integrals.ml b/gaussian_integrals/lib/two_electron_integrals.ml index e851489..0ad9d44 100644 --- a/gaussian_integrals/lib/two_electron_integrals.ml +++ b/gaussian_integrals/lib/two_electron_integrals.ml @@ -4,6 +4,7 @@ open Qcaml_common open Qcaml_linear_algebra open Qcaml_gaussian_basis +open Qcaml_operators open Constants let cutoff = integrals_cutoff @@ -18,7 +19,7 @@ module type Two_ei_structure = sig val name : string val class_of_contracted_shell_pair_couple : - basis:Basis.t -> Cspc.t -> float Zmap.t + ?operator:'a Operator.t -> Cspc.t -> float Zmap.t end @@ -28,34 +29,6 @@ module Make(T : Two_ei_structure) = struct let class_of_contracted_shell_pair_couple = T.class_of_contracted_shell_pair_couple - let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) ~basis shell_pairs = - List.rev_map (fun pair -> - match Cspc.make ~cutoff pair pair with - | Some cspc -> - let cls = class_of_contracted_shell_pair_couple ~basis 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.rev_map fst - - -(* TODO - let filter_contracted_shell_pair_couples - ?(cutoff=integrals_cutoff) shell_pair_couples = - List.rev_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.rev_map fst - -*) - - let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = let to_powers x = let open Zkey in @@ -95,7 +68,7 @@ module Make(T : Two_ei_structure) = struct - let of_basis basis = + let of_basis ?operator basis = let n = Bs.size basis and shell = Bs.contracted_shells basis @@ -113,7 +86,6 @@ module Make(T : Two_ei_structure) = struct let shell_pairs = Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~basis ~cutoff in Printf.printf "%d significant shell pairs computed in %f seconds\n" @@ -156,7 +128,7 @@ module Make(T : Two_ei_structure) = struct match cspc with | Some cspc -> let cls = - class_of_contracted_shell_pair_couple ~basis cspc + class_of_contracted_shell_pair_couple ?operator cspc in store_class ~cutoff eri_array cspc cls | None -> () diff --git a/gaussian_integrals/lib/two_electron_integrals.mli b/gaussian_integrals/lib/two_electron_integrals.mli index a1457f4..60f2709 100644 --- a/gaussian_integrals/lib/two_electron_integrals.mli +++ b/gaussian_integrals/lib/two_electron_integrals.mli @@ -8,6 +8,7 @@ open Qcaml_common open Qcaml_gaussian_basis open Qcaml_linear_algebra +open Qcaml_operators module type Two_ei_structure = sig @@ -15,7 +16,7 @@ val name : string (** Name of the kind of integrals, for printing purposes. *) val class_of_contracted_shell_pair_couple : - basis:Basis.t -> Contracted_shell_pair_couple.t -> float Zmap.t + ?operator:'a 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. *) @@ -27,12 +28,7 @@ module Make : functor (T : Two_ei_structure) -> sig include module type of Four_idx_storage - val filter_contracted_shell_pairs : - ?cutoff:float -> basis:Basis.t -> - Contracted_shell_pair.t list -> Contracted_shell_pair.t list - (** Uses Schwartz screening on contracted shell pairs. *) - - val of_basis : Basis.t -> t + val of_basis : ?operator:'a Operator.t -> Basis.t -> t (** Compute all ERI's for a given {!Basis.t}. *) end diff --git a/gaussian_integrals/lib/two_electron_rr.ml b/gaussian_integrals/lib/two_electron_rr.ml index 6b684c7..4e6c482 100644 --- a/gaussian_integrals/lib/two_electron_rr.ml +++ b/gaussian_integrals/lib/two_electron_rr.ml @@ -302,7 +302,7 @@ let rec hvrr_two_e -let contracted_class_shell_pair_couple ~basis ~zero_m shell_pair_couple : float Zmap.t = +let contracted_class_shell_pair_couple ?operator ~zero_m shell_pair_couple : float Zmap.t = let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in @@ -341,7 +341,7 @@ let contracted_class_shell_pair_couple ~basis ~zero_m shell_pair_couple : float let norm_pq_sq = Co.dot center_pq center_pq in let expo_p_inv = Psp.exponent_inv sp_ab in let expo_q_inv = Psp.exponent_inv sp_cd in - let zero = Zp.zero basis zero_m in + let zero = Zp.zero ?operator zero_m in let zero_m_array = zero_m { zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; @@ -412,7 +412,7 @@ let contracted_class_shell_pair_couple ~basis ~zero_m shell_pair_couple : float -let contracted_class_atomic_shell_pair_couple ~basis ~zero_m atomic_shell_pair_couple : float Zmap.t = +let contracted_class_atomic_shell_pair_couple ?operator ~zero_m atomic_shell_pair_couple : float Zmap.t = let maxm = Am.to_int (Aspc.ang_mom atomic_shell_pair_couple) in @@ -455,7 +455,7 @@ let contracted_class_atomic_shell_pair_couple ~basis ~zero_m atomic_shell_pair_c let norm_pq_sq = Co.dot center_pq center_pq in let expo_q_inv = Psp.exponent_inv sp_cd in - let zero = Zp.zero basis zero_m in + let zero = Zp.zero ?operator zero_m in let zero_m_array = zero_m { zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; diff --git a/gaussian_integrals/lib/two_electron_rr_vectorized.ml b/gaussian_integrals/lib/two_electron_rr_vectorized.ml index 0b1f001..dc01a62 100644 --- a/gaussian_integrals/lib/two_electron_rr_vectorized.ml +++ b/gaussian_integrals/lib/two_electron_rr_vectorized.ml @@ -577,7 +577,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) -let contracted_class_shell_pairs ~basis ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = +let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = let sp = Csp.shell_pairs shell_p and sq = Csp.shell_pairs shell_q @@ -652,7 +652,7 @@ let contracted_class_shell_pairs ~basis ~zero_m ?schwartz_p ?schwartz_q shell_p Co.dot center_pq center_pq in - let zero = Zp.zero basis zero_m in + let zero = Zp.zero ?operator zero_m in let zero_m_array = zero_m {zero with @@ -783,7 +783,7 @@ let contracted_class_shell_pairs ~basis ~zero_m ?schwartz_p ?schwartz_q shell_p let norm_pq_sq = x *. x +. y *. y +. z *. z in - let zero = Zp.zero basis zero_m in + let zero = Zp.zero ?operator zero_m in zero_m {zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; center_pq = Coordinate.make Coordinate.{x ; y ; z} ; diff --git a/gaussian_integrals/lib/zero_m_parameters.ml b/gaussian_integrals/lib/zero_m_parameters.ml index dd9249d..f44f45b 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_gaussian_basis +open Qcaml_operators -type t = +type 'a t = { expo_p_inv : float ; expo_q_inv : float ; @@ -10,14 +10,14 @@ type t = center_pq : Coordinate.t ; center_pa : Coordinate.t ; center_qc : Coordinate.t ; - zero_m_func : t -> float array ; - basis : Basis.t ; + zero_m_func : 'a t -> float array ; + operator : 'a Operator.t option; } -let zero basis zero_m_func = +let zero ?operator zero_m_func = { zero_m_func ; - basis ; + operator ; maxm=0 ; expo_p_inv = 0.; expo_q_inv = 0.; diff --git a/operators/lib/dune b/operators/lib/dune new file mode 100644 index 0000000..2ca47b9 --- /dev/null +++ b/operators/lib/dune @@ -0,0 +1,10 @@ +; name = name of the supermodule that will wrap all source files as submodules +; public_name = name of the library for ocamlfind and opam +(library + (name qcaml_operators) + (public_name qcaml.operators) + (libraries + str + qcaml.common + ) + (synopsis "Parameterized operators")) diff --git a/operators/lib/f12_operator.ml b/operators/lib/f12_operator.ml new file mode 100644 index 0000000..90b1267 --- /dev/null +++ b/operators/lib/f12_operator.ml @@ -0,0 +1,76 @@ +(** Type for f12 correlation factors *) + +type t = + { + expo_s : float ; + gaussian : Gaussian_operator.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 + in + let gaussian = Gaussian_operator.make coef_g expo_sg in + { expo_s ; gaussian } + + +(* -1/expo_s *. exp (-expo_s r) *) +let gaussian_geminal expo_s = + let coef_g = + [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] + |> Array.map (fun x -> -. x /. expo_s) + and expo_sg = + [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] + in + make_gaussian_corr_factor expo_s coef_g expo_sg + + + +(* exp (-expo_s r) *) +let simple_gaussian_geminal expo_s = + let coef_g = + [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] + and expo_sg = + [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] + in + make_gaussian_corr_factor expo_s coef_g expo_sg + + + +(** r12 * exp ( -expo_s * r) *) +let gaussian_geminal_times_r12 expo_s = + let coef_g = + [| 0.2454 ; 0.2938 ; 0.1815 ; 0.11281 ; 0.07502 ; 0.05280 |] + and expo_sg = + [| 0.1824 ; 0.7118; 2.252 ; 6.474 ; 19.66 ; 77.92 |] + in make_gaussian_corr_factor expo_s coef_g expo_sg + + +(* exp (-expo_s r) *) +let simple_gaussian_geminal' expo_s = + let coef_g = + [| + -3.4793465193721626604883567779324948787689208984375 ; + -0.00571703486454788484955047422886309504974633455276489257812 ; + 4.14878218728681513738365538301877677440643310546875 ; + 0.202874298181392742623785352407139725983142852783203125 ; + 0.0819187742387294803858566183407674543559551239013671875 ; + 0.04225945671351955673644695821167260874062776565551757812 ; + |] + and expo_sg = + [| + 0.63172472556807146570889699432882480323314666748046875; + 26.3759196683467962429858744144439697265625; + 0.63172102793029016876147352377302013337612152099609375; + 7.08429025944207335641067402320913970470428466796875; + 42.4442841447001910637482069432735443115234375; + 391.44036073596890901171718724071979522705078125 ; + |] + in make_gaussian_corr_factor expo_s coef_g expo_sg + + + + + + diff --git a/operators/lib/gaussian_operator.ml b/operators/lib/gaussian_operator.ml new file mode 100644 index 0000000..0f47cd5 --- /dev/null +++ b/operators/lib/gaussian_operator.ml @@ -0,0 +1,37 @@ +(** Representation for two-electron operators expressed in a Gaussian basis set. *) + +type t = +{ + coef_g : float array; + expo_g : float array; + expo_g_inv : float array; +} + +let make coef_g 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 ; + 0.564546514023 ; 0.404228610699 ; 0.297458536575 ; 0.223321219537 ; + 0.169933732064 ; 0.130190978230 ; 0.099652303426 ; 0.075428246546 ; + 0.0555635614051 ; 0.0386791283055 ; 0.0237550435652 ; 0.010006278387 ; + |] + and expo_g = + [| 84135.654509 ; 2971.58727634 ; 474.716025959 ; 130.676724560 ; + 47.3938388887 ; 20.2078651631 ; 9.5411021938 ; 4.8109546955 ; + 2.52795733067 ; 1.35894103210 ; 0.73586710268 ; 0.39557629706 ; + 0.20785895177 ; 0.104809693858 ; 0.049485682527 ; 0.021099788990 ; + 0.007652472186 ; 0.0021065225215 ; 0.0003365204879 ; 0.0000118855674 |] + in make coef_g expo_g + + + + + diff --git a/operators/lib/operator.ml b/operators/lib/operator.ml new file mode 100644 index 0000000..86f04e3 --- /dev/null +++ b/operators/lib/operator.ml @@ -0,0 +1,6 @@ +type 'a t = +| F12 of F12_operator.t +| Gaussian of Gaussian_operator.t + +let of_f12 f = F12 f +let of_gaussian g = Gaussian g diff --git a/operators/lib/operator.mli b/operators/lib/operator.mli new file mode 100644 index 0000000..82af419 --- /dev/null +++ b/operators/lib/operator.mli @@ -0,0 +1,9 @@ +type 'a t = +| F12 of F12_operator.t +| Gaussian of Gaussian_operator.t + +val of_f12 : F12_operator.t -> F12_operator.t t +(** Creates a F12 operator *) + +val of_gaussian: Gaussian_operator.t -> Gaussian_operator.t t +(** Creates a Gaussian operator *)