diff --git a/Basis/AOBasis.mli b/Basis/AOBasis.mli
index 1369dfb..c680f4f 100644
--- a/Basis/AOBasis.mli
+++ b/Basis/AOBasis.mli
@@ -1,6 +1,6 @@
 (** Data structure for Atomic Orbitals. *)
 
-type t =
+type t = private
 {
   basis      : Basis.t ;                      (* One-electron basis set *)
   overlap    : Overlap.t lazy_t;              (* Overlap matrix *)
diff --git a/Nuclei/Element.ml b/Nuclei/Element.ml
index 50275e1..c0b8a43 100644
--- a/Nuclei/Element.ml
+++ b/Nuclei/Element.ml
@@ -185,3 +185,21 @@ let mass c =
   end
   |> Mass.of_float
 
+
+let small_core = function
+  |  X   ->  0     |  H   ->  0     |  He  ->  0     |  Li  ->  2
+  |  Be  ->  2     |  B   ->  2     |  C   ->  2     |  N   ->  2   
+  |  O   ->  2     |  F   ->  2     |  Ne  ->  2     |  Na  ->  10
+  |  Mg  ->  10    |  Al  ->  10    |  Si  ->  10    |  P   ->  10  
+  |  S   ->  10    |  Cl  ->  10    |  Ar  ->  10    |  K   ->  10  
+  |  Ca  ->  10    |  Sc  ->  10    |  Ti  ->  10    |  V   ->  10
+  |  Cr  ->  10    |  Mn  ->  10    |  Fe  ->  10    |  Co  ->  10   
+  |  Ni  ->  10    |  Cu  ->  10    |  Zn  ->  10    |  Ga  ->  10   
+  |  Ge  ->  10    |  As  ->  10    |  Se  ->  10    |  Br  ->  10   
+  |  Kr  ->  10    |  Rb  ->  28    |  Sr  ->  28    |  Y   ->  28   
+  |  Zr  ->  28    |  Nb  ->  28    |  Mo  ->  28    |  Tc  ->  28   
+  |  Ru  ->  28    |  Rh  ->  28    |  Pd  ->  28    |  Ag  ->  28   
+  |  Cd  ->  28    |  In  ->  28    |  Sn  ->  28    |  Sb  ->  28   
+  |  Te  ->  28    |  I   ->  28    |  Xe  ->  28    |  Pt  ->  60
+
+
diff --git a/Nuclei/Element.mli b/Nuclei/Element.mli
index 7bf9ee7..9b025d8 100644
--- a/Nuclei/Element.mli
+++ b/Nuclei/Element.mli
@@ -63,4 +63,10 @@ val vdw_radius : t -> NonNegativeFloat.t
 val mass : t -> Mass.t
 (** Atomic mass of the elements, in atomic units. *)
 
+val small_core : t -> int
+(** Number of electrons in the small core model (all except the outermost two shells). *)
 
+(* TODO
+val large_core : t -> int
+(** Number of electrons in the large core model (all except the outermost shell). *)
+*)
diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml
index 1c4c14f..3fca48f 100644
--- a/Nuclei/Nuclei.ml
+++ b/Nuclei/Nuclei.ml
@@ -123,3 +123,8 @@ let to_t2_string t =
     |> Array.to_list )
   |> String.concat "\n"
 
+
+let small_core a = 
+  Array.fold_left (fun accu (e,_) -> accu + (Element.small_core e)) 0 a
+
+
diff --git a/Nuclei/Nuclei.mli b/Nuclei/Nuclei.mli
index abc9f59..0e27337 100644
--- a/Nuclei/Nuclei.mli
+++ b/Nuclei/Nuclei.mli
@@ -27,5 +27,8 @@ val repulsion : t -> float
 val charge : t -> Charge.t
 (** Sum of the charges of the nuclei. *)
 
+val small_core : t -> int
+(** Number of core electrons in the small core model. *)
+
 val to_xyz_string  : t -> string
 val to_t2_string  : t -> string
diff --git a/SCF/Fock.ml b/SCF/Fock.ml
index ea69b93..b55c041 100644
--- a/SCF/Fock.ml
+++ b/SCF/Fock.ml
@@ -29,57 +29,31 @@ let make ~density ao_basis =
   and m_J = Array.make_matrix nBas nBas 0.
   and m_K = Array.make_matrix nBas nBas 0.
   in
-(*
+  (*
   let permutations  i j k l  =
-    let result =
-      [  [| i ; j ; k ; l |] ; 
-         [| i ; l ; k ; j |] ; 
-         [| k ; j ; i ; l |] ; 
-         [| k ; l ; i ; j |] ; 
-         [| j ; i ; l ; k |] ; 
-         [| j ; k ; l ; i |] ; 
-         [| l ; i ; j ; k |] ; 
-         [| l ; k ; j ; i |] ] 
-    in
-    if i<>k && j<>l && i<>j && k<>l then
-      result
-    else
-      List.map Zkey.of_int_array result
-      |> List.sort_uniq Pervasives.compare 
-      |> List.map Zkey.to_int_array
+      [  [ i ; j ; k ; l ] ; 
+         [ k ; j ; i ; l ] ; 
+         [ i ; l ; k ; j ] ; 
+         [ k ; l ; i ; j ] ; 
+         [ j ; i ; l ; k ] ; 
+         [ j ; k ; l ; i ] ; 
+         [ l ; i ; j ; k ] ; 
+         [ l ; k ; j ; i ]
+         ] 
   in
-  let sum = ref 0 in
   ERI.to_stream m_G
   |> Stream.iter (fun { ERI.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } ->
      permutations i_r1 j_r2 k_r1 l_r2 
      |> List.iter ( fun ijkl ->
-        let mu     = ijkl.(0)
-        and nu     = ijkl.(2)
-        and lambda = ijkl.(1)
-        and sigma  = ijkl.(3)
-        in
-        incr sum;
-        let p = m_P.{lambda,sigma} in
-        if abs_float p > epsilon then
-          let m_Jnu = m_J.(nu-1) in
-          m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. p *. value )
-     );
-
-  Printf.printf "%d %d\n%!" !sum (nBas*nBas*nBas*nBas);
-*)
-  for sigma = 1 to nBas do
-    for nu = 1 to nBas do
-      let m_Jnu = m_J.(nu-1) in
-      for lambda = 1 to nBas do
-        let p = m_P.{lambda,sigma} in
-        for mu = 1 to nBas do
-          m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  p *. 
-                        ERI.get_phys m_G mu lambda nu sigma
-        done
-      done
-    done
-  done;
-    (*
+        match ijkl with
+        | mu :: lambda :: nu :: sigma :: [] ->
+          let p = m_P.{lambda,sigma} in
+          if abs_float p > epsilon then
+            let m_Jnu = m_J.(nu-1) in
+            m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. p *. value 
+        | _ -> assert false
+    ));
+  *)
   for sigma = 1 to nBas do
     for nu = 1 to nBas do
       let m_Jnu = m_J.(nu-1) in
@@ -103,7 +77,6 @@ let make ~density ao_basis =
       m_J.(mu-1).(nu-1) <- m_J.(nu-1).(mu-1);
     done
   done;
-     *)
 
   for nu = 1 to nBas do
     let m_Knu = m_K.(nu-1) in
diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml
index de54bd8..17edd2b 100644
--- a/Utils/FourIdxStorage.ml
+++ b/Utils/FourIdxStorage.ml
@@ -173,7 +173,12 @@ let to_stream d =
   and k = ref 1
   and l = ref 1
   in
-  let f_dense _ =
+  let f i k = 
+    let p, r = 
+      if i <= k then i, k else k, i
+    in p+ (r*r-r)/2
+  in
+  let rec f_dense _ =
     i := !i+1;
     if !i > !k then begin
       i := 1;
@@ -181,7 +186,7 @@ let to_stream d =
       if !j > !l then begin
         j := 1;
         k := !k + 1;
-        if !k > !l then begin
+        if !k > d.size then begin
           k := 1;
           l := !l + 1;
         end;