diff --git a/CI/CI.ml b/CI/CI.ml
index ab1ef95..2a87baa 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -290,6 +290,7 @@ let create_matrix_spin ?(nmax=2) f det_space =
                     let singles = 
                       List.filter (fun (i,d,det_j) -> d < 2) doubles
                     in
+
                     let triples = 
                       List.map (fun (i,_,det_j) -> (i,det_j)) triples
                     in
diff --git a/CI/CIMatrixElementF12.ml b/CI/CIMatrixElementF12.ml
deleted file mode 100644
index 47019f1..0000000
--- a/CI/CIMatrixElementF12.ml
+++ /dev/null
@@ -1,147 +0,0 @@
-open Lacaml.D
-
-module De = Determinant
-module Ex = Excitation
-module Sp = Spindeterminant
-
-type t = float list
-
-
-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 diag_element =
-    let mo_a = Sp.to_list kia
-    and mo_b = Sp.to_list kib
-    in
-    fun one_e two_e ->
-      let one =
-        (List.fold_left (fun accu i -> accu +. one_e i i Spin.Alfa) 0. mo_a)
-        +.
-        (List.fold_left (fun accu i -> accu +. one_e i i Spin.Beta) 0. mo_b)
-      in
-      let two =
-        let rec aux_same spin accu = function
-          | [] -> accu
-          | i :: rest -> 
-            let new_accu = 
-              List.fold_left (fun accu j -> accu +. two_e i j i j spin spin) accu rest
-            in
-            (aux_same [@tailcall]) spin new_accu rest
-        in
-        let rec aux_opposite accu other = function
-          | [] -> accu
-          | i :: rest ->
-            let new_accu =
-              List.fold_left (fun accu j -> accu +. two_e i j i j Spin.Alfa Spin.Beta) accu other
-            in
-            (aux_opposite [@tailcall]) new_accu other rest
-        in
-        (aux_same Spin.Alfa 0. mo_a)  +.  (aux_same Spin.Beta 0. mo_b)  +.
-        (aux_opposite 0. mo_a mo_b)
-      in
-      one +. two
-  in
-  let result =
-    match degree_a, degree_b with
-    | 1, 1 -> (* alpha-beta double *)
-      begin
-        let ha, pa, phase_a = Ex.single_of_spindet kia kja in
-        let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
-        let s1 = 
-          match phase_a, phase_b with
-          | Phase.Pos, Phase.Pos
-          | Phase.Neg, Phase.Neg -> fun _ two_e ->    two_e ha hb pa pb Spin.Alfa Spin.Beta
-          | Phase.Neg, Phase.Pos                       
-          | Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
-        in
-        let kk = Determinant.double_excitation Spin.Alfa ha pb Spin.Beta hb pa ki in
-        let kka = De.alfa kk and kkb = De.beta kk in
-        match Spindeterminant.(degree kia kka, degree kib kkb) with
-        | (1,1) ->
-          let s2 = 
-            begin
-              let ha, pa, phase_a = Ex.single_of_spindet kia kka in
-              let hb, pb, phase_b = Ex.single_of_spindet kib kkb in
-              match phase_a, phase_b with
-              | Phase.Pos, Phase.Pos
-              | Phase.Neg, Phase.Neg -> fun _ two_e ->    two_e ha hb pa pb Spin.Alfa Spin.Beta
-              | Phase.Neg, Phase.Pos                       
-              | Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
-            end
-          in fun x two_e -> 0.75 *. (s1 x two_e) +. 0.25 *. (s2 x two_e)  
-        | _ -> fun x two_e -> 0.
-      end
-
-    | 2, 0 -> (* alpha double *)
-      begin
-        let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in
-        match phase with 
-        | Phase.Pos -> fun _ two_e ->    two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
-        | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
-      end
-
-    | 0, 2 -> (* beta double *)
-      begin
-        let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in
-        match phase with 
-        | Phase.Pos -> fun _ two_e ->    two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
-        | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
-      end
-
-    | 1, 0    (* alpha single *)
-    | 0, 1    (* beta single *)
-           -> fun _ _ -> 0.
-    | 0, 0 -> (* diagonal element *)
-      diag_element
-
-    | _ -> assert false
-  in
-  List.map (fun (one_e, two_e) -> result one_e two_e) integrals
-
-
-let make integrals ki kj =
-  let degree_a, degree_b =
-    De.degrees ki kj 
-  in
-  if degree_a+degree_b > 2 then
-    List.map (fun _ -> 0.) integrals
-  else
-    non_zero integrals degree_a degree_b ki kj
-    
-
-
-
-let make_s2 ki kj =
-  let degree_a = De.degree_alfa ki kj in
-  let kia = De.alfa ki in
-  let kja = De.alfa kj in
-  if degree_a > 1 then 0.
-  else
-    let degree_b = De.degree_beta ki kj in
-    let kib = De.beta ki in
-    let kjb = De.beta kj in
-    match degree_a, degree_b with
-    | 1, 1 -> (* alpha-beta double *)
-      let ha, pa, phase_a = Ex.single_of_spindet kia kja in
-      let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
-      if ha = pb && hb = pa then
-        begin
-          match phase_a, phase_b with
-          | Phase.Pos, Phase.Pos
-          | Phase.Neg, Phase.Neg -> -1.
-          | Phase.Neg, Phase.Pos                       
-          | Phase.Pos, Phase.Neg -> 1.
-        end
-      else 0.
-    | 0, 0 -> 
-      let ba = Sp.bitstring kia and bb = Sp.bitstring kib in
-      let tmp = Bitstring.logxor ba bb in
-      let n_a = Bitstring.logand ba tmp |> Bitstring.popcount in
-      let n_b = Bitstring.logand bb tmp |> Bitstring.popcount in
-      let s_z = 0.5 *. float_of_int (n_a - n_b) in
-      float_of_int n_a +. s_z *. (s_z -. 1.)
-    | _ -> 0.
-
-
diff --git a/CI/F12CI.ml b/CI/F12CI.ml
index 5ee8586..2b241ca 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -4,7 +4,7 @@ module Ds = DeterminantSpace
 module De = Determinant
 module Sp = Spindeterminant
 
-type t = 
+type t =
   {
     mo_basis  : MOBasis.t ;
     det_space : DeterminantSpace.t ;
@@ -20,245 +20,6 @@ let mo_class t  =  Ds.mo_class @@ det_space t
 let eigensystem t =  Lazy.force t.eigensystem
 
 
-(*
-
-let single_matrices hf12_integrals density =
-    let nocc = Mat.dim1 density in
-    let nvir = Mat.dim2 density in
-    let { HF12.
-          simulation ; aux_basis ;
-          hf12 ; hf12_anti ;
-          hf12_single ; hf12_single_anti } = hf12_integrals
-    in
-    let f12 = MOBasis.f12_ints   aux_basis in
-    let eri = MOBasis.two_e_ints aux_basis in
-
-    let d = Mat.as_vec density in
-    let v_s = Mat.create nocc nvir in
-    let v_o = Mat.create nocc nvir in
-    let h_o, h_s, f_o, f_s = 
-      Mat.create nocc nvir ,
-      Mat.create nocc nvir ,
-      Mat.create nocc nvir ,
-      Mat.create nocc nvir ,
-    in
-    for a=1 to nvir do
-      for m=1 to occ do
-       for u=1 to nocc do
-        for t=1 to nocc do
-          let hmtau = ERI.get_phys eri m t a u
-          and hmtua = ERI.get_phys eri m t u a
-	  in
-          v_o.{t,u} <- hmtau;
-          v_s.{t,u} <- hmtau -. hmtua;
-        done
-       done;
-       h_o.{m,a} <- dot d_o @@ Mat.as_vec v_o;
-       h_s.{m,a} <- dot d_s @@ Mat.as_vec v_s
-       for u=1 to nocc do
-        for t=1 to nocc do
-          let fmtau = ERI.get_phys f12 m t a u
-          and fmtua = ERI.get_phys f12 m t u a
-	  in
-          v_o.{t,u} <- 0.375 *. fmtau +. 0.125 *. fmtua;
-          v_s.{t,u} <- 0.25 *, (fmtau -. fmtua);
-        done
-       done;
-       f_o.{m,a} <- dot d_o @@ Mat.as_vec v_o;
-       f_s.{m,a} <- dot d_s @@ Mat.as_vec v_s
-      done
-    done;
-*)
-  
-(*---
-let hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj =
-  let integrals = [
-
-    let { HF12.
-          simulation ; aux_basis ;
-          hf12 ; hf12_anti ;
-          hf12_single ; hf12_single_anti } = hf12_integrals
-    in
-
-    let kia = De.alfa ki and kib = De.beta ki in
-    let kja = De.alfa kj and kjb = De.beta kj in
-
-    let mo_a =
-      Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja)
-      |> Bitstring.to_list
-      |> Array.of_list
-    and mo_b =
-      Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb)
-      |> Bitstring.to_list
-      |> Array.of_list
-    in
-
-    let aux_mos = 
-      Util.list_range (MOBasis.size mo_basis) (MOBasis.size aux_basis)
-      |> Array.of_list
-    in
-
-    let one_e _ _ _ = 0. in
-
-    let h = MOBasis.ee_ints aux_basis in
-    let two_e_h i j k l s s' =
-      if s' <> s then
-        ERI.get_phys h i j k l
-      else
-        (ERI.get_phys h i j k l) -.  (ERI.get_phys h i j l k)
-    in
-
-    let f = MOBasis.f12_ints aux_basis in
-    let two_e_f i j k l s s' =
-      let fijkl = F12.get_phys f i j k l
-      and fijlk = F12.get_phys f i j l k
-      in
-      if s' <> s then
-        0.325 *. fijkl +. 0.125 *. fijlk
-      else
-        0.25 *. (fijkl -. fijlk)
-    in
-
-    (*
-    let mo_of_s = function
-    | Spin.Alfa -> mo_a
-    | Spin.Beta -> mo_b
-    in
-    *)
-
-    let two_e i j k l s s' =
-      (
-      if s = s' then
-	hf12_anti.{i,j,k,l} -. (
-	(Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +.
-	(Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) )
-      else
-	hf12.{i,j,k,l}  -. (
-	(Array.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +.
-	(Array.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) )
-      )
-      (*
-      +. Array.fold_left ( fun accu a -> accu +.
-         (Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s s) +.
-          Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s @@ Spin.other s) )
-         *. (two_e_f a j k l s s') ) 0. aux_mos 
-         *)
-    in
-
-    let three_e i j k l m n s s' s'' =
-         Array.fold_left (fun accu a -> accu +. two_e_h i j l a s   s'  *. two_e_f a k m n s'  s'') 0. aux_mos
-      -. Array.fold_left (fun accu a -> accu +. two_e_h j i m a s'  s   *. two_e_f a k l n s   s'') 0. aux_mos
-      +. Array.fold_left (fun accu a -> accu +. two_e_h j k m a s'  s'' *. two_e_f a i n l s'' s  ) 0. aux_mos
-      -. Array.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s'  *. two_e_f a i m l s'  s  ) 0. aux_mos
-      +. Array.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s   *. two_e_f a j l m s   s' ) 0. aux_mos
-      -. Array.fold_left (fun accu a -> accu +. two_e_h i k l a s   s'' *. two_e_f a j n m s'' s' ) 0. aux_mos
-    in
-    (one_e, two_e, Some three_e)
-  ]
-  in
-  CIMatrixElement.non_zero integrals deg_a deg_b ki kj
-  |> List.hd
-
-
-let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
-
-  if Parallel.master then
-    Printf.printf "Building matrix\n%!";
-  let det_space = 
-    ci.CI.det_space 
-  in
-  let mo_basis =
-    Ds.mo_basis det_space
-  in
-
-  let m_HF = 
-    
-    let f =
-      match Ds.determinants det_space with
-      | Ds.Arbitrary _ -> CI.create_matrix_arbitrary
-      | Ds.Spin      _ -> CI.create_matrix_spin_computed ~nmax:3
-    in
-    f (fun deg_a deg_b ki kj ->
-         hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj
-      ) det_space 
-
-  in
-
-  Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes)
-
---- *)
-
-
-let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
-  let integrals = [
-
-    let { HF12.
-          simulation ; aux_basis ;
-          hf12 ; hf12_anti ;
-          hf12_single ; hf12_single_anti } = hf12_integrals
-    in
-
-    let kia = De.alfa ki and kib = De.beta ki in
-    let 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 one_e _ _ _ = 0. in
-
-    let two_e i j k l s s' =
-      if s = s' then
-	hf12_anti.{i,j,k,l} -. (
-	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +.
-	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) 
-        )
-      else
-	hf12.{i,j,k,l}  -. (
-	(List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +.
-	(List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) 
-        )
-    in
-
-    let h = MOBasis.ee_ints aux_basis in
-    let two_e_h i j k l s s' =
-      if s' <> s then
-        ERI.get_phys h l k j i
-      else
-        (ERI.get_phys h l k j i) -.  (ERI.get_phys h k l j i)
-    in
-
-    let f = MOBasis.f12_ints aux_basis in
-    let two_e_f i j k l s s' =
-      if s' <> s then
-        F12.get_phys f i j k l
-      else
-        (F12.get_phys f i j k l) -.  (F12.get_phys f i j l k)
-    in
-
-    let mo_of_s = function
-    | Spin.Alfa -> mo_a
-    | Spin.Beta -> mo_b
-    in
-
-    let three_e i j k l m n s s' s'' =
-         List.fold_left (fun accu a -> accu +. two_e_h i j l a s   s'  *. two_e_f a k m n s'  s'') 0. (mo_of_s s' )
-      -. List.fold_left (fun accu a -> accu +. two_e_h j i m a s'  s   *. two_e_f a k l n s   s'') 0. (mo_of_s s  )
-      +. List.fold_left (fun accu a -> accu +. two_e_h j k m a s'  s'' *. two_e_f a i n l s'' s  ) 0. (mo_of_s s'')
-      -. List.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s'  *. two_e_f a i m l s'  s  ) 0. (mo_of_s s' )
-      +. List.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s   *. two_e_f a j l m s   s' ) 0. (mo_of_s s  )
-      -. List.fold_left (fun accu a -> accu +. two_e_h i k l a s   s'' *. two_e_f a j n m s'' s' ) 0. (mo_of_s s'')
-    in
-
-    (one_e, two_e, Some three_e)
-  ]
-  in
-  CIMatrixElement.non_zero integrals deg_a deg_b ki kj
-  |> List.hd
 
 
 let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
@@ -269,6 +30,11 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
     ci.CI.det_space
   in
 
+  let { HF12.
+        simulation ; aux_basis ;
+        f_0 ; f_1 ; f_2 ; f_3 } = hf12_integrals
+  in
+
   let m_HF =
 
     let f =
@@ -277,7 +43,12 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
       | Ds.Spin      _ -> CI.create_matrix_spin_computed ~nmax:3
     in
     f (fun deg_a deg_b ki kj ->
-         hf_ij_non_zero hf12_integrals deg_a deg_b ki kj
+         match deg_a + deg_b with
+         | 0 -> f_0 ki
+         | 1 -> f_1 ki kj
+         | 2 -> f_2 ki kj
+         | 3 -> f_3 ki kj
+         | _ -> assert false
       ) det_space
 
   in
@@ -285,10 +56,59 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
   Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes)
 
 
+let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l
+
+let array_3_init d1 d2 d3 f =
+  let result =
+      Bigarray.(Array3.create Float64 fortran_layout) d1 d2 d3
+  in
+  for k=1 to d3 do
+      for j=1 to d2 do
+          for i=1 to d1 do
+            result.{i,j,k} <- f i j k
+          done
+      done
+  done;
+  result
+
+let array_4_init d1 d2 d3 d4 f =
+  let result =
+      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4 |]
+  in
+  for l=1 to d4 do
+    for k=1 to d3 do
+      for j=1 to d2 do
+        for i=1 to d1 do
+           result.{i,j,k,l} <- f i j k l
+        done
+      done
+    done
+  done;
+  result
+
+let array_5_init d1 d2 d3 d4 d5 f =
+  let result =
+      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4;d5 |]
+  in
+  for m=1 to d5 do
+    for l=1 to d4 do
+      for k=1 to d3 do
+        for j=1 to d2 do
+          for i=1 to d1 do
+             result.{i,j,k,l,m} <- f i j k l m
+          done
+        done
+      done
+    done
+  done;
+  result
+
+
+
 let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () =
 
-  let det_space = 
-    DeterminantSpace.fci_of_mo_basis  mo_basis  ~frozen_core 
+  let det_space =
+    DeterminantSpace.fci_of_mo_basis  mo_basis  ~frozen_core
   in
 
   let ci = CI.make ~n_states:state det_space in
@@ -302,20 +122,21 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
     Parallel.broadcast (lazy x)
   in
 
+
   let eigensystem = lazy (
     let m_H =
       Lazy.force ci.CI.m_H
     in
 
-      
-    let rec iteration ~state psi = 
+
+    let rec iteration ~state psi =
 (*
 Format.printf "%a@." DeterminantSpace.pp_det_space @@ CI.det_space ci;
 Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi;
 *)
       let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
 
-      let delta = 
+      let delta =
          (* delta_i = {% $\sum_j c_j H_{ij}$ %} *)
          dressing_vector ~frozen_core hf12_integrals psi ci
          |> Matrix.to_mat
@@ -339,7 +160,7 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi;
       delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00;
 
 
-      let eigenvectors, eigenvalues = 
+      let eigenvectors, eigenvalues =
 
         let delta = lacpy delta in
         Mat.scal f delta;
@@ -357,8 +178,8 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi;
             )
         in
 
-        let matrix_prod c = 
-          let w = 
+        let matrix_prod c =
+          let w =
             Matrix.mm ~transa:`T m_H c
             |> Matrix.to_mat
           in
@@ -382,7 +203,7 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi;
           ))
 
       in
-      let eigenvectors = 
+      let eigenvectors =
         Conventions.rephase @@ Util.remove_epsilons eigenvectors
       in
 
@@ -396,11 +217,11 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi;
           (Mat.to_col_vecs eigenvectors).(0) )
       in
       if Parallel.master then
-        Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{state} 
+        Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{state}
           +. Simulation.nuclear_repulsion simulation);
 
 (*
-      let cabs_singles = 
+      let cabs_singles =
         let f =
           Fock.make_rhf ~density ~ao_basis:large_ao_basis
         in
diff --git a/MOBasis/HF12.ml b/MOBasis/HF12.ml
index 30794c3..2e9a6db 100644
--- a/MOBasis/HF12.ml
+++ b/MOBasis/HF12.ml
@@ -10,21 +10,61 @@ type q = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarra
 type t = {
   simulation : Simulation.t ;
   aux_basis  : MOBasis.t ;
-  hf12     : q ; (* hf12.{i,j,k,l} = \sum_{ab} <ij|H|ab> <ab|F12|kl>
-                                   + \sum_{am} <ij|H|am> <am|F12|kl>
-                                   + \sum_{bm} <ij|H|mb> <mb|F12|kl>
-                 *)
-  hf12_anti: q ; (* hf12_anti.{i,j,k,l} = 
-                      \sum_{ab} (<ij|H|ab> - <ij|H|ba>) (<ab|F12|kl> - <ab|F12|lk>)
-                    + \sum_{am} (<ij|H|am> - <ij|H|ma>) (<am|F12|kl> - <am|F12|lk>)
-                    + \sum_{bm} (<ij|H|mb> - <ij|H|bm>) (<mb|F12|kl> - <mb|F12|lk>)
-                    *)
-  hf12_single     : q ; (* hf12.{m,i,j,k,l} = \sum_{a} <ij|H|am> <am|F12|kl> *)
-  hf12_single_anti: q ; (* hf12_anti.{m,i,j,k,l} = 
-                           \sum_{ab} (<ij|H|am> - <ij|H|ma>) (<am|F12|kl> - <am|F12|lk>) *)
+  f_0 : Determinant.t -> float ;
+  f_1 : Determinant.t -> Determinant.t -> float ;
+  f_2 : Determinant.t -> Determinant.t -> float ;
+  f_3 : Determinant.t -> Determinant.t -> float ;
 }
 
 
+let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l
+
+let array_3_init d1 d2 d3 f =
+  let result =
+      Bigarray.(Array3.create Float64 fortran_layout) d1 d2 d3
+  in
+  for k=1 to d3 do
+      for j=1 to d2 do
+          for i=1 to d1 do
+            result.{i,j,k} <- f i j k
+          done
+      done
+  done;
+  result
+
+let array_4_init d1 d2 d3 d4 f =
+  let result =
+      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4 |]
+  in
+  for l=1 to d4 do
+    for k=1 to d3 do
+      for j=1 to d2 do
+        for i=1 to d1 do
+           result.{i,j,k,l} <- f i j k l
+        done
+      done
+    done
+  done;
+  result
+
+let array_5_init d1 d2 d3 d4 d5 f =
+  let result =
+      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4;d5 |]
+  in
+  for m=1 to d5 do
+    for l=1 to d4 do
+      for k=1 to d3 do
+        for j=1 to d2 do
+          for i=1 to d1 do
+             result.{i,j,k,l,m} <- f i j k l m
+          done
+        done
+      done
+    done
+  done;
+  result
+
+
 let make ~simulation ~mo_basis ~aux_basis_filename () =
   
   let f12 = Util.of_some @@ Simulation.f12 simulation in                                      
@@ -53,209 +93,768 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   let aux_num = MOBasis.size aux_basis in
 
   (* Fire calculation of F12 and ERI *)
-  let f12 = 
-    MOBasis.f12_ints aux_basis
-  in
-  let eri =
-    MOBasis.two_e_ints aux_basis
-  in
+  ignore @@  MOBasis.f12_ints aux_basis ;
+  ignore @@  MOBasis.two_e_ints aux_basis ;
 
   (* Compute the <ij|QHF|kl> integrals *)
   if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!";
 
-  let hf12, hf12_anti, hf12_single, hf12_single_anti =
-    let create4 n =
-      Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n |] 
-    in
-    let create5 n =
-      Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n ; n |] 
-    in
-    create4 mo_num ,
-    create4 mo_num ,
-    create5 mo_num ,
-    create5 mo_num 
+  let mos_cabs =
+    Util.list_range (mo_num+1) aux_num
   in
 
-
-  let h_s  = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num aux_num aux_num in
-  let f_s  = Bigarray.Array3.create Float64 Bigarray.fortran_layout aux_num aux_num mo_num in
-  let h_o  = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num aux_num aux_num in
-  let f_o  = Bigarray.Array3.create Float64 Bigarray.fortran_layout aux_num aux_num mo_num in
-  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
-        h_s.{i, a, b} <- 0. ;
-        h_o.{i, a, b} <- 0. 
-     done
-    done
-  done;
-
-  for k=1 to mo_num do
-    for b=1 to mo_num do
-      for a=1 to mo_num do
-        f_s.{a, b, k} <- 0. ;
-        f_o.{a, b, k} <- 0.
-     done
-    done
-  done;
-
-  let task (j,l) = 
-    
-      let h i a b = 
-         let ijab = ERI.get_phys eri i j a b
-         and ijba = ERI.get_phys eri i j b a
-         in
-         h_s.{i, a, b} <- ijab -. ijba ;
-         h_o.{i, a, b} <- ijab
-      and f a b k =
-         let abkl = F12.get_phys f12 a b k l
-         and ablk = F12.get_phys f12 a b l k
-         in
-         f_s.{a, b, k} <- 0.25 *. (abkl -. ablk) ;
-         f_o.{a, b, k} <- 0.375 *. abkl +. 0.125 *. ablk
-      in
-
-      for a=mo_num+1 to aux_num do
-        for b=mo_num+1 to aux_num do
-         for i=1 to mo_num do
-            h i a b
-         done
-        done
-      done;
-
-      for k=1 to mo_num do
-        for b=mo_num+1 to aux_num do
-          for a=mo_num+1 to aux_num do
-            f a b k
-         done
-        done
-      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 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' =
-        Bigarray.(reshape (genarray_of_array3 h_o)) [| mo_num ; aux_num*aux_num |]
-        |> Bigarray.array2_of_genarray
-      in
-      let f_o' =
-        Bigarray.(reshape (genarray_of_array3 f_o)) [| aux_num*aux_num ; mo_num |] 
-        |> Bigarray.array2_of_genarray
-      in
-      let h_s' =
-        Bigarray.(reshape (genarray_of_array3 h_s)) [| mo_num ; aux_num*aux_num |]
-        |> Bigarray.array2_of_genarray
-      in
-      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 ~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
+  let mos_in   =
+    Util.list_range 1 mo_num
   in
 
-  let tasks = 
-    let rec next accu = function
-      | _, 0 -> accu
-      | 0, l -> next accu  (mo_num, l-1)
-      | j, l -> next ((j,l) :: accu)  ((j-1), l)
+  let mos_a k =
+    Determinant.alfa k
+    |> Spindeterminant.to_list
+  in
+
+  let mos_b k =
+    Determinant.beta k
+    |> Spindeterminant.to_list
+  in
+
+  let h_one =
+    let h =
+      MOBasis.one_e_ints aux_basis
+    in fun i j _ -> h.{i,j}
+  in
+
+  let h_two =
+    let two_e_ints = MOBasis.two_e_ints aux_basis in
+    let h2 i j k l (s:Spin.t) (s':Spin.t) =
+      if s' <> s then
+	  ERI.get_phys two_e_ints i j k l
+      else
+	  (ERI.get_phys two_e_ints i j k l) -.
+	  (ERI.get_phys two_e_ints i j l k)
     in
-    next [] (mo_num, mo_num)
-    |> Stream.of_list
+    h2
   in
 
+  let f_two =
+    let f12_ints = MOBasis.f12_ints aux_basis in
+    let f2 i j k l (s:Spin.t) (s':Spin.t) =
+      if s' <> s then
+	  0.375 *. F12.get_phys f12_ints i j k l +.
+	  0.125 *. F12.get_phys f12_ints i j l k
+      else
+	  0.25 *. (
+	  (F12.get_phys f12_ints i j k l) -.
+	  (F12.get_phys f12_ints i j l k) )
+    in
+    f2
+  in
 
-  Farm.run ~f:task ~ordered:true tasks
-  |> 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
-        hf12.{i,j,k,l}      <- hf_o.{i,k} ;
-        hf12_anti.{i,j,k,l} <- hf_s.{i,k} ;
-        for m=1 to mo_num do
-          hf12_single.{m,i,j,k,l}      <- hfm_o.{m,i,k} ;
-          hf12_single_anti.{m,i,j,k,l} <- hfm_s.{m,i,k}
-        done
-      done
-    done );
+  let f_one = fun _ _ _ -> 0. in
 
 
 
-  (*
-  for l=1 to mo_num do
-    for k=1 to mo_num do
-      for j=1 to mo_num do
-        for i=1 to mo_num do
-          Printf.printf "%d %d %d %d %e\n" i j k l result.{i,j,k,l} 
-        done
-      done
-    done
-  done;
-  Printf.printf "%!";
-  *)
+  (* Pre-compute dressed integrals *)
+  let m_0111_1H_1F =
+    Vec.init mo_num (fun i ->
+	  sum mos_cabs (fun a ->
+	      h_one i a Spin.Alfa *. f_one a i Spin.Alfa ))
+  in
+
+  let m_0111_1H_2Fa, m_0111_2Ha_2Fa =
+
+      let m_0122_Haa =
+	array_3_init mo_num mo_num mo_num (fun i j k ->
+	    sum mos_cabs (fun a ->
+	      h_two i j a k Spin.Alfa Spin.Alfa *. f_two a k i j Spin.Alfa Spin.Alfa
+	      ) )
+      in
+
+      let m_0111_1H_2Fa =
+	Mat.init_cols mo_num mo_num (fun i j ->
+	    sum mos_cabs (fun a ->
+	      h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Alfa +.
+	      h_two i j a j Spin.Alfa Spin.Alfa *. f_one a i Spin.Alfa
+	    ) +.
+	    if i < j then 0. else
+		begin
+		    sum mos_cabs (fun a ->
+			sum mos_cabs (fun b -> if b >= a then 0. else
+			  h_two i j a b Spin.Alfa Spin.Alfa *. f_two a b i j Spin.Alfa Spin.Alfa
+			)
+		    ) +.
+		    sum mos_in (fun k -> m_0122_Haa.{i,j,k})
+		end
+	)
+      in
+
+      let m_0111_2Ha_2Fa =
+	array_3_init mo_num mo_num mo_num (fun i j k ->
+	  sum mos_cabs (fun a ->
+	    h_two i j a j Spin.Alfa Spin.Alfa *.
+	    f_two a k i k Spin.Alfa Spin.Alfa
+	  ) -. if i < j then 0. else m_0122_Haa.{i,j,k}
+	)
+      in m_0111_1H_2Fa, m_0111_2Ha_2Fa
+  in
+
+  let m_0111_1H_2Fb, m_0111_2Hb_2Fb =
+      let m_0122_Hab =
+	array_3_init mo_num mo_num mo_num (fun i j k ->
+	    sum mos_cabs (fun a ->
+	      h_two i j a k Spin.Alfa Spin.Beta *. f_two a k i j Spin.Alfa Spin.Beta
+	    ) )
+    in
+
+    let m_0111_1H_2Fb =
+      Mat.init_cols mo_num mo_num (fun i j ->
+          sum mos_cabs (fun a ->
+            h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Beta +.
+            h_two i j a j Spin.Alfa Spin.Beta *. f_one a i Spin.Alfa +.
+            h_one j a Spin.Alfa *. f_two a i j i Spin.Alfa Spin.Beta +.
+            h_two j i a i Spin.Alfa Spin.Beta *. f_one a j Spin.Alfa
+          ) +.
+        sum mos_in (fun k -> m_0122_Hab.{i,j,k} +. m_0122_Hab.{j,i,k} ) +.
+        sum mos_cabs (fun a ->
+            sum mos_cabs (fun b ->
+            h_two i j a b Spin.Alfa Spin.Beta *. f_two a b i j Spin.Alfa Spin.Beta
+            )
+          )
+        )
+    in
+
+    let m_0111_2Hb_2Fb =
+      array_3_init mo_num mo_num mo_num (fun i j k ->
+        sum mos_cabs (fun a ->
+          h_two k i a i Spin.Alfa Spin.Beta *.
+          f_two a j k j Spin.Alfa Spin.Beta +.
+          h_two i k a k Spin.Alfa Spin.Beta *.
+          f_two a j i j Spin.Alfa Spin.Alfa
+        ) -. m_0122_Hab.{k,i,j}
+      )
+    in
+    m_0111_1H_2Fb, m_0111_2Hb_2Fb
+  in
+
+  let m_0111_2Ha_2Fb =
+    array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+	h_two i j a j Spin.Alfa Spin.Alfa *.
+	f_two a k i k Spin.Alfa Spin.Beta
+      )
+    )
+  in
+
+  let f_0 ki =
+
+    let mos_i, mos_i' = mos_a ki, mos_b ki in
+    let same = (mos_i = mos_i') in
+
+    (* Alpha *)
+    let a =
+	sum mos_i  (fun i -> m_0111_1H_1F.{i})
+    in
+
+    let b =
+	if same then a else
+	sum mos_i' (fun i -> m_0111_1H_1F.{i})
+    in
+
+    let aa =
+	sum mos_i (fun j ->
+	  sum mos_i (fun i -> m_0111_1H_2Fa.{i,j} ))
+    in
+
+    let bb =
+      if same then aa else
+      sum mos_i' (fun j ->
+	sum mos_i' (fun i -> m_0111_1H_2Fa.{i,j} ))
+    in
+
+    let ab =
+      sum mos_i' (fun j ->
+	sum mos_i (fun i -> m_0111_1H_2Fb.{i,j}  ))
+    in
+
+    let aaa =
+	sum mos_i (fun k ->
+	  sum mos_i (fun j ->
+	    sum mos_i (fun i -> m_0111_2Ha_2Fa.{i,j,k} )))
+    in
+
+    let bbb =
+      if same then aaa else
+	sum mos_i' (fun k ->
+	  sum mos_i' (fun j ->
+	    sum mos_i' (fun i -> m_0111_2Ha_2Fa.{i,j,k} )))
+    in
+
+    let baa =
+	sum mos_i' (fun k ->
+	  sum mos_i (fun j ->
+	    sum mos_i (fun i ->
+	      m_0111_2Ha_2Fb.{i,j,k} +. m_0111_2Hb_2Fb.{i,j,k}
+	)))
+    in
+
+    let bba =
+	sum mos_i  (fun k ->
+	  sum mos_i' (fun j ->
+	    sum mos_i' (fun i ->
+	      m_0111_2Ha_2Fb.{i,j,k} +. m_0111_2Hb_2Fb.{j,i,k}
+	)))
+    in
+
+    a +. b +. aa +. bb +. ab +. aaa +. baa +. bba +. bbb
+  in
+
+  let m_1111_1H_1F =
+    Mat.init_cols mo_num mo_num (fun i k ->
+	  sum mos_cabs (fun a -> h_one i a Spin.Alfa *. f_one a k Spin.Alfa ))
+  in
+
+  let m_1111_2Ha_2Fa =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      if l=i then
+	  sum mos_cabs (fun a ->
+	    h_two j l a l Spin.Alfa Spin.Alfa *.
+	    f_two a i j k Spin.Alfa Spin.Alfa
+	  )
+      else
+	  sum mos_cabs (fun a ->
+	    h_two i j a j Spin.Alfa Spin.Alfa *.
+	    f_two a l k l Spin.Alfa Spin.Alfa +.
+	    h_two j l a l Spin.Alfa Spin.Alfa *.
+	    f_two a i j k Spin.Alfa Spin.Alfa )
+    )
+  in
+
+  let m_1111_2Hb_2Fa =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      if l=i then
+	  sum mos_cabs (fun a ->
+	    h_two j l a l Spin.Alfa Spin.Beta *.
+	    f_two a i j k Spin.Alfa Spin.Beta
+	  )
+      else
+	  sum mos_cabs (fun a ->
+	    h_two i j a j Spin.Alfa Spin.Beta *.
+	    f_two a l k l Spin.Alfa Spin.Alfa +.
+	    h_two j l a l Spin.Alfa Spin.Beta *.
+	    f_two a i j k Spin.Alfa Spin.Beta
+	  )
+    )
+  in
+
+  let m_1111_2Ha_2Fb =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j a j Spin.Alfa Spin.Alfa *.
+	f_two a l k l Spin.Alfa Spin.Beta +.
+	h_two j l a l Spin.Alfa Spin.Beta *.
+	f_two a i j k Spin.Alfa Spin.Alfa
+      )
+    )
+  in
+
+  let m_1111_2Hb_2Fb =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j a j Spin.Alfa Spin.Beta *.
+	f_two a l k l Spin.Alfa Spin.Beta +.
+	h_two j l a l Spin.Alfa Spin.Alfa *.
+	f_two a i j k Spin.Alfa Spin.Beta
+      )
+    )
+  in
+
+  let m_1121_2Ha_2Fa =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j k a Spin.Alfa Spin.Alfa *.
+	f_two l a l j Spin.Alfa Spin.Alfa
+      )
+    )
+  in
+
+  let m_1121_2Hb_2Fa =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j k a Spin.Alfa Spin.Beta *.
+	f_two l a l j Spin.Alfa Spin.Alfa
+      )
+    )
+  in
+
+  let m_1121_2Ha_2Fb =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j k a Spin.Alfa Spin.Alfa *.
+	f_two l a l j Spin.Alfa Spin.Beta
+      )
+    )
+  in
+
+  let m_1121_2Hb_2Fb =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j k a Spin.Alfa Spin.Beta *.
+	f_two l a l j Spin.Alfa Spin.Beta
+      )
+    )
+  in
+
+  let m_1122_va =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two j i a l Spin.Alfa Spin.Alfa *.
+	f_two l a k j Spin.Alfa Spin.Alfa
+      )
+    )
+  in
+
+  let m_1111_1H_2Fa =
+    array_3_init mo_num mo_num mo_num (fun j i k ->
+	sum mos_in (fun l -> m_1122_va.{l,j,i,k}) +.
+	sum mos_cabs (fun a ->
+	  h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Alfa +.
+	  h_two i j a j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.
+	  h_one j a Spin.Alfa *. f_two a i j k Spin.Alfa Spin.Alfa +.
+	  h_two i j k a Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa +.
+	  sum mos_cabs (fun b -> if b > a then 0. else
+	    h_two i j a b Spin.Alfa Spin.Alfa *.
+	    f_two a b k j Spin.Alfa Spin.Alfa )
+	)
+      )
+  in
+
+  let m_1122_v2 =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j l a Spin.Alfa Spin.Beta *.
+	f_two l a k j Spin.Alfa Spin.Beta
+      )
+    )
+  in
+
+  let m_1122_v3 =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two i j a l Spin.Alfa Spin.Beta *.
+	f_two a l k j Spin.Alfa Spin.Beta
+      )
+    )
+  in
+
+  let m_1111_1H_2Fb =
+    array_3_init mo_num mo_num mo_num (fun j i k ->
+	sum mos_in (fun l -> m_1122_v2.{l,j,i,k} +. m_1122_v3.{l,j,i,k}) +.
+	sum mos_cabs (fun a ->
+	  h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Beta +.
+	  h_two i j a j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.
+	  h_one j a Spin.Beta *. f_two a i j k Spin.Alfa Spin.Beta +.
+	  h_two i j k a Spin.Alfa Spin.Beta *. f_one a j Spin.Beta +.
+	  sum mos_cabs (fun b ->
+	    h_two i j a b Spin.Alfa Spin.Beta *.
+	    f_two a b k j Spin.Alfa Spin.Beta
+	  )
+	)
+      )
+  in
+
+  let m_1122_oa =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+    if l > j then
+      sum mos_cabs (fun a ->
+	h_two j l a k Spin.Alfa Spin.Alfa *.
+	f_two i a l j Spin.Alfa Spin.Alfa
+      )
+    else 0.
+    )
+  in
+
+  let m_1122_o =
+    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
+      sum mos_cabs (fun a ->
+	h_two l j k a Spin.Alfa Spin.Beta *.
+	f_two i a l j Spin.Alfa Spin.Beta
+      )
+    )
+  in
+
+  let f_1 ki kj =
+      let i, k, s, phase =
+	    match Excitation.of_det ki kj with
+	    | Excitation.(Single (phase, { hole ; particle ; spin })) ->
+	      hole, particle, spin, phase
+	    | _ -> assert false
+      in
+
+      let mos_novirt, mos_novirt' =
+	let alfa =
+	  let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
+	  Bitstring.to_list (Bitstring.logor i j)
+	in
+	let beta  =
+	  let i = Spindeterminant.bitstring @@ Determinant.beta ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.beta kj in
+	  Bitstring.to_list (Bitstring.logor i j)
+	in
+	match s with
+	| Spin.Alfa -> alfa, beta
+	| Spin.Beta -> beta, alfa
+      in
+
+      let mos_ij, mos_ij' =
+	let alfa =
+	  let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
+	  Bitstring.to_list (Bitstring.logand i j)
+	in
+	let beta  =
+	  let i = Spindeterminant.bitstring @@ Determinant.beta ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.beta kj in
+	  Bitstring.to_list (Bitstring.logand i j)
+	in
+	match s with
+	| Spin.Alfa -> alfa, beta
+	| Spin.Beta -> beta, alfa
+      in
+
+      let mos_i, mos_i' =
+	match s with
+	| Spin.Alfa -> mos_a ki, mos_b ki
+	| Spin.Beta -> mos_b ki, mos_a ki
+      in
+
+      let mos_j, mos_j' =
+	match s with
+	| Spin.Alfa -> mos_a kj, mos_b kj
+	| Spin.Beta -> mos_b kj, mos_a kj
+      in
+
+      let result =
+      m_1111_1H_1F.{i,k} +.
+      sum mos_ij (fun j ->
+	    m_1111_1H_2Fa.{j,i,k}
+	+. sum mos_i  (fun l -> m_1111_2Ha_2Fa.{l,j,i,k})
+	+. sum mos_i' (fun l -> m_1111_2Ha_2Fb.{l,j,i,k})
+	+. sum mos_j  (fun l -> m_1121_2Ha_2Fa.{l,j,i,k})
+	+. sum mos_j' (fun l -> m_1121_2Ha_2Fb.{l,j,i,k})
+	-. sum mos_novirt (fun l -> m_1122_va.{l,j,i,k})
+	-. sum mos_ij (fun l -> m_1122_oa.{l,j,i,k} )
+	) +.
+      sum mos_ij' (fun j ->
+	    m_1111_1H_2Fb.{j,i,k}
+	+. sum mos_i  (fun l -> m_1111_2Hb_2Fa.{l,j,i,k})
+	+. sum mos_i' (fun l -> m_1111_2Hb_2Fb.{l,j,i,k})
+	+. sum mos_j  (fun l -> m_1121_2Hb_2Fb.{l,j,i,k})
+	+. sum mos_j' (fun l -> m_1121_2Hb_2Fa.{l,j,i,k})
+	-. sum mos_novirt  (fun l ->  m_1122_v2.{l,j,i,k})
+	-. sum mos_novirt' (fun l ->  m_1122_v3.{l,j,i,k})
+	-. sum mos_ij    (fun l ->  m_1122_o.{l,j,i,k})
+	)
+      in
+      match phase with
+      | Phase.Pos -> result
+      | Phase.Neg -> -. result
+  in
+
+  let m_2112_1H_2Fa =
+    array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+      sum mos_cabs (fun a ->
+	  h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.
+	  h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa +.
+	  h_two i j a l Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.
+	  h_two i j k a Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa +.
+	  sum mos_in  (fun m -> -. h_two i j a m Spin.Alfa Spin.Alfa *.
+				  f_two m a k l Spin.Alfa Spin.Alfa) +.
+	  sum mos_cabs (fun b -> if b >= a then 0. else
+	      h_two i j a b Spin.Alfa Spin.Alfa *. f_two a b k l Spin.Alfa Spin.Alfa
+	  )
+	)
+      )
+  in
+
+  let m_2112_1H_2Fb =
+    array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+      sum mos_cabs (fun a ->
+	  h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.
+	  h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Beta +.
+	  h_two i j a l Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.
+	  h_two i j k a Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa +.
+	  sum mos_in (fun m ->
+	    h_two i j a m Spin.Alfa Spin.Beta *.  f_two a m k l Spin.Alfa Spin.Beta +.
+	    h_two i j m a Spin.Alfa Spin.Beta *.  f_two m a k l Spin.Alfa Spin.Beta ) +.
+	  sum mos_cabs (fun b ->
+	      h_two i j a b Spin.Alfa Spin.Beta *. f_two a b k l Spin.Alfa Spin.Beta
+	  )
+	)
+      )
+  in
+
+  let m_2112_2Ha_2Fa =
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	  h_two i n a n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.
+	  h_two j n a n Spin.Alfa Spin.Alfa *. f_two a i l k Spin.Alfa Spin.Alfa
+	)
+      )
+  in
+
+  let m_2112_2Hb_2Fa =
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	  h_two i n a n Spin.Alfa Spin.Beta *. f_two a j k l Spin.Alfa Spin.Alfa +.
+	  h_two j n a n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Alfa
+	)
+      )
+  in
+
+  let m_2112_2Ha_2Fb =
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	  h_two i n a n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.
+	  h_two j n a n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Beta )
+      )
+  in
+
+  let m_2121_2Ha_2Fa =
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	  h_two i j a l Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Alfa +.
+	  h_two j i a k Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Alfa
+	)
+      )
+  in
+
+  let m_2121_2Hb_2Fa =
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	  h_two i j a l Spin.Alfa Spin.Beta *. f_two a n k n Spin.Alfa Spin.Alfa +.
+	  h_two j i a k Spin.Alfa Spin.Beta *. f_two a n l n Spin.Alfa Spin.Beta
+	)
+      )
+  in
+
+  let m_2121_2Ha_2Fb =
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	  h_two i j a l Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Beta +.
+	  h_two j i a k Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Beta
+	)
+      )
+  in
+
+  let m_2122_2Ha_2Fa_ij =
+    let s = Spin.Alfa in
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	    h_two i n a k s s  *.  f_two j a n l s s
+	  +. h_two i n a l s s  *.  f_two j a k n s s
+	  -. h_two j n a k s s  *.  f_two i a n l s s
+	  -. h_two j n a l s s  *.  f_two i a k n s s
+	)
+      )
+  in
+
+  let m_2122_2Hb_2Fb_ij =
+    let s, s' = Spin.(Alfa, Beta) in
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	    h_two n i a k s  s  *.  f_two a j n l s  s'
+	+. h_two n j a l s  s' *.  f_two i a k n s  s
+	-. h_two n j k a s  s' *.  f_two i a n l s  s'
+      )
+    )
+  in
+
+  let m_2122_2Hb_2Fb_ij2 =
+    let s, s' = Spin.(Alfa, Beta) in
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	-. h_two i n a l s  s' *.  f_two a j k n s  s' +.
+	(if n < j then
+		    h_two i n k a s  s' *.  f_two j a l n s' s'
+		+. h_two n j a l s' s' *.  f_two i a k n s  s'
+	else
+		-. h_two i n k a s  s' *.  f_two j a n l s' s'
+		-. h_two j n a l s' s' *.  f_two i a k n s  s'
+	)
+      )
+    )
+  in
+
+  let m_2122_2Ha_2Fa_ij2 =
+    let s, s' = Spin.(Alfa, Beta) in
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a ->
+	  h_two i n k a s s' *.  f_two j a l n s s'
+	+. h_two j n l a s s' *.  f_two i a k n s s'
+	-. h_two i n l a s s' *.  f_two j a k n s s'
+	-. h_two j n k a s s' *.  f_two i a l n s s'
+	)
+      )
+  in
+
+  let m_2122_2Ha_2Fa_nv =
+    let s = Spin.Alfa in
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a -> h_two i j a n s  s  *.  f_two n a k l s  s ) )
+  in
+
+  let m_2122_2Hb_2Fb_nv =
+    let s, s' = Spin.(Alfa, Beta) in
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a -> -. h_two i j a n s s'  *.  f_two a n k l s s' ) )
+  in
+
+  let m_2122_2Hb_2Fb_nv2 =
+    let s, s' = Spin.(Alfa, Beta) in
+    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
+      sum mos_cabs (fun a -> -. h_two i j n a s s'  *.  f_two n a k l s s' ) )
+  in
+
+  let f_2 ki kj =
+      let i, j, k, l, s, s', phase =
+	    match Excitation.of_det ki kj with
+	    | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
+	      hole, hole', particle, particle', spin, spin', phase
+	    | _ -> assert false
+      in
+      let mos_i, mos_i' =
+	match s with
+	| Spin.Alfa -> mos_a ki, mos_b ki
+	| Spin.Beta -> mos_b ki, mos_a ki
+      in
+
+      let mos_j, mos_j' =
+	match s with
+	| Spin.Alfa -> mos_a kj, mos_b kj
+	| Spin.Beta -> mos_b kj, mos_a kj
+      in
+
+      let mos_ij, mos_ij' =
+	let alfa =
+	  let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
+	  Bitstring.to_list (Bitstring.logand i j)
+	in
+	let beta  =
+	  let i = Spindeterminant.bitstring @@ Determinant.beta ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.beta kj in
+	  Bitstring.to_list (Bitstring.logand i j)
+	in
+	match s with
+	| Spin.Alfa -> alfa, beta
+	| Spin.Beta -> beta, alfa
+      in
+
+      let mos_novirt, mos_novirt' =
+	let alfa =
+	  let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
+	  Bitstring.to_list (Bitstring.logor i j)
+	in
+	let beta  =
+	  let i = Spindeterminant.bitstring @@ Determinant.beta ki in
+	  let j = Spindeterminant.bitstring @@ Determinant.beta kj in
+	  Bitstring.to_list (Bitstring.logor i j)
+	in
+	match s with
+	| Spin.Alfa -> alfa, beta
+	| Spin.Beta -> beta, alfa
+      in
+
+      let result =
+	if s = s' then
+	  m_2112_1H_2Fa.{i,j,k,l} +.
+	  sum mos_i      (fun n -> m_2112_2Ha_2Fa.{n,i,j,k,l} ) +.
+	  sum mos_i'     (fun n -> m_2112_2Hb_2Fa.{n,i,j,k,l} ) +.
+	  sum mos_j      (fun n -> m_2121_2Ha_2Fa.{n,i,j,k,l} ) +.
+	  sum mos_j'     (fun n -> m_2121_2Ha_2Fb.{n,i,j,k,l} ) +.
+	  sum mos_ij     (fun n -> m_2122_2Ha_2Fa_ij.{n,i,j,k,l} ) +.
+	  sum mos_ij'    (fun n -> m_2122_2Ha_2Fa_ij2.{n,i,j,k,l} ) +.
+	  sum mos_novirt (fun n -> m_2122_2Ha_2Fa_nv.{n,i,j,k,l} )
+	else
+	  m_2112_1H_2Fb.{i,j,k,l} +.
+	  sum mos_i      (fun n -> m_2112_2Ha_2Fb.{n,i,j,k,l} ) +.
+	  sum mos_i'     (fun n -> m_2112_2Ha_2Fb.{n,j,i,l,k} ) +.
+	  sum mos_j      (fun n -> m_2121_2Hb_2Fa.{n,i,j,k,l} ) +.
+	  sum mos_j'     (fun n -> m_2121_2Hb_2Fa.{n,j,i,l,k} ) +.
+	  sum mos_novirt'(fun n -> m_2122_2Hb_2Fb_nv.{n,i,j,k,l} ) +.
+	  sum mos_novirt (fun n -> m_2122_2Hb_2Fb_nv2.{n,i,j,k,l} )+.
+	  sum mos_ij     (fun n -> m_2122_2Hb_2Fb_ij.{n,i,j,k,l} ) +.
+	  sum mos_ij'    (fun n -> m_2122_2Hb_2Fb_ij2.{n,i,j,k,l} )
+      in
+
+      match phase with
+      | Phase.Pos -> result
+      | Phase.Neg -> -. result
+
+  in
+
+  let f_3 ki kj =
+      let i, j, m, k, l, n, s1, s2, s3, phase =
+	    match Excitation.of_det ki kj with
+	    | Excitation.(Triple (phase,
+		{ hole=h1 ; particle=p1 ; spin=s1 },
+		{ hole=h2 ; particle=p2 ; spin=s2 },
+		{ hole=h3 ; particle=p3 ; spin=s3 }) ) ->
+		      h1, h2, h3, p1, p2, p3, s1, s2, s3, phase
+	    | _ -> assert false
+      in
+
+      let result =
+	let open Spin in
+	match s1, s2, s3 with
+	| Alfa, Alfa, Alfa
+	| Beta, Beta, Beta ->
+	    sum mos_cabs (fun a ->
+		  h_two i j a k s1 s2 *.  f_two m a l n s3 s3
+	      +. h_two i j a n s1 s2 *.  f_two m a k l s3 s2
+	      +. h_two i m a l s1 s3 *.  f_two j a k n s2 s3
+	      +. h_two j m a k s2 s3 *.  f_two i a l n s1 s3
+	      +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2
+	      -. h_two i j a l s1 s2 *.  f_two m a k n s3 s3
+	      -. h_two i m a k s1 s3 *.  f_two j a l n s2 s3
+	      -. h_two i m a n s1 s3 *.  f_two j a k l s2 s2
+	      -. h_two j m a l s2 s3 *.  f_two i a k n s1 s3 )
+	| Alfa, Alfa, Beta
+	| Beta, Beta, Alfa ->
+	    sum mos_cabs (fun a ->
+		  h_two i j a l s1 s2 *.  f_two a m k n s1 s3
+	      +. h_two i m k a s1 s3 *.  f_two j a l n s2 s3
+	      +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2
+	      +. h_two j m l a s2 s3 *.  f_two i a k n s1 s3
+	      -. h_two i j a k s1 s2 *.  f_two a m l n s1 s3
+	      -. h_two i m a n s1 s3 *.  f_two j a k l s2 s2
+	      -. h_two i m l a s1 s3 *.  f_two j a k n s2 s3
+	      -. h_two j m k a s2 s3 *.  f_two i a l n s1 s3
+	    )
+	| Alfa, Beta, Beta
+	| Beta, Alfa, Alfa ->
+	    sum mos_cabs (fun a ->
+                 h_two i j a l s1 s2 *.  f_two a m k n s1 s3
+	      +. h_two i m a n s1 s3 *.  f_two a j k l s1 s2
+	      +. h_two i m k a s1 s3 *.  f_two j a l n s2 s3
+	      +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2
+	      -. h_two i j a n s1 s2 *.  f_two a m k l s1 s2
+	      -. h_two i j k a s1 s2 *.  f_two m a l n s2 s3
+	      -. h_two i m a l s1 s3 *.  f_two a j k n s1 s3
+	      -. h_two j m a l s2 s3 *.  f_two i a k n s1 s3
+	    )
+	| Beta, Alfa, Beta
+	| Alfa, Beta, Alfa ->  assert false (*TODO *)
+      in
+      match phase with
+      | Phase.Pos -> result
+      | Phase.Neg -> -. result
+  in
 
   let result = 
     { simulation ; aux_basis ;
-      hf12 ; hf12_anti ;
-      hf12_single ; hf12_single_anti
+      f_0 ; f_1 ; f_2 ; f_3 
     } 
   in
   Parallel.broadcast (lazy result)