open Lacaml.D

type t = 
  {
    mo_basis  : MOBasis.t ;
    aux_basis : MOBasis.t ;
    det_space : DeterminantSpace.t ;
    ci        : CI.t ;
    f12_amplitudes : Mat.t;
  }

let ci t        =  t.ci
let mo_basis t  =  t.mo_basis
let det_space t =  t.det_space
let mo_class t  =  DeterminantSpace.mo_class @@ det_space t


let f12_integrals mo_basis =
  let two_e_ints  = MOBasis.two_e_ints mo_basis in
  ( (fun i j _ -> 0.),
    (fun i j k l s s' ->
       if s' = Spin.other s then
         ERI.get_phys two_e_ints i j k l
       else
         (ERI.get_phys two_e_ints i j k l) -.
         (ERI.get_phys two_e_ints i j l k)
    ) )




let h_ij mo_basis ki kj =
  let integrals =
    List.map (fun f -> f mo_basis)
      [ CI.h_integrals ]
  in
  CIMatrixElement.make integrals ki kj 
  |> List.hd


let f_ij mo_basis ki kj =
  let integrals =
    List.map (fun f -> f mo_basis)
      [ f12_integrals ]
  in
  CIMatrixElement.make integrals ki kj 
  |> List.hd


let dressing_vector ci f12_amplitudes =

  let mo_basis  = DeterminantSpace.mo_basis ci.CI.det_space in

  let i_o1_alfa = h_ij mo_basis in
  let alfa_o2_i = f_ij mo_basis in

  let w_alfa _ _ = 1. in

  let mo_class = CI.mo_class ci in

  let list_holes = List.concat
    [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
  and list_particles = MOClass.auxiliary_mos mo_class
  in
  CI.second_order_sum  ci  list_holes  list_particles
      i_o1_alfa  alfa_o2_i  w_alfa  f12_amplitudes
  |> Vec.of_list





let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =

  let mo_num = MOBasis.size mo_basis in

  (* Add auxiliary basis set *)
  let s =
    let charge        = Charge.to_int @@ Simulation.charge simulation 
    and multiplicity  = Electrons.multiplicity @@ Simulation.electrons simulation
    and nuclei        = Simulation.nuclei simulation
    in
    let general_basis = 
       Basis.general_basis @@ Simulation.basis simulation
    in
    GeneralBasis.combine [
         general_basis ; GeneralBasis.read aux_basis_filename
         ]
    |> Basis.of_nuclei_and_general_basis nuclei
    |> Simulation.make ~charge ~multiplicity ~nuclei 
  in

  let aux_basis = 
    MOBasis.of_mo_basis s mo_basis
  in

  let det_space = 
     DeterminantSpace.fci_f12_of_mo_basis  aux_basis  ~frozen_core  mo_num
  in

  let ci = CI.make det_space in

  let ci_coef, ci_energy = Parallel.broadcast ci.eigensystem in

  let f12_amplitudes = 
    (* While in a sequential region, initiate the parallel
       4-idx transformation to avoid nested parallel jobs
    *)
    ignore @@ MOBasis.f12_ints mo_basis;

    let f = fun ki kj ->
         if ki <> kj then
            f_ij mo_basis ki kj
         else
            f_ij mo_basis ki kj +. 1.
    in
    let m_F = 
      CI.create_matrix_spin f det_space
      |> Lazy.force
    in
    Matrix.ax_eq_b (Matrix.dense_of_sparse m_F) (Matrix.dense_of_mat ci_coef)
    |> Matrix.to_mat
  in


  { mo_basis ; aux_basis ; det_space ; ci ; f12_amplitudes }