diff --git a/CI/CIMatrixElement.ml b/CI/CIMatrixElement.ml
index e5b65d7..cac190c 100644
--- a/CI/CIMatrixElement.ml
+++ b/CI/CIMatrixElement.ml
@@ -11,6 +11,7 @@ let non_zero integrals degree_a degree_b ki kj =
   let kia = De.alfa ki and kib = De.beta ki
   and kja = De.alfa kj and kjb = De.beta kj
   in
+
   let single h p spin same opposite = 
     let same_spin_mo_list = 
       Sp.to_list same
@@ -24,6 +25,7 @@ let non_zero integrals degree_a degree_b ki kj =
         List.fold_left (fun accu i -> accu +. two_e h i p i spin (Spin.other spin) ) 0. opposite_spin_mo_list
       in (one_e h p spin) +. same_spin +. opposite_spin
   in
+
   let diag_element =
     let mo_a = Sp.to_list kia
     and mo_b = Sp.to_list kib
@@ -56,6 +58,7 @@ let non_zero integrals degree_a degree_b ki kj =
       in
       one +. two
   in
+
   let result =
     match degree_a, degree_b with
     | 1, 1 -> (* alpha-beta double *)
diff --git a/CI/F12CI.ml b/CI/F12CI.ml
index 4363527..f46d311 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -1,6 +1,8 @@
 open Lacaml.D
 
 module Ds = DeterminantSpace
+module De = Determinant
+module Sp = Spindeterminant
 
 type t = 
   {
@@ -22,12 +24,28 @@ let eigensystem t =  Lazy.force t.eigensystem
 let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
   let integrals = [
     let one_e _ _ _ = 0. in
-    let hf12_s, hf12_o = hf12_integrals in
+    let hf12_s, hf12_o, hf12m_s, hf12m_o = hf12_integrals in
+    let kia = De.alfa ki and kib = De.beta ki
+    and kja = De.alfa kj and kjb = De.beta kj
+    in
+    let mo_a =
+      Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja)
+      |> Bitstring.to_list
+    and mo_b =
+      Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb)
+      |> Bitstring.to_list
+    in
     let two_e i j k l s s' =
       if s' = Spin.other s then
-	hf12_o.{i,j,k,l} 
+	hf12_o.{i,j,k,l}  -. 1. *. (
+	(List.fold_left (fun accu m -> accu +. hf12m_o.{m,i,j,k,l}) 0. mo_a) +.
+	(List.fold_left (fun accu m -> accu +. hf12m_o.{m,j,i,l,k}) 0. mo_b) 
+        )
       else
-	hf12_s.{i,j,k,l}
+	hf12_s.{i,j,k,l} -. 1. *. (
+	(List.fold_left (fun accu m -> accu +. hf12m_s.{m,i,j,k,l}) 0. mo_a) +.
+	(List.fold_left (fun accu m -> accu +. hf12m_s.{m,j,i,l,k}) 0. mo_b) 
+        )
     in
     (one_e, two_e)
   ]
@@ -37,6 +55,7 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
 
 
 
+
 let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
 
   if Parallel.master then
diff --git a/MOBasis/HF12.ml b/MOBasis/HF12.ml
index 10c352c..d514ec9 100644
--- a/MOBasis/HF12.ml
+++ b/MOBasis/HF12.ml
@@ -6,6 +6,8 @@ module Fis = FourIdxStorage
 
 type t = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
        * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
+       * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
+       * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
 
 
 let make ~simulation ~mo_basis ~aux_basis_filename () =
@@ -45,9 +47,11 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   (* Compute the <ij|QHF|kl> integrals *)
   if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!";
 
-  let result_s, result_o =
+  let result_s, result_o, resultm_s, resultm_o =
     Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] , 
-    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] 
+    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] , 
+    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num ; mo_num |] , 
+    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num ; mo_num |] 
   in
 
 
@@ -58,6 +62,9 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   let hf_s = Mat.create mo_num mo_num in
   let hf_o = Mat.create mo_num mo_num in
 
+  let hfm_s  = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num mo_num mo_num in
+  let hfm_o  = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num mo_num mo_num in
+
   for a=1 to mo_num do
     for b=1 to mo_num do
      for i=1 to mo_num do
@@ -103,43 +110,79 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
       done;
 
 (*
-      for a=1 to mo_num do
-        for b=mo_num+1 to aux_num do
-         for i=1 to mo_num do
-          if i <> a then
-            h i a b
-         done
+      let h i a b = 
+         h_s.{i, a, b} <- 0. ;
+         h_o.{i, a, b} <- 0.
+      and f a b k =
+         f_s.{a, b, k} <- 0. ;
+         f_o.{a, b, k} <- 0.
+      in
+*)
+
+      for m=1 to mo_num do
+        for a=mo_num+1 to aux_num do
+          for i=1 to mo_num do
+              h i a m ;
+              h i m a
+          done
         done
       done;
 
       for k=1 to mo_num do
-        for b=mo_num+1 to aux_num do
-          for a=1 to mo_num do
-          if k <> a then
-            f a b k
-         done
+        for m=1 to mo_num do
+          for a=mo_num+1 to aux_num do
+            f a m k ;
+            f m a k
+          done
         done
       done;
-*)
-      let h_o =
+
+      let h_o' =
         Bigarray.(reshape (genarray_of_array3 h_o)) [| mo_num ; aux_num*aux_num |]
         |> Bigarray.array2_of_genarray
       in
-      let f_o =
+      let f_o' =
         Bigarray.(reshape (genarray_of_array3 f_o)) [| aux_num*aux_num ; mo_num |] 
         |> Bigarray.array2_of_genarray
       in
-      let h_s =
+      let h_s' =
         Bigarray.(reshape (genarray_of_array3 h_s)) [| mo_num ; aux_num*aux_num |]
         |> Bigarray.array2_of_genarray
       in
-      let f_s =
+      let f_s' =
         Bigarray.(reshape (genarray_of_array3 f_s)) [| aux_num*aux_num ; mo_num |] 
         |> Bigarray.array2_of_genarray
       in
-      let hf_s = gemm ~alpha:1.0 ~c:hf_s h_s f_s in
-      let hf_o = gemm ~alpha:1.0 ~c:hf_o h_o f_o in
-      hf_s, hf_o, j, l
+      let hf_s = gemm ~c:hf_s h_s' f_s' in
+      let hf_o = gemm ~c:hf_o h_o' f_o' in
+
+      let () = 
+        for m=1 to mo_num do
+          let h_o' =
+            Mat.init_cols mo_num aux_num (fun i a -> h_o.{i,m,a})
+          in
+          let f_o' =
+            Mat.init_cols aux_num mo_num (fun a k -> f_o.{m,a,k})
+          in
+          let h_s' =
+            Mat.init_cols mo_num aux_num (fun i a -> h_s.{i,m,a})
+          in
+          let f_s' =
+            Mat.init_cols aux_num mo_num (fun a k -> f_s.{m,a,k})
+          in
+          let r_s, r_o =
+            gemm h_s' f_s' , 
+            gemm h_o' f_o'
+          in
+          for k = 1 to mo_num do
+            for i = 1 to mo_num do
+              hfm_s.{m,i,k} <- r_s.{i,k};
+              hfm_o.{m,i,k} <- r_o.{i,k}
+            done
+          done
+        done
+      in
+      hf_s, hf_o, hfm_s, hfm_o, j, l
   in
 
   let tasks = 
@@ -154,14 +197,15 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
 
 
   Farm.run ~f:task ~ordered:true tasks
-  |> Stream.iter (fun (hf_s, hf_o, j, l) ->
-    (*
-      Printf.printf "%d %d\n" j l ;
-    *)
+  |> Stream.iter (fun (hf_s, hf_o, hfm_s, hfm_o, j, l) ->
     for k=1 to mo_num do
       for i=1 to mo_num do
-        result_s.{i,k,j,l} <- hf_s.{i,k} ;
-        result_o.{i,k,j,l} <- hf_o.{i,k}
+        result_s.{i,j,k,l} <- hf_s.{i,k} ;
+        result_o.{i,j,k,l} <- hf_o.{i,k} ;
+        for m=1 to mo_num do
+          resultm_s.{m,i,j,k,l} <- hfm_s.{m,i,k} ;
+          resultm_o.{m,i,j,k,l} <- hfm_o.{m,i,k}
+        done
       done
     done );
 
@@ -180,7 +224,7 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   Printf.printf "%!";
   *)
 
-  Parallel.broadcast (lazy (result_s, result_o) )
+  Parallel.broadcast (lazy (result_s, result_o, resultm_s, resultm_o) )