open Lacaml.D
open Simulation
open Constants

type t = 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_F = Mat.add m_T m_V 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_F.{mu,nu} <- m_F.{mu,nu} +.  p *. 
                          (m_G.{mu,lambda,nu,sigma} -. 0.5 *. m_G.{mu,lambda,sigma,nu})
          done
      done
    done
  done;
  for nu = 1 to nBas do
    for mu = 1 to nu do
      m_F.{nu,mu} <- m_F.{mu,nu}
    done
  done;
  m_F