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/F12CI.ml b/CI/F12CI.ml
index f76f8d2..9beb0f4 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
+          and 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)
+            0.375 *. ijkl +.  0.125 *. ijlk
+            (*
+            0.25 *. ijkl
+            *)
           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
     ) )
 
@@ -204,7 +209,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
@@ -243,7 +248,7 @@ Printf.printf "det space\n%!";
     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 +268,26 @@ 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 = 
          dressing_vector ~frozen_core aux_basis psi ci
       in
 
-      let f = 1.0 /. psi.{1,1} in
+      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. 
+            (Matrix.get delta i state) *. psi.{i,state} *. f) 0. 
       in
       let delta = Matrix.to_mat delta in
-      delta.{1,1} <- delta.{1,1} -. 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 +298,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 state 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 +325,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/Davidson.ml b/Utils/Davidson.ml
index bd52d7f..d72d482 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 = (* Number of requested states *)
-    match guess with
-    | Some vectors -> (Mat.dim2 vectors) * n_states
-    | None -> n_states
-  in
+  let m = n_states in
 
   (* Create guess vectors u, with randomly initialized unknown vectors. *)
-  let init_vectors =
+  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]
@@ -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/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;
-  *)
   ()