diff --git a/HartreeFock/Fock.ml b/HartreeFock/Fock.ml
index bb015b6..62749b3 100644
--- a/HartreeFock/Fock.ml
+++ b/HartreeFock/Fock.ml
@@ -4,7 +4,13 @@ open Constants
 
 
 
-type t = Mat.t
+type t =
+  {
+    fock     : Mat.t ;
+    core     : Mat.t ;
+    coulomb  : Mat.t ;
+    exchange : Mat.t ;
+  }
 
 let make ~density simulation =
   let m_P  = density
@@ -15,24 +21,40 @@ let make ~density simulation =
   let nBas = Mat.dim1 m_T
   in
 
-  let m_F = Mat.add m_T m_V in
+  let m_Hc = Mat.add m_T m_V 
+  and m_J = Mat.make0 nBas nBas 
+  and m_K = Mat.make0 nBas nBas 
+  in
   for sigma = 1 to nBas do
     for nu = 1 to nBas do
       for lambda = 1 to nBas do
         let p = m_P.{lambda,sigma} in
         if abs_float p > epsilon then
           for mu = 1 to nu do
-            m_F.{mu,nu} <- m_F.{mu,nu} +.  p *. 
-                          (ERI.get_phys m_G mu lambda nu sigma
-                           -. 0.5 *. ERI.get_phys m_G mu lambda sigma nu)
+            m_J.{mu,nu} <- m_J.{mu,nu} +.  p *. 
+                          ERI.get_phys m_G mu lambda nu sigma
+          done
+      done
+    done
+  done;
+  for nu = 1 to nBas do
+    for sigma = 1 to nBas do
+      for lambda = 1 to nBas do
+        let p = 0.5 *. m_P.{lambda,sigma} in
+        if abs_float p > epsilon then
+          for mu = 1 to nu do
+            m_K.{mu,nu} <- m_K.{mu,nu} -.  p *. 
+                          ERI.get_phys m_G mu lambda sigma nu
           done
       done
     done
   done;
   for nu = 1 to nBas do
     for mu = 1 to nu do
-      m_F.{nu,mu} <- m_F.{mu,nu}
+      m_J.{nu,mu} <- m_J.{mu,nu};
+      m_K.{nu,mu} <- m_K.{mu,nu};
     done
   done;
-  m_F
+  { fock = Mat.add m_Hc @@ Mat.add m_J m_K ;
+    core = m_Hc ; coulomb = m_J ; exchange = m_K }
 
diff --git a/HartreeFock/HartreeFock.ml b/HartreeFock/HartreeFock.ml
index 4673eac..da1420a 100644
--- a/HartreeFock/HartreeFock.ml
+++ b/HartreeFock/HartreeFock.ml
@@ -7,3 +7,5 @@ let make ?guess simulation =
   else
     invalid_arg "UHF or ROHF not implemented"
 
+
+let to_string = HartreeFock_type.to_string
diff --git a/HartreeFock/HartreeFock_type.ml b/HartreeFock/HartreeFock_type.ml
index 0d60f72..7e87868 100644
--- a/HartreeFock/HartreeFock_type.ml
+++ b/HartreeFock/HartreeFock_type.ml
@@ -1,11 +1,44 @@
 type t =
   {
-    guess        : Guess.t;
-    eigenvectors : Lacaml.D.Mat.t ;
-    eigenvalues  : Lacaml.D.Vec.t ;
-    energy       : float ;
-    iterations   : float array;
+    guess             : Guess.t;
+    eigenvectors      : Lacaml.D.Mat.t ;
+    eigenvalues       : Lacaml.D.Vec.t ;
+    nocc              : int;
+    energy            : float ;
+    nuclear_repulsion : float ;
+    kin_energy        : float ;
+    eN_energy         : float ;
+    coulomb_energy    : float ;
+    exchange_energy   : float ;
+    iterations        : (float * float * float) array;
   }
 
 
+let to_string hf_calc =
+  let sprintf = Printf.sprintf in
+  "
+  =====================================================
+                     Hartree-Fock                     
+  =====================================================
+      #       HF energy     Convergence    HOMO-LUMO   
+   ---------------------------------------------------" ::
+  Array.to_list (Array.mapi (fun i (e, c, g) ->
+  sprintf "    %5d   %13.8f   %11.4e   %11.4f  " (i+1) e c g) hf_calc.iterations)
+
+  @ [
+"  =====================================================";
+     sprintf "%20s energy  %16.10f" "One-electron" (hf_calc.kin_energy +. hf_calc.eN_energy) ;
+     sprintf "%20s energy  %16.10f" "Kinetic " hf_calc.kin_energy ;
+     sprintf "%20s energy  %16.10f" "Potential" hf_calc.eN_energy ;
+"   ---------------------------------------------------";
+     sprintf "%20s energy  %16.10f" "Two-electron" (hf_calc.coulomb_energy +. hf_calc.exchange_energy) ;
+     sprintf "%20s energy  %16.10f" "Coulomb" hf_calc.coulomb_energy ;
+     sprintf "%20s energy  %16.10f" "Exchange" hf_calc.exchange_energy ;
+"   ---------------------------------------------------";
+     sprintf "%20s energy  %16.10f" "Electronic" (hf_calc.energy -. hf_calc.nuclear_repulsion);
+     sprintf "%17s repulsion  %16.10f" "Nuclear" hf_calc.coulomb_energy ;
+     sprintf "%20s energy  %16.10f" "Hartree-Fock" hf_calc.energy ;
+"  =====================================================" ]
+  |> String.concat "\n"
+
 
diff --git a/HartreeFock/HartreeFock_type.mli b/HartreeFock/HartreeFock_type.mli
index 3a7527b..3d4fa01 100644
--- a/HartreeFock/HartreeFock_type.mli
+++ b/HartreeFock/HartreeFock_type.mli
@@ -2,11 +2,20 @@
 
 type t = 
   {
-    guess        : Guess.t;         (** Initial guess *)
-    eigenvectors : Lacaml.D.Mat.t ; (** Final eigenvectors *)
-    eigenvalues  : Lacaml.D.Vec.t ; (** Final eigenvalues  *)
-    energy       : float ;          (** Final energy *)
-    iterations   : float array;     (** Energies of all iterations *)
+    guess             : Guess.t;         (** Initial guess *)
+    eigenvectors      : Lacaml.D.Mat.t ; (** Final eigenvectors *)
+    eigenvalues       : Lacaml.D.Vec.t ; (** Final eigenvalues  *)
+    nocc              : int   ;          (** Number of occupied MOs *)
+    energy            : float ;          (** Final energy *)
+    nuclear_repulsion : float ;          (** Nucleus-Nucleus potential energy *)
+    kin_energy        : float ;          (** Kinetic energy *)
+    eN_energy         : float ;          (** Electron-nucleus potential energy *)
+    coulomb_energy    : float ;          (** Electron-Electron potential energy *)
+    exchange_energy   : float ;          (** Exchange energy *)
+    iterations        : (float * float * float) array;
+                                         (** Energy, convergence and HOMO-LUMO gap of all iterations *)
   }
 
+val to_string : t -> string
+(** Results of a Hartree-Fock calculation pretty-printed in a string. *)
 
diff --git a/HartreeFock/RHF.ml b/HartreeFock/RHF.ml
index 481a11f..e072328 100644
--- a/HartreeFock/RHF.ml
+++ b/HartreeFock/RHF.ml
@@ -10,6 +10,10 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
     simulation.electrons.Electrons.n_alpha
   in
 
+  let nuclear_repulsion = 
+    simulation.nuclear_repulsion
+  in
+
   (* Initial guess *)
   let guess = 
     Guess.make ~guess simulation 
@@ -21,18 +25,14 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
   in
 
 
-  (* Core Hamiltonian *)
-  let m_Hc =
-    let m_T = Lazy.force simulation.kin_ints
-    and m_V = Lazy.force simulation.eN_ints
-    in Mat.add m_T m_V
-  in
-
   (* Overlap matrix *)
   let m_S =
     Lazy.force  simulation.overlap
   in
 
+  let m_T = Lazy.force simulation.kin_ints
+  and m_V = Lazy.force simulation.eN_ints
+  in
 
   (* SCF iterations *)
   let rec loop nSCF iterations m_C =
@@ -43,8 +43,11 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
       in
 
       (* Fock matrix in AO basis *)
-      let m_F =
-        Fock.make  ~density:m_P  simulation
+      let m_F, m_Hc, m_J, m_K =
+        let x = 
+          Fock.make  ~density:m_P  simulation
+        in
+        x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange
       in
 
       (* Fock matrix in MO basis *)
@@ -64,7 +67,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
 
       (* Hartree-Fock energy *)
       let energy =
-        simulation.nuclear_repulsion +. 0.5 *.
+        nuclear_repulsion +. 0.5 *.
           Mat.gemm_trace  m_P  (Mat.add  m_Hc  m_F)
       in
 
@@ -81,24 +84,34 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
       in
 
       let converged =
-        nSCF = max_scf  ||  (abs_float commutator) < threshold_SCF
+        nSCF = max_scf  || (abs_float commutator) < threshold_SCF
       in
 
-      Printf.printf "%d  %16.10f  %10.4e\n%!" nSCF energy commutator;
+      let gap =
+        eigenvalues.{nocc+1} -. eigenvalues.{nocc};
+      in
+
+      Printf.printf "%d  %16.10f  %11.4e  %10.4f\n%!" nSCF energy commutator gap;
 
       if not converged then
-          loop (nSCF+1) (energy :: iterations) m_C
+          loop (nSCF+1) ( (energy, commutator, gap) :: iterations) m_C
       else
           let iterations = 
-            List.rev iterations
+            List.rev ( (energy, commutator, gap) :: iterations )
             |> Array.of_list
           in
           { HartreeFock_type.
+            nocc;
             guess ;
             eigenvectors = m_C ;
             eigenvalues ;
-            energy = iterations.(Array.length iterations - 1) ;
+            energy ;
+            nuclear_repulsion;
             iterations ;
+            kin_energy = Mat.gemm_trace m_P m_T;
+            eN_energy = Mat.gemm_trace m_P m_V;
+            coulomb_energy  = 0.5 *. Mat.gemm_trace m_P m_J;
+            exchange_energy = 0.5 *. Mat.gemm_trace m_P m_K;
           }
   in
 
diff --git a/Utils/Util.mli b/Utils/Util.mli
index fdd1b92..b94d358 100644
--- a/Utils/Util.mli
+++ b/Utils/Util.mli
@@ -58,7 +58,10 @@ val array_product : float array -> float
 
 
 (** {2 Extension of the List module} *)
+
 val list_some : 'a option list -> 'a list
+(** Filters out all [None] elements of the list, and returns the elements without
+    the [Some]. *)
 
 
 (** {2 Linear algebra } *)
@@ -72,11 +75,12 @@ val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
 val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat
 (** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %},
     and the last argument contains the vectors {% $\mathbf{C}$ %} to orthogonalize.
+
 {%
-\begin{align}
-\mathbf{U\, D\, V^\dag} & = \mathbf{S} \\
-\mathbf{C_\bot} & = \mathbf{C\, U\, D^{-1/2}}
-\end{align}
+$$
+\mathbf{C_\bot} = \mathbf{C\, U\, D^{-1/2}}, \;
+\mathbf{U\, D\, V^\dag} = \mathbf{S}
+$$
 %}
     
     *)
diff --git a/run_hartree_fock.ml b/run_hartree_fock.ml
index 144da8d..fba70fe 100644
--- a/run_hartree_fock.ml
+++ b/run_hartree_fock.ml
@@ -39,11 +39,9 @@ let run ~out =
     Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file 
   in
 
-  let _ =
-    HartreeFock.make s 
-  in
-  Printf.printf "Done.\n%!";
-  ()
+  HartreeFock.make s 
+  |> HartreeFock.to_string
+  |> print_endline
 
 
 let () =