open Util
open Constants

exception NullPair

(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *)
let chop f g =
  if (abs_float f) < cutoff then 0.
  else f *. (g ())



(** Horizontal and Vertical Recurrence Relations (HVRR) *)
let hvrr_one_e
    (angMom_a, angMom_b) (totAngMom_a, totAngMom_b)
    (maxm, zero_m_array) (expo_inv_p) (center_ab, center_pa, center_pc)
    map
  =

  let totAngMom_a = Angular_momentum.to_int totAngMom_a
  and totAngMom_b = Angular_momentum.to_int totAngMom_b
  in
  let maxm = totAngMom_a+totAngMom_b in
  let maxsze = maxm+1 in
  let empty = Array.make maxsze 0. in

  (** Vertical recurrence relations *)
  let rec vrr angMom_a totAngMom_a =
    let ax,ay,az = angMom_a in
    if (ax < 0) || (ay < 0) || (az < 0) then
      empty
    else
      match totAngMom_a with
      | 0 -> zero_m_array
      | _ ->
        let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in

        try Zmap.find map key with
        | Not_found -> 
          let result = 
            let am, amm, amxyz, xyz  =
              match angMom_a with
              | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0
              | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1
              | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2
            in
            if amxyz < 0 then empty else
              let f1 = Coordinate.coord center_pa xyz
              and f2 = expo_inv_p *.  (Coordinate.coord center_pc xyz)
              in
              if amxyz < 1 then
                let v1 =
                  vrr am (totAngMom_a-1)
                in
                Array.init maxsze (fun m ->
                    if m = maxm then 0. else (f1 *. v1.(m) ) -. f2 *. v1.(m+1) )
              else
                let v3 =
                  vrr amm (totAngMom_a-2)
                in
                let v1 =
                  vrr am (totAngMom_a-1)
                in
                let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
                Array.init maxsze (fun m ->
                    (if m = maxm then 0. else
                       (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
                    +. f3 *. (v3.(m) +.  if m = maxm then 0. else
                                expo_inv_p *. v3.(m+1))
                  )
          in Zmap.add map key result;
          result


  (** Horizontal recurrence relations *)
  and hrr angMom_a angMom_b totAngMom_a totAngMom_b =

    let bx,by,bz = angMom_b in
    if (bx < 0) || (by < 0) || (bz < 0) then 0.
    else
      match totAngMom_b with
      | 0 -> (vrr angMom_a totAngMom_a).(0)
      | _ ->
        let angMom_ax, angMom_ay, angMom_az = angMom_a
        and angMom_bx, angMom_by, angMom_bz = angMom_b in
        let bxyz, xyz =
          match angMom_b with
          | (_,0,0) -> angMom_bx, 0
          | (_,_,0) -> angMom_by, 1
          | (_,_,_) -> angMom_bz, 2
        in
        if (bxyz < 1) then 0. else
          let ap, bm =
            match xyz with
            | 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz)
            | 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz)
            | _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1)
          in
          let h1 =
            hrr ap bm (totAngMom_a+1) (totAngMom_b-1) 
          in
          let f2 =
            (Coordinate.coord center_ab xyz)
          in
          if (abs_float f2 < cutoff) then h1 else
            let h2 =
              hrr angMom_a bm totAngMom_a (totAngMom_b-1)
            in
            h1 +. f2 *. h2

  in
  hrr
    (angMom_a.(0),angMom_a.(1),angMom_a.(2))
    (angMom_b.(0),angMom_b.(1),angMom_b.(2))
    totAngMom_a totAngMom_b









(** Computes all the one-electron integrals of the contracted shell pair *)
let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =

  let shell_a = shell_p.(0).Shell_pair.shell_a
  and shell_b = shell_p.(0).Shell_pair.shell_b
  in
  let maxm = 
    let open Angular_momentum in
    (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) 
  in

  (* Pre-computation of integral class indices *)
  let class_indices = 
    Angular_momentum.zkey_array
      (Angular_momentum.Doublet
         Contracted_shell.(totAngMom shell_a, totAngMom shell_b))
  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 *)

  for ab=0 to (Array.length shell_p - 1)
  do
    try
      begin
        let coef_prod = shell_p.(ab).Shell_pair.coef in
        let norm_coef_scale_p = shell_p.(ab).Shell_pair.norm_coef_scale in

        (** Screening on the product of coefficients *)
        if (abs_float coef_prod) < 1.e-4*.cutoff then
          raise NullPair;


        let expo_pq_inv =
          shell_p.(ab).Shell_pair.expo_inv 
        in

        let center_ab = 
          shell_p.(ab).Shell_pair.center_ab 
        in
        let center_pa = 
          Coordinate.(center_ab |- shell_a.Contracted_shell.center) 
        in

        for c=0 to Array.length geometry - 1 do
          let element, nucl_coord = geometry.(c) in
          let charge = Element.to_charge element |> Charge.to_float in
          let center_pc =
            Coordinate.(shell_p.(ab).Shell_pair.center |- nucl_coord ) 
          in
          let norm_pq_sq = 
            Coordinate.dot center_pc center_pc
          in

          let zero_m_array = 
            zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
          in
          match Contracted_shell.(totAngMom shell_a, totAngMom shell_b) with
          | Angular_momentum.(S,S) -> 
            let integral =
              zero_m_array.(0)
            in
            contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge
          | _ ->
            let map = Zmap.create (2*maxm) in
            let norm_coef_scale = norm_coef_scale_p in
            (* Compute the integral class from the primitive shell quartet *)
            class_indices
            |> Array.iteri (fun i key ->
                let (angMomA,angMomB) =
                  let a = Zkey.to_int_array Zkey.Kind_6 key in
                  ( [| a.(0) ;  a.(1) ;  a.(2) |], 
                    [| a.(3) ;  a.(4) ;  a.(5) |])
                in
                let norm = norm_coef_scale.(i) in
                let coef_prod = coef_prod *. norm in
                let integral = 
                  hvrr_one_e (angMomA, angMomB)
                    (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b)
                    (maxm, zero_m_array) shell_p.(ab).Shell_pair.expo_inv
                    (center_ab, center_pa, center_pc)
                    map
                in
                contracted_class.(i) <- contracted_class.(i) -. coef_prod *. integral *. charge
              )
      done
  end
  with NullPair -> ()
done;
let result = 
  Zmap.create (Array.length contracted_class)
in
Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
result