open Util
open Lacaml.D

module Si = Simulation
module El = Electrons
module Ao = AOBasis
module Ov = Overlap

let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1)
  ?threshold_SCF:(threshold_SCF=1.e-5) simulation =

  (* Number of occupied MOs *)
  let nocc =
    simulation.Si.electrons.El.n_alpha
  in

  let nuclear_repulsion = 
    simulation.Si.nuclear_repulsion
  in
  
  let ao_basis = 
    simulation.Si.ao_basis
  in

  (* Initial guess *)
  let guess = 
    Guess.make ~nocc ~guess ao_basis
  in

  (* Orthogonalization matrix *)
  let m_X = 
    Lazy.force ao_basis.Ao.ortho
  in


  (* Overlap matrix *)
  let m_S =
    Lazy.force  ao_basis.Ao.overlap
    |> Ov.matrix
  in

  let m_T = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix
  and m_V = Lazy.force ao_basis.Ao.eN_ints  |> NucInt.matrix
  in

  (* Level shift in MO basis *)
  let m_LSmo =
    Array.init (Mat.dim2 m_X) (fun i ->
      if i > nocc then level_shift else 0.)
    |> Vec.of_array
    |> Mat.of_diag
  in

  (* SCF iterations *)
  let rec loop nSCF iterations energy_prev m_C diis =

      (* Density matrix over nocc occupied MOs *)
      let m_P =
        gemm  ~alpha:2.  ~transb:`T  ~k:nocc  m_C  m_C
      in

      (* Fock matrix in AO basis *)
      let m_F, m_Hc, m_J, m_K =
        let x = 
          Fock.make  ~density:m_P  ao_basis
        in
        x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange
      in
(*
debug_matrix "Fock" m_F;
debug_matrix "Coulomb" m_J;
debug_matrix "Exchange" m_K;
debug_matrix "HCore" m_Hc;
*)
      (* Add level shift in AO basis *)
      let m_F =
        let m_SC =
          gemm  m_S  m_C
        in
        gemm  m_SC (gemm  m_LSmo  m_SC  ~transb:`T)
        |> Mat.add m_F
      in


      (* Fock matrix in orthogonal basis *)
      let m_F_ortho =
        xt_o_x m_F m_X
      in

      let error_fock = 
        let fps =
          gemm  m_F  (gemm  m_P  m_S)
        and spf =
          gemm  m_S  (gemm  m_P  m_F)
        in
        xt_o_x  (Mat.sub  fps  spf) m_X
      in

      let diis = 
        DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
      in

      let m_F_diis =
        let x = 
          Bigarray.genarray_of_array1 (DIIS.next diis)
        in
        Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
      in


      (* MOs in orthogonal MO basis *)
      let m_C', eigenvalues =
        diagonalize_symm  m_F_diis
      in

      (* MOs in AO basis *)
      let m_C =
        gemm  m_X  m_C'
      in

      (* Hartree-Fock energy *)
      let energy =
        nuclear_repulsion +. 0.5 *.
          Mat.gemm_trace  m_P  (Mat.add  m_Hc  m_F)
      in

      (* Convergence criterion *)
      let error =
        error_fock 
        |> Mat.as_vec 
        |> amax
        |> abs_float
      in

      let converged =
        nSCF = max_scf  || error < threshold_SCF
      in

      let gap =
        eigenvalues.{nocc+1} -. eigenvalues.{nocc};
      in

      let () = 
        match energy_prev with
        | Some energy_prev ->
          Printf.eprintf "%3d  %16.10f  %16.10f %11.4e  %10.4f\n%!" nSCF energy (energy -. energy_prev) error gap
        | None -> 
          Printf.eprintf "%3d  %16.10f  %16s %11.4e  %10.4f\n%!" nSCF energy "" error gap
      in

      if not converged then
          loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C diis
      else
          let iterations = 
            List.rev ( (energy, error, gap) :: iterations )
            |> Array.of_list
          in
          { HartreeFock_type.
            nocc;
            guess ;
            eigenvectors = m_C ;
            eigenvalues ;
            energy ;
            nuclear_repulsion;
            iterations ;
            kin_energy = Mat.gemm_trace m_P m_T;
            eN_energy = Mat.gemm_trace m_P m_V;
            coulomb_energy  = 0.5 *. Mat.gemm_trace m_P m_J;
            exchange_energy = 0.5 *. Mat.gemm_trace m_P m_K;
          }
  in


  (* Guess coefficients *)
  let m_H = 
    match guess with
    | Guess.Hcore m_H -> m_H 
    | Guess.Huckel m_H -> m_H 
  in
  let m_Hmo =
    xt_o_x  m_H  m_X
  in
  let m_C', _ =
    diagonalize_symm  m_Hmo
  in
  let m_C =
    gemm  m_X  m_C'
  in

  let diis = DIIS.make () in
  loop 1 [] None m_C diis