QCaml/HartreeFock/Fock.ml

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 *. 
                          (ERI.get_phys m_G mu lambda nu sigma
                           -. 0.5 *. 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_F.{nu,mu} <- m_F.{mu,nu}
    done
  done;
  m_F