open Lacaml.D
open Simulation
open Constants



type t =
  {
    fock     : Mat.t ;
    core     : Mat.t ;
    coulomb  : Mat.t ;
    exchange : Mat.t ;
  }

let make ~density simulation =
  let m_P  = density
  and m_T  = Lazy.force simulation.kin_ints
  and m_V  = Lazy.force simulation.eN_ints
  and m_G  = Lazy.force simulation.ee_ints
  in
  let nBas = Mat.dim1 m_T
  in

  let m_Hc = Mat.add m_T m_V 
  and m_J = Mat.make0 nBas nBas 
  and m_K = Mat.make0 nBas nBas 
  in
  for sigma = 1 to nBas do
    for nu = 1 to nBas do
      for lambda = 1 to nBas do
        let p = m_P.{lambda,sigma} in
        if abs_float p > epsilon then
          for mu = 1 to nu do
            m_J.{mu,nu} <- m_J.{mu,nu} +.  p *. 
                          ERI.get_phys m_G mu lambda nu sigma
          done
      done
    done
  done;
  for nu = 1 to nBas do
    for sigma = 1 to nBas do
      for lambda = 1 to nBas do
        let p = 0.5 *. m_P.{lambda,sigma} in
        if abs_float p > epsilon then
          for mu = 1 to nu do
            m_K.{mu,nu} <- m_K.{mu,nu} -.  p *. 
                          ERI.get_phys m_G mu lambda sigma nu
          done
      done
    done
  done;
  for nu = 1 to nBas do
    for mu = 1 to nu do
      m_J.{nu,mu} <- m_J.{mu,nu};
      m_K.{nu,mu} <- m_K.{mu,nu};
    done
  done;
  { fock = Mat.add m_Hc @@ Mat.add m_J m_K ;
    core = m_Hc ; coulomb = m_J ; exchange = m_K }