Merge branch 'master' of gitlab.com:scemama/QCaml

Basis/ERI_lr.ml

@@ -0,0 +1,64 @@
+(** Long-range electron-electron repulsion integrals.
+    See Eq(52) in 10.1039/b605188j
+*)
+
+open Constants
+open Util
+
+module Csp = ContractedShellPair
+module Cspc = ContractedShellPairCouple
+
+module T = struct
+
+  let name = "Long-range electron repulsion integrals"
+
+  open Zero_m_parameters
+
+  let zero_m z =
+    let mu_erf = 
+      match z.mu_erf with
+      | Some x -> x
+      | None   -> 0.5  (*TODO invalid_arg "mu_erf is None"*)
+    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 > integrals_cutoff then
+        z.norm_pq_sq *. fG
+        else 0.
+    in
+    let maxm = z.maxm in
+    let result = 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 = two_over_sq_pi *. (sqrt fG) in
+    aux f 0 maxm;
+    result
+
+  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_pair_couple
+        ~zero_m shell_pair_couple
+    else
+      TwoElectronRRVectorized.contracted_class_shell_pairs
+        ~zero_m shell_p shell_q
+
+end
+
+module M = TwoElectronIntegrals.Make(T) 
+include M
+