open Util
open Constants
open Powers
open Coordinate

let debug=false

let cutoff2 = cutoff *. cutoff

exception NullQuartet

(** Horizontal and Vertical Recurrence Relations (HVRR) *)
let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
    zero_m_array
    (expo_b, expo_d)
    (expo_inv_p, expo_inv_q)
    (center_ab, center_cd, center_pq)
    map_1d map_2d map_1d' map_2d'
  =

  (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *)
  if  angMom_a.tot + angMom_b.tot < angMom_c.tot + angMom_d.tot then
    hvrr_two_e (angMom_c, angMom_d, angMom_a, angMom_b)
      zero_m_array
      (expo_d, expo_b)
      (expo_inv_q, expo_inv_p)
      (center_cd, center_ab, (Coordinate.neg center_pq) )
      map_1d' map_2d' map_1d map_2d

  else

    let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in
    let maxsze = maxm+1 in

    let empty = Array.make (maxm+1) 0.  in

    let get_xyz angMom = 
      match angMom with
      | { y=0 ; z=0 ; _ } -> X
      | {       z=0 ; _ } -> Y
      | _                 -> Z
    in


    (** Vertical recurrence relations *)
    let rec vrr0 angMom_a = 

      match angMom_a.tot with 
      | 0 -> zero_m_array
      | _ ->
        let key = Zkey.of_powers (Zkey.Three angMom_a) in

        try Zmap.find map_1d key with
        | Not_found ->
          let result = 
            let xyz = get_xyz angMom_a in
            let am    = Powers.decr xyz angMom_a in
            let amxyz = Powers.get xyz am in

            if amxyz >= 0 then
                let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) 
                and f2 = expo_b *. expo_inv_p *.  (Coordinate.get xyz center_ab)
                in
                let result = Array.create_float maxsze in
                if amxyz < 1 then
                    begin
                      let v1 = 
                        vrr0 am
                      in
                      for m=0 to maxm-1 do
                          result.(m) <- f1 *. v1.(m+1)  -. f2 *. v1.(m) 
                      done;
                      result.(maxm) <- -. f2 *. v1.(maxm) 
                    end
                else
                    begin
                      let amm = Powers.decr xyz am in
                      let v3 = vrr0 amm in
                      let v1 = vrr0 am in
                      let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
                      for m=0 to maxm-1 do
                        result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) 
                          +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1))
                      done;
                      result.(maxm) <- f3 *. v3.(maxm) 
                    end;
                result
            else
                empty
          in Zmap.add map_1d key result;
          result


    and vrr angMom_a angMom_c =

      match (angMom_a.tot, angMom_c.tot) with
      | (i,0) -> if (i>0) then vrr0 angMom_a
                 else zero_m_array 
      | (_,_) ->
        let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in

        try Zmap.find map_2d key with
        | Not_found -> 
          let result = 
            let xyz = get_xyz angMom_c in
            let cm = Powers.decr xyz angMom_c in
            let cmxyz = Powers.get xyz cm in
            let axyz  = Powers.get xyz angMom_a in

            if cmxyz >= 0 then
                let f1 =
                  -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) 
                and f2 = 
                  expo_inv_q *. (Coordinate.get xyz center_pq)
                in
                let result = Array.make maxsze 0. in
                if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then
                    begin
                      let v1 = 
                        vrr angMom_a cm 
                      in
                      for m=0 to maxm-1 do
                        result.(m) <- f1 *. v1.(m) -.  f2 *. v1.(m+1) ;
                      done;
                      result.(maxm) <- f1 *. v1.(maxm) ;
                    end;
                if cmxyz > 0 then 
                    begin
                      let f3 =
                        (float_of_int cmxyz) *. expo_inv_q *. 0.5 
                      in
                      if (abs_float f3 > cutoff) ||
                          (abs_float (f3 *. expo_inv_q) > cutoff) then 
                        begin
                          let v3 = 
                            let cmm = Powers.decr xyz cm in
                            vrr angMom_a cmm 
                          in
                          for m=0 to maxm-1 do
                            result.(m) <- result.(m) +.
                                f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1))
                          done;
                          result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm)
                        end
                    end;
                if axyz > 0 && cmxyz >= 0 then 
                    begin
                      let am = Powers.decr xyz angMom_a in
                      let f5 = 
                        (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 
                      in
                      if (abs_float f5 > cutoff) then 
                        let v5 =
                          vrr am cm 
                        in
                        for m=0 to maxm-1 do
                          result.(m) <- result.(m) -.  f5 *. v5.(m+1)
                        done
                    end;
                result
            else
                empty
          in Zmap.add map_2d key result;
          result





    (** Horizontal recurrence relations *)
    and hrr0 angMom_a angMom_b angMom_c =

      match angMom_b.tot with
      | 0 -> (vrr angMom_a angMom_c).(0)
      | 1 ->
        let xyz  = get_xyz angMom_b in
        let ap = Powers.incr xyz angMom_a in
        let v1 = vrr ap angMom_c in
        let f2 = Coordinate.get xyz center_ab in
        if (abs_float f2 < cutoff) then v1.(0) else
            let v2 = vrr angMom_a angMom_c in
            v1.(0) +. f2 *. v2.(0)
      | _ ->
        let xyz  = get_xyz angMom_b in
        let bxyz = Powers.get xyz angMom_b in
        if bxyz > 0 then 
            let ap = Powers.incr xyz angMom_a in
            let bm = Powers.decr xyz angMom_b in
            let h1 = hrr0 ap bm angMom_c in
            let f2 = Coordinate.get xyz center_ab in
            if abs_float f2 < cutoff then h1 else
                let h2 = hrr0 angMom_a bm angMom_c in
                h1 +. f2 *. h2
        else 0.


    and hrr angMom_a angMom_b angMom_c angMom_d =

      match (angMom_b.tot, angMom_d.tot) with
      | (_,0) ->
        if (angMom_b.tot = 0) then 
            (vrr angMom_a angMom_c).(0)
        else
            hrr0 angMom_a angMom_b angMom_c
      | (_,_) ->
        let xyz  = get_xyz angMom_d in
        let cp = Powers.incr xyz angMom_c in
        let dm = Powers.decr xyz angMom_d in 
        let h1 = hrr angMom_a angMom_b cp dm in
        let f2 = Coordinate.get xyz center_cd in
        if abs_float f2 < cutoff then h1 else
            let h2 = hrr angMom_a angMom_b angMom_c dm in
            h1 +. f2 *. h2

    in
    hrr angMom_a angMom_b angMom_c angMom_d




let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =

  let shell_a = shell_p.ContractedShellPair.shell_a
  and shell_b = shell_p.ContractedShellPair.shell_b
  and shell_c = shell_q.ContractedShellPair.shell_a
  and shell_d = shell_q.ContractedShellPair.shell_b
  and sp = shell_p.ContractedShellPair.shell_pairs
  and sq = shell_q.ContractedShellPair.shell_pairs
  in
  let maxm = 
    shell_p.ContractedShellPair.totAngMomInt +
    shell_q.ContractedShellPair.totAngMomInt 
  in

  (* Pre-computation of integral class indices *)
  let class_indices = 
    Angular_momentum.zkey_array
      (Angular_momentum.Quartet
         Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
                           totAngMom shell_c, totAngMom shell_d))
  in

  let contracted_class = 
    Array.make (Array.length class_indices) 0.;
  in

  let monocentric =
      shell_p.ContractedShellPair.monocentric &&
      shell_q.ContractedShellPair.monocentric &&
      Contracted_shell.center shell_p.ContractedShellPair.shell_a = 
      Contracted_shell.center shell_q.ContractedShellPair.shell_a
  in

  (* Compute all integrals in the shell for each pair of significant shell pairs *)

  let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
  let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in
  for ab=0 to (Array.length sp - 1) do
    let cab = shell_p.ContractedShellPair.coef.(ab)  in
    let b = sp.(ab).ShellPair.j in
    for cd=0 to (Array.length shell_q.ContractedShellPair.shell_pairs - 1) do
      let coef_prod =
        cab *. shell_q.ContractedShellPair.coef.(cd)
      in
      (** Screening on the product of coefficients *)
      try
        if (abs_float coef_prod) < 1.e-3*.cutoff then
          raise NullQuartet;

        let center_pq = 
          sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center
        in
        let norm_pq_sq =
          Coordinate.dot center_pq center_pq
        in

        let expo_pq_inv =
          shell_p.ContractedShellPair.expo_inv.(ab) +.
          shell_q.ContractedShellPair.expo_inv.(cd) 
        in

        let zero_m_array = 
          zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
        in
        begin
          match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
                                  totAngMom shell_c, totAngMom shell_d) with
          | Angular_momentum.(S,S,S,S) ->
            let integral =
              zero_m_array.(0)
            in
            contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
          | _ -> 
            let d = shell_q.ContractedShellPair.shell_pairs.(cd).ShellPair.j in
            let map_1d = Zmap.create (4*maxm) in
            let map_2d = Zmap.create (Array.length class_indices) in
            let map_1d' = Zmap.create (4*maxm) in
            let map_2d' = Zmap.create (Array.length class_indices) in
            let norm_coef_scale =
              Array.to_list norm_coef_scale_p
              |> List.map (fun v1 ->
                  Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) 
              |> Array.concat
            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 ~kind:Zkey.Kind_12 key with
                  | Zkey.Twelve x -> x
                  | _             -> assert false
                in
                try
                  if monocentric then
                    begin
                      if (  ((1 land angMom_a.x+angMom_b.x+angMom_c.x+angMom_d.x)=1) ||
                            ((1 land angMom_a.y+angMom_b.y+angMom_c.y+angMom_d.y)=1) ||
                            ((1 land angMom_a.z+angMom_b.z+angMom_c.z+angMom_d.z)=1) 
                          ) then
                        raise NullQuartet
                    end;

                  (*
                  (* Schwartz screening *)
                  if (maxm > 2) then
                  (
                    let schwartz_p = 
                      let key = Zkey.of_int_tuple (Zkey.Twelve
                      (angMomA, angMomB, angMomA, angMomB) )
                      in
                      match schwartz_p with
                      | None -> 1.
                      | Some schwartz_p -> Zmap.find schwartz_p key
                    in
                    if schwartz_p < cutoff then raise NullQuartet;
                    let schwartz_q = 
                      let key = Zkey.of_int_tuple (Zkey.Twelve
                      (angMomC, angMomD, angMomC, angMomD) ) 
                      in
                      match schwartz_q with
                      | None -> 1.
                      | Some schwartz_q -> Zmap.find schwartz_q key
                    in
                    if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet;
                    );
                  *)


                  let norm =  norm_coef_scale.(i) in
                  let coef_prod = coef_prod *. norm in
                  let integral = 
                    hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
                      zero_m_array
                      (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
                      (shell_p.ContractedShellPair.expo_inv.(ab), 
                       shell_q.ContractedShellPair.expo_inv.(cd) )
                      (sp.(ab).ShellPair.center_ab, sq.(cd).ShellPair.center_ab, center_pq)
                      map_1d map_2d map_1d' map_2d'
                  in
                  contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
                with NullQuartet -> ()
              )
        end
      with NullQuartet -> ()
    done
  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



(** Computes all the two-electron integrals of the contracted shell quartet *)
let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =

  let shell_p = ContractedShellPair.create ~cutoff shell_a shell_b
  and shell_q = ContractedShellPair.create ~cutoff shell_c shell_d
  in
  contracted_class_shell_pairs ~zero_m shell_p shell_q