diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index 4d58689..c517f84 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -6,6 +6,16 @@ open Bigarray
 
 type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t
 
+let get ~r1 ~r2 t =
+  let i,k = r1
+  and j,l = r2
+  in
+  t.{i,k,j,l}
+
+let get_chem t i j k l = get ~r1:(i,j) ~r2:(k,l) t
+let get_phys t i j k l = get ~r1:(i,k) ~r2:(j,l) t
+
+
 module Am  = AngularMomentum
 module As  = AtomicShell
 module Asp = AtomicShellPair
diff --git a/Basis/ERI.mli b/Basis/ERI.mli
new file mode 100644
index 0000000..640e9d5
--- /dev/null
+++ b/Basis/ERI.mli
@@ -0,0 +1,29 @@
+type t 
+val get : r1:int * int -> r2:int * int -> t -> float
+val get_chem : t -> int -> int -> int -> int -> float
+val get_phys : t -> int -> int -> int -> int -> float
+
+val zero_m : maxm:int -> expo_pq_inv:float -> norm_pq_sq:float -> float array
+
+val contracted_class_shell_pairs :
+  ?schwartz_p:float Zmap.t ->
+  ?schwartz_q:float Zmap.t ->
+  ContractedShellPair.t ->
+  ContractedShellPair.t -> float Zmap.t
+
+val contracted_class_shell_pairs_vec :
+  ?schwartz_p:float Zmap.t ->
+  ?schwartz_q:float Zmap.t ->
+  ContractedShellPair.t ->
+  ContractedShellPair.t -> float Zmap.t
+
+val contracted_class_atomic_shell_pairs :
+  ?schwartz_p:float Zmap.t ->
+  ?schwartz_q:float Zmap.t ->
+  AtomicShellPair.t ->
+  AtomicShellPair.t -> float Zmap.t
+
+val of_basis : Basis.t -> t
+
+val to_file : filename:string -> t -> unit
+
diff --git a/Basis/PrimitiveShellPairCouple.ml b/Basis/PrimitiveShellPairCouple.ml
index 70dc88c..cf07d73 100644
--- a/Basis/PrimitiveShellPairCouple.ml
+++ b/Basis/PrimitiveShellPairCouple.ml
@@ -69,3 +69,5 @@ let norm_scales t =
     norm_coef_scale_p_list
   |> Array.concat
 
+
+
diff --git a/HartreeFock/Fock.ml b/HartreeFock/Fock.ml
index b62184d..62646eb 100644
--- a/HartreeFock/Fock.ml
+++ b/HartreeFock/Fock.ml
@@ -2,6 +2,8 @@ open Lacaml.D
 open Simulation
 open Constants
 
+
+
 type t = Mat.t
 
 let make ~density simulation =
@@ -21,7 +23,8 @@ let make ~density simulation =
         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})
+                          (ERI.get_chem m_G mu lambda nu sigma
+                           -. 0.5 *. ERI.get_chem m_G mu lambda sigma nu)
           done
       done
     done