diff --git a/Basis/F12.ml b/Basis/F12.ml
index 2540ea8..3f0c79f 100644
--- a/Basis/F12.ml
+++ b/Basis/F12.ml
@@ -373,9 +373,6 @@ let of_basis_parallel f12 basis =
       Fis.create ~size:n `Dense
   in
 
-  (*
-  let lambda_inv = -. 1. /. f12.expo_s in
-  *)
   Farm.run ~ordered:false ~f input_stream
   |> Stream.iter (fun l ->
       Array.iter (fun (i_c,j_c,k_c,l_c,value) ->
diff --git a/CI/CI.ml b/CI/CI.ml
index c192bec..2f19a68 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -403,25 +403,34 @@ let create_matrix_spin_computed f det_space =
     let h123 = ref (fun _ -> 0.)  in
 
     let g i =
+      (*
+        i      : index of the i-th determinant. 1 <= i <= ndet
+        i_prev : index of the i-th determinant in the previous function call.
+                  1 <= i_prev <= ndet
+        i_a    : index of the i_a-th alpha determinant. 0 <= i_a < n_alfa
+        i_b    : index of the i_b-th beta determinant. 0 <= i_b < n_beta
+        j0     : index - 1 of the first determinant with the same alfa component
+                  0 <= j0 < n_beta*(n_alfa-1)
+        j1     : index - 1 of the next determinant with the 1st beta determinant
+                  n_beta <= j1 <= ndet
+        j_a    : index of the j_a-th alpha determinant. 0 <= j_a < n_alfa
+        j_b    : index of the j_b-th beta determinant. 0 <= j_b < n_beta
+      *)
       if i <> !i_prev then
         begin
           i_prev := i;
           let i_a = (i-1)/n_beta in
           let h1 = h i_a in
           let i_b = i - i_a*n_beta - 1  in
-          let j0 = ref max_int in
-          let j1 = ref min_int in
+          let j0 = ref (2*ndet) in
+          let j1 = ref (2*ndet) in
           let j_a = ref 0 in
           result := fun j ->
-              (* The following test checks if !j0 < j < !j1 *)
-              if  (!j1 - j) lor (j - !j0) > 0 then
-                begin
-                  let j_b = j - !j0 - 1 in
-                  !h123 j_b
-                end
+              if  (!j0 < j) && (j <= !j1) then 
+                  ()
               else
                 begin
-                  if j < !j1 + n_beta then
+                  if (!j1 < j) && (j <= !j1 + n_beta) then 
                     begin
                       incr j_a;
                       j0 := !j1;
@@ -433,9 +442,9 @@ let create_matrix_spin_computed f det_space =
                     end;
                   j1 := !j0 + n_beta;
                   h123 := h1 !j_a i_b;
-                  let j_b = j - !j0 - 1 in
-                  !h123 j_b
-                end
+                end;
+                let j_b = j - !j0 - 1 in
+                !h123 j_b
         end;
       !result
     in
diff --git a/CI/CIMatrixElementF12.ml b/CI/CIMatrixElementF12.ml
new file mode 100644
index 0000000..3b7c7c6
--- /dev/null
+++ b/CI/CIMatrixElementF12.ml
@@ -0,0 +1,147 @@
+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 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 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 f76f8d2..c41442a 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -24,11 +24,16 @@ let f12_integrals mo_basis =
           0.
        else
         begin
+          let ijkl = F12.get_phys two_e_ints i j k l
+          in
+          let ijlk = F12.get_phys two_e_ints i j l k
+          in
           if s' = Spin.other s then
-            0.5 *. (F12.get_phys two_e_ints i j k l)
+            (* Minus sign because we swap spin variables
+               instead of orbital variables *)
+            0.375 *. ijkl +. 0.125 *. ijlk 
           else
-            0.25 *. ((F12.get_phys two_e_ints i j k l) -.
-                      (F12.get_phys two_e_ints i j l k))
+            0.25 *. (ijkl -. ijlk) 
         end
     ) )
 
@@ -44,14 +49,6 @@ let h_ij mo_basis ki kj =
   |> List.hd
 
 
-let hf_ij mo_basis ki kj =
-  let integrals =
-    List.map (fun f -> f mo_basis)
-      [ CI.h_integrals ; f12_integrals ]
-  in
-  CIMatrixElement.make integrals ki kj 
-
-
 let f_ij mo_basis ki kj =
   let integrals =
     List.map (fun f -> f mo_basis)
@@ -60,8 +57,11 @@ let f_ij mo_basis ki kj =
   CIMatrixElement.make integrals ki kj 
   |> List.hd
 
+let hf_ij mo_basis ki kj =
+  [ h_ij mo_basis ki kj ; f_ij mo_basis ki kj ]
 
-let is_internal det_space =
+
+let is_a_double det_space =
     let mo_class      = DeterminantSpace.mo_class det_space in
     let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
     let m l =
@@ -83,8 +83,8 @@ let is_internal det_space =
       and b = Bitstring.logand aux_mask beta
       in
       match Bitstring.popcount a + Bitstring.popcount b with
-      | 1 | 2 -> false
-      | _ -> true
+      | 2 -> true
+      | _ -> false
 
 
 let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
@@ -107,13 +107,12 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
         |> DeterminantSpace.determinant_stream
 
       in
-      (* Select only singly and doubly excited determinants
-          wrt FCI space *)
+      (* Select only doubly excited determinants wrt FCI space *)
       Stream.from (fun _ ->
         try
           let rec result () = 
             let ki = Stream.next s in
-            if is_internal ci.CI.det_space ki then
+            if not (is_a_double ci.CI.det_space ki) then
                 result ()
             else
                 Some ki
@@ -168,17 +167,6 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
         if result = [] then None else Some result
         )
     in
-    let result = 
-      let x = 
-        try [ Stream.next out_dets_stream ]
-        with Stream.Failure -> failwith "Auxiliary basis set does not produce any excited determinant"
-      in
-      let m_H_aux, m_F_aux = make_h_and_f x in
-      let m_HF =
-        gemm m_H_aux m_F_aux ~transb:`T 
-      in
-      gemm m_HF f12_amplitudes
-    in
 
     let iteration input =
       Printf.printf ".%!";
@@ -188,6 +176,14 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
       in
       gemm m_HF f12_amplitudes
     in
+
+    let result = 
+      let x = 
+        try [ Stream.next out_dets_stream ]
+        with Stream.Failure -> failwith "Auxiliary basis set does not produce any excited determinant"
+      in
+      iteration x
+    in
  
     input_stream
     |> Farm.run ~ordered:false ~f:iteration
@@ -204,7 +200,7 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
 
 
 
-let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
+let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () =
 
   let f12 = Util.of_some @@ Simulation.f12 simulation in
   let mo_num = MOBasis.size mo_basis in
@@ -230,20 +226,17 @@ Printf.printf "Add aux basis\n%!";
     MOBasis.of_mo_basis s mo_basis
   in
   let () = 
-Printf.printf "F12 ints\n%!";
     ignore @@ MOBasis.f12_ints aux_basis
   in
   let () =
-Printf.printf "2e ints\n%!";
     ignore @@ MOBasis.two_e_ints aux_basis
   in
 
-Printf.printf "det space\n%!";
   let det_space = 
     DeterminantSpace.fci_f12_of_mo_basis  aux_basis  ~frozen_core  mo_num
   in
 
-  let ci = CI.make det_space in
+  let ci = CI.make ~n_states:state det_space in
 
   let ci_coef, ci_energy =
     let x = Lazy.force ci.eigensystem in
@@ -263,28 +256,31 @@ Printf.printf "det space\n%!";
       Lazy.force ci.CI.m_H
     in
 
-    let rec iteration ?(state=1) psi = 
+      
+    let rec iteration ~state psi = 
+      let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
       let delta = 
+         (* delta_i = {% $\sum_j c_j H_{ij}$ %} *)
          dressing_vector ~frozen_core aux_basis psi ci
+         |> Matrix.to_mat
       in
 
-      let f = 1.0 /. psi.{1,1} in
+      Printf.printf "Cmax : %e\n" psi.{column_idx,state};
+      Printf.printf "Norm : %e\n" (sqrt (gemm ~transa:`T delta delta).{state,state});
+
+      let f = 1.0 /. psi.{column_idx,state} in
       let delta_00 =
-        Util.list_range 2 (Mat.dim1 psi)
-        |> List.fold_left (fun accu i -> accu +.
-            (Matrix.get delta i 1) *. psi.{i,1} *. f) 0. 
+        (* Delta_00 = {% $\sum_{j \ne x} delta_j c_j / c_x$ %} *)
+        f *. ( (gemm ~transa:`T delta psi).{state,state} -.
+          delta.{column_idx,state} *. psi.{column_idx,state} )
       in
-      let delta = Matrix.to_mat delta in
-      delta.{1,1} <- delta.{1,1} -. delta_00;
+      Printf.printf "Delta_00 : %e %e\n" delta.{column_idx,state} delta_00;
+      delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00;
 
-(*------
-TODO SINGLE STATE HERE
-*)
-      let n_states = ci.CI.n_states in
       let diagonal =
         Vec.init (Matrix.dim1 m_H) (fun i ->
-          if i = 1 then
-            Matrix.get m_H i i +. delta.{1,1} *. f 
+          if i = column_idx then
+            Matrix.get m_H i i +. delta.{column_idx,state} *. f 
           else
             Matrix.get m_H i i
           )
@@ -295,24 +291,26 @@ TODO SINGLE STATE HERE
           Matrix.mm ~transa:`T m_H c
           |> Matrix.to_mat
         in
-        let c11 = Matrix.get c 1 1 in
+        let c11 = Matrix.get c column_idx state in
         Util.list_range 1 (Mat.dim1 w)
         |> List.iter (fun i ->
           let dci = 
-            delta.{i,1} *. f ;
+            delta.{i,state} *. f ;
           in
-          w.{i,1} <- w.{i,1} +. dci *.  c11;
-          if (i <> 1) then
-            w.{1,1} <- w.{1,1} +. dci *. (Matrix.get c i 1);
+          w.{i,state} <- w.{i,state} +. dci *.  c11;
+          if (i <> column_idx) then
+            w.{column_idx,state} <- w.{column_idx,state} +. dci *. (Matrix.get c i state);
           );
         Matrix.dense_of_mat w
       in
 
       let eigenvectors, eigenvalues = 
         Parallel.broadcast (lazy (
-          Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
+          Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states:state diagonal matrix_prod
           ))
       in
+      Vec.iter (fun energy -> Printf.printf "%g\t" energy) eigenvalues;
+      print_newline ();
 
       let conv =
         1.0 -. abs_float ( dot
@@ -320,18 +318,18 @@ TODO SINGLE STATE HERE
           (Mat.to_col_vecs eigenvectors).(0) )
       in
       if Parallel.master then
-        Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift 
+        Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{state} +. e_shift 
           +. Simulation.nuclear_repulsion simulation);
 
       if conv > threshold then
-        iteration eigenvectors
+        iteration ~state eigenvectors
       else
         let eigenvalues =
           Vec.map (fun x -> x +. e_shift) eigenvalues
         in
         eigenvectors, eigenvalues
     in
-    iteration ci_coef
+    iteration ~state ci_coef
 
   )
   in
diff --git a/Utils/Cholesky.ml b/Utils/Cholesky.ml
index de3d8a6..586f304 100644
--- a/Utils/Cholesky.ml
+++ b/Utils/Cholesky.ml
@@ -52,6 +52,10 @@ let full_ldl m_A =
 
 
 let pivoted_ldl threshold m_A =
+(** {% $P A P^\dagger = L D L^\dagger$. %}
+ Input : Matrix $A$
+ Output : Matrices $L, D, P$.
+*)
 
   let n = Mat.dim1 m_A in
   assert (Mat.dim2 m_A = n);
@@ -63,14 +67,14 @@ let pivoted_ldl threshold m_A =
 
   let v_D_inv  = Vec.make0 n in
 
-  let compute_d j =
-    let l_jk = 
-      Mat.col m_Lt j
+  let compute_d i =
+    let l_ik = 
+      Mat.col m_Lt i
     in
-    let l_jk__d_k = 
-      Vec.mul ~n:(j-1) l_jk v_D
+    let l_ik__d_k = 
+      Vec.mul ~n:(i-1) l_ik v_D
     in
-    m_A.{pi.(j),pi.(j)} -. dot ~n:(j-1) l_jk  l_jk__d_k
+    m_A.{pi.(i),pi.(i)} -. dot ~n:(i-1) l_ik  l_ik__d_k
   in
 
   let compute_l i =
@@ -86,8 +90,9 @@ let pivoted_ldl threshold m_A =
       v_D_inv.{j} *. ( m_A.{pi.(j),pi.(i)} -. dot ~n:(j-1) l_ik__d_k l_jk )
   in
 
-  let pivot i =
-    let rec aux (pos,value) i = 
+
+  let maxloc v i =
+    let rec aux pos value i =
       if i > n then
         pos
       else if v_D.{i} > value then
@@ -95,23 +100,36 @@ let pivoted_ldl threshold m_A =
       else
         aux pos value (i+1)
     in
-    let j = aux i v_D.{i} (i+1) in
-    let p_i, p_j = pi.(i), pi.(j) in
-    pi.(i) <- pj;
-    pi.(j) <- pi;
+    aux i v.{i} (i+1)
   in
 
-  for i=1 to n do
-    pivot i;
-    for j=1 to (i-1) do
-      m_Lt.{j,i} <- compute_l i j;
-    done;
-    let d_i = compute_d i in
-    v_D.{i} <- d_i;
-    v_D_inv.{i} <- 1. /. d_i;
-  done;
-  m_Lt, v_D
 
+  let pivot i =
+    let j = maxloc v_D i in
+    let p_i, p_j = pi.(i), pi.(j) in
+    pi.(i) <- p_j;
+    pi.(j) <- p_i;
+  in                                                                                          
+
+
+
+
+  let () = 
+    try
+      for i=1 to n do
+        pivot i;
+        for j=1 to (i-1) do
+          m_Lt.{j,i} <- compute_l i j;
+        done;
+        let d_i = compute_d i in
+        if abs_float d_i < threshold then
+          raise Exit;
+        v_D.{i} <- d_i;
+        v_D_inv.{i} <- 1. /. d_i;
+      done
+    with Exit -> ()
+  in
+  m_Lt, v_D, pi
 
 
 
@@ -121,6 +139,11 @@ let make_ldl ?(threshold=Constants.epsilon) m_A =
 
 
 
+
+
+
+
+
 let test_case () =
 
   let matrix_diff m_A m_B = 
@@ -132,7 +155,10 @@ let test_case () =
     dot v v
   in
 
-  let test_full () =
+  let m_A = Mat.random 1000 1000 in
+  let m_A = Mat.add m_A (Mat.transpose_copy m_A) in
+
+  let test_full_small () =
     let m_A =
       Mat.of_array [| [|   4. ;  12. ; -16. |] ;
                       [|  12. ;  37. ; -43. |] ;
@@ -154,11 +180,61 @@ let test_case () =
     in
 
     let m_Lt, v_D = full_ldl m_A in
-
     Alcotest.(check (float 1.e-15)) "full L" 0. (matrix_diff m_Lt m_Lt_ref);
-    Alcotest.(check (float 1.e-15)) "full D" 0. (vector_diff v_D v_D_ref)
+    Alcotest.(check (float 1.e-15)) "full D" 0. (vector_diff v_D v_D_ref);
+    let m_D = Mat.of_diag v_D in
+    let m_B = gemm ~transa:`T m_Lt @@ gemm m_D m_Lt in
+    Alcotest.(check (float 1.e-15)) "full L" 0. (matrix_diff m_A m_B);
+    ()
+  in
+
+  let test_full () =
+    let m_Lt, v_D = full_ldl m_A in
+    let m_D = Mat.of_diag v_D in
+    let m_B = gemm ~transa:`T m_Lt @@ gemm m_D m_Lt in
+    Alcotest.(check (float 1.e-15)) "full D" 0. (matrix_diff m_A m_B);
+    ()
+  in
+
+  let test_pivoted () =
+    let m_Lt, v_D, pi = pivoted_ldl 0. m_A  in
+    let n = Mat.dim1 m_A in
+    let m_P = Mat.make0 n n in
+    for i=1 to n do
+      m_P.{i,pi.(i)} <- 1.
+    done;
+    let m_D = Mat.of_diag v_D in
+    let m_B =
+      gemm ~transa:`T m_P  @@
+      gemm ~transa:`T m_Lt @@
+      gemm m_D @@
+      gemm m_Lt m_P
+    in
+    Alcotest.(check (float 1.e-14)) "pivoted D" 0. (matrix_diff m_A m_B);
+    ()
+  in
+
+  let test_truncated () =
+    let m_Lt, v_D, pi = pivoted_ldl 0.001 m_A  in
+    let n = Mat.dim1 m_Lt in
+    let m_P = Mat.make0 n n in
+    for i=1 to n do
+      m_P.{i,pi.(i)} <- 1.
+    done;
+    let m_D = Mat.of_diag v_D in
+    let m_B =
+      gemm ~transa:`T m_P  @@
+      gemm ~transa:`T m_Lt @@
+      gemm m_D @@
+      gemm m_Lt m_P
+    in
+    Alcotest.(check (float 1.e-3)) "full D" 0. (matrix_diff m_A m_B);
+    ()
   in
   [
+    "Small", `Quick, test_full_small;
     "Full", `Quick, test_full;
+    "Pivoted", `Quick, test_pivoted ;
+    "Truncated", `Quick, test_truncated ;
   ]
 
diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml
index bd52d7f..0af4e15 100644
--- a/Utils/Davidson.ml
+++ b/Utils/Davidson.ml
@@ -14,15 +14,10 @@ let make
   (* Size of the matrix to diagonalize *)
   let n = Vec.dim diagonal in
 
+  let m = n_states in
 
-  let m = (* Number of requested states *)
-    match guess with
-    | Some vectors -> (Mat.dim2 vectors) * n_states
-    | None -> n_states
-  in
-
-  (* Create guess vectors u, with randomly initialized unknown vectors. *)
-  let init_vectors =
+  (* Create guess vectors u, with unknown vectors initialized to unity. *)
+  let init_vectors m =
     let init_vector k =
       Vector.sparse_of_assoc_list n [ (k,1.0) ]
     in
@@ -31,6 +26,20 @@ let make
     |> Matrix.to_mat
   in
 
+  let guess =
+    match guess with
+    | Some vectors -> vectors 
+    | None -> init_vectors m
+  in
+
+  let guess =
+    if Mat.dim2 guess = n_states then guess
+    else
+      (Mat.to_col_vecs_list guess) @
+        (Mat.to_col_vecs_list (init_vectors (m-(Mat.dim2 guess))) )
+      |> Mat.of_col_vecs_list
+  in
+
   let pick_new u =
     Mat.to_col_vecs_list u
     |> Util.list_pack m
@@ -38,19 +47,14 @@ let make
     |> List.hd
   in
 
-  let u_new =
-   match guess with
-   | Some vectors -> Mat.to_col_vecs_list vectors 
-   | None -> Mat.to_col_vecs_list init_vectors
-  in
+  let u_new = Mat.to_col_vecs_list guess in
 
-  let rec iteration u u_new w iter =
+  let rec iteration u u_new w iter macro =
     (* u is a list of orthonormal vectors, on which the operator has
        been applied : w = op.u
-       u_new is a list of vector which will increase the size of the
+       u_new is a list of vectors which will increase the size of the
        space.
     *)
-
     (* Orthonormalize input vectors u_new *)
     let u_new_ortho =
       List.concat [u ; u_new]
@@ -59,7 +63,7 @@ let make
       |> pick_new
     in
 
-    (* Apply the operator the m last vectors *)
+    (* Apply the operator to the m last vectors *)
     let w_new =
       matrix_prod (
         u_new_ortho
@@ -101,17 +105,29 @@ let make
     (* Compute the residual as proposed new vectors *)
     let u_proposed = 
       Mat.init_cols n m (fun i k ->
-        (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. 
-          (max (diagonal.{i} -. lambda.{k}) 0.01) )
+          let delta = diagonal.{i} -. lambda.{k} in
+          let delta = 
+            if abs_float delta > 0.01 then delta
+            else if delta < 0. then -.0.01
+                 else 0.01
+          in
+        (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. delta )
       |> Mat.to_col_vecs_list
     in
 
 
     let residual_norms = List.map nrm2 u_proposed in
-    let residual_norm  = List.fold_left (fun accu i -> max accu i) 0. residual_norms in
+    let residual_norm  =
+      List.fold_left (fun accu i -> accu +. i *. i) 0. residual_norms
+      |> sqrt
+    in
 
     if Parallel.master then
-      Printf.printf "%3d %16.10f  %16.8e%!\n" iter lambda.{1} residual_norm;
+      begin
+        Printf.printf "%3d " iter;
+        Vec.iteri (fun i x -> if (i<=m) then Printf.printf "%16.10f  " x) lambda;
+        Printf.printf "%16.8e%!\n" residual_norm
+      end;
 
     (* Make new vectors sparse *)
     let u_proposed = 
@@ -125,20 +141,20 @@ let make
     in
 
 
-    if residual_norm > threshold then
-      let u_next, w_next, iter =
+    if residual_norm > threshold && macro > 0 then
+      let u_next, w_next, iter, macro =
         if iter = n_iter then
           m_new_U |> pick_new, 
           m_new_W |> pick_new,
-          0
+          0, (macro-1)
         else
-          u_next, w_next, iter
+          u_next, w_next, (iter+1), macro
       in
-      iteration u_next u_proposed w_next (iter+1)
+      iteration u_next u_proposed w_next iter macro
     else
       (m_new_U |> pick_new |> Mat.of_col_vecs_list), lambda
 
   in
-  iteration [] u_new [] 1
+  iteration [] u_new [] 1 5
 
 
diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml
index c77494b..f3ede9a 100644
--- a/Utils/FourIdxStorage.ml
+++ b/Utils/FourIdxStorage.ml
@@ -10,6 +10,7 @@ type array2 = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Ar
 
 type storage_t = 
   | Dense  of array2
+  | Svd of Vec.t * array2 * array2
   | Sparse of (int, float) Hashtbl.t
 
 type t =
@@ -86,6 +87,7 @@ let unsafe_get_four_index ~r1 ~r2 t =
   else
     match t.four_index with
     | Dense  a -> unsafe_get a (dense_index i k t.size) (sym_index j l)
+    | Svd _ -> assert false
     | Sparse a -> let key = key_of_indices ~r1 ~r2 in
                   try Hashtbl.find a key
                   with Not_found -> 0.
@@ -174,6 +176,7 @@ let unsafe_set_four_index ~r1 ~r2 ~value t =
                 end
   | Sparse a -> let key = key_of_indices ~r1 ~r2 in
                 Hashtbl.replace a key value
+  | Svd _ -> assert false
     
 
 let set_four_index ~r1 ~r2 ~value t = 
@@ -301,6 +304,7 @@ let get_chem_all_ij d ~k ~l  =
           Bigarray.reshape_2 result d.size d.size
       | Sparse a ->
           Mat.init_cols d.size d.size (fun i j -> get_chem d i j k l)
+      | Svd _ -> assert false
     
 
     
@@ -419,16 +423,12 @@ t
 
      
 
-let four_index_transform coef source =
+let four_index_transform_dense_sparse ds coef source =
+
+  let mo_num = Mat.dim2 coef in
+  let destination = create ~size:mo_num ds in
 
   let ao_num = Mat.dim1 coef in
-  let mo_num = Mat.dim2 coef in
-  let destination =
-      match source.four_index with
-      | Dense _  -> create ~size:mo_num `Dense
-      | Sparse _ -> create ~size:mo_num `Sparse
-  in
-
   let mo_num_2  = mo_num * mo_num in
   let ao_num_2  = ao_num * ao_num in
   let ao_mo_num = ao_num * mo_num in
@@ -512,3 +512,181 @@ let four_index_transform coef source =
   if Parallel.master then Printf.eprintf "\n%!";
   broadcast destination
 
+
+let svd_of_dense t =
+  let four_index =
+    let f =
+      match t.four_index with
+      | Dense a -> a
+      | _ -> invalid_arg "svd_of_dense expects a Dense storage"
+    in
+    let size = t.size in
+    let b = Mat.create (size*size) (size*size) in
+    for k = 1 to size do
+      for l = 1 to size do
+        let ac = sym_index k l in
+        lacpy f ~n:1 ~ac ~b ~bc:(dense_index k l size)
+        |> ignore;
+        lacpy f ~n:1 ~ac ~b ~bc:(dense_index l k size)
+        |> ignore
+      done
+    done;
+    let d, u, v = gesdd b in
+    let rank =
+      let w = copy d in
+      scal (1. /. d.{1}) w;
+      let rec aux i = 
+        if i = Vec.dim w then i else
+        if w.{i} < 1.e-15 then i else
+          aux (i+1)
+      in aux 1
+    in
+    let d = copy d ~n:rank in
+    let u = lacpy ~n:rank u in
+    let v = lacpy ~m:rank v in
+    Svd (d,u,v)
+  in
+  { t with four_index }
+
+
+let dense_of_svd t = 
+  let four_index =
+    let d, u, v =
+      match t.four_index with
+      | Svd (d, u, v) -> d, u, v
+      | _ -> invalid_arg "dense_of_svd expects a Svd storage"
+    in
+    let f = 
+      gemm u @@ gemm (Mat.of_diag d) v 
+    in
+    let size = t.size in
+    let b = Mat.create (size*size) ((size*(size+1))/2) in
+    for k = 1 to size do
+      for l = 1 to k do
+        let bc = sym_index k l in
+        lacpy f ~ac:(dense_index l k size) ~n:1 ~bc ~b
+        |> ignore
+      done
+    done;
+    Dense b
+  in
+  { t with four_index }
+
+
+let four_index_transform_svd coef source =
+  let t0 = Unix.gettimeofday () in
+  let svd = svd_of_dense source in
+  Printf.printf "Computed SVD in %f seconds\n%!" (Unix.gettimeofday () -. t0);
+
+  let u0, d0, v0 = 
+    match svd.four_index with
+    | Svd (d, u, v) -> u, d, (Mat.transpose_copy v)
+    | _ -> assert false
+  in
+
+  let t0 = Unix.gettimeofday () in
+
+  let ao_num = Mat.dim1 coef in
+  let mo_num = Mat.dim2 coef in
+  let svd_num = Vec.dim d0 in
+
+  let destination = create ~size:mo_num `Dense in
+
+  let mo_num_2  = mo_num * mo_num in
+
+  if Parallel.master then Printf.eprintf "4-idx transformation \n%!";
+
+  let u_vecs = Mat.to_col_vecs u0
+  and v_vecs = Mat.to_col_vecs v0
+  in
+
+
+  Printf.printf "%d %d %d %d\n" (Mat.dim1 u0) (Mat.dim2 u0)  (Mat.dim1 v0) (Mat.dim2 v0);
+  let task a =
+    let o =
+      Bigarray.reshape_2 (Bigarray.genarray_of_array1 u_vecs.(a-1)) ao_num ao_num
+    in
+
+    let u =
+      let p = gemm ~transa:`T coef @@ gemm o coef in
+      Bigarray.reshape_1 (Bigarray.genarray_of_array2 p) mo_num_2
+    in
+
+    let o =
+      Bigarray.reshape_2 (Bigarray.genarray_of_array1 v_vecs.(a-1)) ao_num ao_num
+    in
+
+    let v =
+      let p = gemm ~transa:`T coef @@ gemm o coef in
+      Bigarray.reshape_1 (Bigarray.genarray_of_array2 p) mo_num_2
+    in
+    (u, v)
+  in
+
+
+  (*
+  let four_index = Svd (d,u'',v'') in
+  let t = make_from_four_index
+
+  { source with
+    size = mo_num;
+    two_index;
+    two_index_anti;
+    three_index;
+    three_index_anti;
+    four_index
+  }
+  *)
+
+  let n = ref 0 in
+  let u_list = ref []
+  and v_list = ref []
+  in
+  Util.list_range 1 svd_num
+  |> Stream.of_list 
+  |> Farm.run ~f:task ~ordered:true ~comm:Parallel.Node.comm
+  |> Stream.iter (fun (u,v) ->
+     if Parallel.master then (incr n ; Printf.eprintf "\r%d / %d%!" !n mo_num);
+     u_list := u :: !u_list;
+     v_list := v :: !v_list);
+
+  let u = Mat.of_col_vecs_list (List.rev !u_list)
+  and v = Mat.of_col_vecs_list (List.rev !v_list)
+  in
+  Printf.printf "Computed 4-idx in %f seconds\n%!" (Unix.gettimeofday () -. t0);
+
+  let result = 
+    let p = gemm u @@ gemm ~transb:`T (Mat.of_diag d0) v in
+    Bigarray.reshape (Bigarray.genarray_of_array2 p) [| mo_num ; mo_num ; mo_num ; mo_num |]
+  in
+  for a=1 to mo_num do
+    for b=a to mo_num do
+      for c=1 to mo_num do
+        for d=c to mo_num do
+          let x = result.{a,b,c,d} in
+          if abs_float x > epsilon then
+            set_chem destination a b c d x;
+        done
+      done
+    done
+  done;
+
+  Printf.printf "Computed 4-idx in %f seconds\n%!" (Unix.gettimeofday () -. t0);
+
+  if Parallel.master then Printf.eprintf "\n%!";
+  broadcast destination
+
+
+
+
+
+
+
+
+
+let four_index_transform coef source =
+  match source.four_index with
+  | Sparse   _ -> four_index_transform_dense_sparse `Sparse coef source
+  | Svd      _ -> assert false (* four_index_transform_svd coef source *)
+  | Dense    _ -> four_index_transform_dense_sparse `Dense  coef source
+
diff --git a/Utils/Vector.ml b/Utils/Vector.ml
index fe1fd55..8905680 100644
--- a/Utils/Vector.ml
+++ b/Utils/Vector.ml
@@ -171,7 +171,7 @@ let axpy ?(threshold=epsilon) ?(alpha=1.) x y =
             in aux new_accu r1 r2
           | _ -> assert false
           end
-        | ({index=i ; value=x}::r1), []  -> aux ({index=i ; value=x}::accu) r1 []
+        | ({index=i ; value=x}::r1), []  -> aux ({index=i ; value=alpha *. x}::accu) r1 []
         | [] , ({index=j ; value=y}::r2) -> aux ({index=j ; value=y}::accu) [] r2
         | [] , [] -> {n ; v=List.rev accu}
       in
@@ -297,6 +297,7 @@ let test_case () =
     Alcotest.(check  bool)  "dense   dense   axpy"  true  (axpy ~alpha:3.  v1    v2    =  v6);
     Alcotest.(check  bool)  "dense   sparse  axpy"  true  (sub ~threshold:1.e-12 (axpy ~alpha:3.  v1    v2_s) v6_s = zero_s);
     Alcotest.(check  bool)  "sparse  dense   axpy"  true  (sub ~threshold:1.e-12 (axpy ~alpha:3.  v1_s  v2) v6_s = zero_s);
+
     Alcotest.(check  bool)  "sparse  sparse  axpy"  true  (sub ~threshold:1.e-12 (axpy ~alpha:3.  v1_s  v2_s) v6_s = zero_s);
   in
   let test_dot () =
diff --git a/run_fci_f12.ml b/run_fci_f12.ml
index 5a374a0..d7e3390 100644
--- a/run_fci_f12.ml
+++ b/run_fci_f12.ml
@@ -1,3 +1,5 @@
+open Lacaml.D
+
 let () =
   let open Command_line in
   begin
@@ -27,6 +29,10 @@ let () =
         { short='e' ; long="expo" ; opt=Optional;
           arg=With_arg "<float>";
           doc="Exponent of the Gaussian geminal"; } ;
+
+        { short='s' ; long="state" ; opt=Optional;
+          arg=With_arg "<int>";
+          doc="Requested state. Default is 1."; } ;
       ]
   end;
 
@@ -49,6 +55,12 @@ let () =
     | None   -> 1.0
   in
 
+  let state =
+    match Command_line.get "state" with
+    | Some x -> int_of_string x
+    | None -> 1
+  in
+
   let multiplicity =
     match Command_line.get "multiplicity" with
     | Some x -> int_of_string x
@@ -76,11 +88,21 @@ let () =
   in
 
   let fcif12 =
-    F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ()
+    F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ~state ()
   in
 
   let ci = F12CI.ci fcif12 in
-  Format.fprintf ppf "FCI energy     : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation);
+  Format.fprintf ppf "FCI energy     : ";
+  Vec.iteri (fun i x -> if i <= state then
+    Format.fprintf ppf "%20.16f@;  " (x +. Simulation.nuclear_repulsion simulation) )
+    (CI.eigenvalues ci);
+  Format.fprintf ppf "@.";
 
   let _, e_cif12 = F12CI.eigensystem fcif12 in
-  Format.fprintf ppf "FCI-F12 energy : %20.16f@." (e_cif12.{1} +. Simulation.nuclear_repulsion simulation);
+  Format.fprintf ppf "FCI-F12 energy : ";
+  Vec.iteri (fun i x -> if i <= state then
+    Format.fprintf ppf "%20.16f@;  " (x +. Simulation.nuclear_repulsion simulation) )
+    e_cif12;
+  Format.fprintf ppf "@."
+
+
diff --git a/run_integrals.ml b/run_integrals.ml
index 5b69c30..98dd70e 100644
--- a/run_integrals.ml
+++ b/run_integrals.ml
@@ -27,8 +27,9 @@ let run ~out =
     | Some x -> x
   in
 
+  let f12 = F12.gaussian_geminal 1.0 in
   let s =
-    Simulation.of_filenames ~nuclei:nuclei_file basis_file 
+    Simulation.of_filenames ~nuclei:nuclei_file ~f12 basis_file 
   in
 
   print_endline @@ Nuclei.to_string @@ Simulation.nuclei s;
@@ -45,10 +46,8 @@ let run ~out =
   NucInt.to_file ~filename:(out_file^".nuc") eN_ints;
   KinInt.to_file ~filename:(out_file^".kin") kin_ints;
   ERI.to_file ~filename:(out_file^".eri") ee_ints;
-  (*
-  let f12_ints  = AOBasis.f12_ints  ao_basis in
+  let f12_ints  = AOBasis.f12_ints ao_basis in
   F12.to_file ~filename:(out_file^".f12") f12_ints;
-  *)
   ()