diff --git a/SCF/RHF.ml b/SCF/RHF.ml
index c6a778d..e121805 100644
--- a/SCF/RHF.ml
+++ b/SCF/RHF.ml
@@ -3,7 +3,7 @@ open Lacaml.D
 open Simulation
 
 let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1)
-  ?threshold_SCF:(threshold_SCF=1.e-6) simulation =
+  ?threshold_SCF:(threshold_SCF=1.e-5) simulation =
 
   (* Number of occupied MOs *)
   let nocc =
@@ -34,8 +34,8 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
   and m_V = Lazy.force simulation.eN_ints
   in
 
-  (* Level shift *)
-  let m_LS =
+  (* Level shift in MO basis *)
+  let m_LSmo =
     Array.init (Mat.dim2 m_X) (fun i ->
       if i > nocc then level_shift else 0.)
     |> Vec.of_array
@@ -43,7 +43,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
   in
 
   (* SCF iterations *)
-  let rec loop nSCF iterations m_C diis =
+  let rec loop nSCF iterations energy_prev m_C diis =
 
       (* Density matrix over nocc occupied MOs *)
       let m_P =
@@ -58,41 +58,50 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
         x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange
       in
 
+      (* Add level shift in AO basis *)
+      let m_F =
+        let m_SC =
+          gemm  m_S  m_C
+        in
+        gemm  m_SC (gemm  m_LSmo  m_SC  ~transb:`T)
+        |> Mat.add m_F
+      in
+
+
+      (* Fock matrix in orthogonal basis *)
+      let m_F_ortho =
+        xt_o_x m_F m_X
+      in
+
       let error_fock = 
         let fps =
           gemm  m_F  (gemm  m_P  m_S)
         and spf =
           gemm  m_S  (gemm  m_P  m_F)
         in
-        Mat.sub  fps  spf
+        xt_o_x  (Mat.sub  fps  spf) m_X
       in
 
       let diis = 
-        DIIS.append ~p:(Mat.as_vec m_F) ~e:(Mat.as_vec error_fock) diis
+        DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
       in
 
       let m_F_diis =
         let x = 
           Bigarray.genarray_of_array1 (DIIS.next diis)
         in
-        Bigarray.reshape_2 x (Mat.dim1 m_F) (Mat.dim2 m_F)
+        Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
       in
 
 
-      (* Fock matrix in MO basis *)
-      let m_Fmo =
-        xt_o_x  m_F_diis  m_C
-        |> Mat.add m_LS
-      in
-
-      (* MOs in old MO basis *)
+      (* MOs in orthogonal MO basis *)
       let m_C', eigenvalues =
-        diagonalize_symm  m_Fmo
+        diagonalize_symm  m_F_diis
       in
 
       (* MOs in AO basis *)
       let m_C =
-        gemm  m_C  m_C'
+        gemm  m_X  m_C'
       in
 
       (* Hartree-Fock energy *)
@@ -103,7 +112,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
 
       (* Convergence criterion *)
       let error =
-        xt_o_x error_fock m_C
+        error_fock 
         |> Mat.as_vec 
         |> amax
         |> abs_float
@@ -117,10 +126,16 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
         eigenvalues.{nocc+1} -. eigenvalues.{nocc};
       in
 
-      Printf.printf "%d  %16.10f  %11.4e  %10.4f\n%!" nSCF energy error gap;
+      let () = 
+        match energy_prev with
+        | Some energy_prev ->
+          Printf.eprintf "%3d  %16.10f  %16.10f %11.4e  %10.4f\n%!" nSCF energy (energy -. energy_prev) error gap
+        | None -> 
+          Printf.eprintf "%3d  %16.10f  %16s %11.4e  %10.4f\n%!" nSCF energy "" error gap
+      in
 
       if not converged then
-          loop (nSCF+1) ( (energy, error, gap) :: iterations) m_C diis
+          loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C diis
       else
           let iterations = 
             List.rev ( (energy, error, gap) :: iterations )
@@ -159,7 +174,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
   in
 
   let diis = DIIS.make () in
-  loop 1 [] m_C diis
+  loop 1 [] None m_C diis
 
 
 
diff --git a/Utils/DIIS.ml b/Utils/DIIS.ml
index fe5c8ba..147a956 100644
--- a/Utils/DIIS.ml
+++ b/Utils/DIIS.ml
@@ -44,7 +44,7 @@ let next diis =
       let a = Mat.make (m+1) (m+1) 1. in
       a.{m+1,m+1} <- 0.;
       ignore @@ lacpy ~b:a (gemm ~transa:`T ~m ~n:m e e);
-      if sycon (lacpy a) > 1.e-10 then a
+      if sycon (lacpy a) > 1.e-14 then a
       else  aux (m-1)
     in
     aux diis.m