open Util

(** (00|00)^m : Fundamental electron repulsion integral
  $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $

  maxm        : Maximum total angular momentum
  expo_pq_inv : $1./p + 1./q$ where $p$ and $q$ are the exponents of $\phi_p$ and $\phi_q$
  norm_pq_sq  : square of the distance between the centers of $\phi_p$ and $\phi_q$
*)
let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq =
  let exp_pq =
    1. /. expo_pq_inv
  in
  let t =
    norm_pq_sq *. exp_pq
  in
  boys_function ~maxm t
  |> Array.mapi (fun m fm ->
      two_over_sq_pi *. (if m mod 2 = 0 then fm else -.fm) *.
      (pow exp_pq m) *. (sqrt exp_pq)
    )


(** Electron-electron repulsion integral *)
let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t =
  TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d