(** 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.range_separation with
      | Some x -> x
      | None   -> 0.5  (* TODO invalid_arg "range_separation 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