diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml
index 2c57789..8f4d47e 100644
--- a/MOBasis/MOBasis.ml
+++ b/MOBasis/MOBasis.ml
@@ -1,5 +1,6 @@
 open Lacaml.D
 open Util
+open Constants
 
 (** One-electron orthogonal basis set, corresponding to Molecular Orbitals. *)
 
@@ -33,8 +34,7 @@ type t =
 
 
 let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = 
-  gemm ~transa:`T mo_coef   @@
-  gemm ao_matrix mo_coef  
+  xt_o_x ~x:mo_coef ~o:ao_matrix
 
 
 let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix =
@@ -56,14 +56,10 @@ let four_index_transform ~mo_coef eri_ao =
   let range_mo = list_range ~start:1 mo_num in
   let range_ao = list_range ~start:1 ao_num in
 
-  let u = 
-      Mat.create mo_num_2 mo_num
-  and o = 
-      Mat.create ao_num ao_num_2 
-  and p = 
-      Mat.create ao_num_2 mo_num
-  and q = 
-      Mat.create ao_mo_num mo_num
+  let u = Mat.create mo_num_2 mo_num
+  and o = Mat.create ao_num ao_num_2 
+  and p = Mat.create ao_num_2 mo_num
+  and q = Mat.create ao_mo_num mo_num
   in
   Printf.eprintf "Transforming %d integrals : %!" mo_num;
   List.iter (fun delta ->
@@ -71,39 +67,40 @@ let four_index_transform ~mo_coef eri_ao =
       Mat.fill u 0.;
 
       List.iter (fun l ->
+          if abs_float mo_coef.{l,delta} > epsilon then
+          begin
+            let jk = ref 0 in
+            List.iter (fun k ->
+                List.iter (fun j ->
+                    incr jk;
+                    ERI.get_chem_all_i eri_ao ~j ~k ~l
+                    |> Array.iteri (fun i x -> o.{i+1,!jk} <- x)
+                ) range_ao
+            ) range_ao;
+            (* o_i_jk *)
 
-          let jk = ref 0 in
-          List.iter (fun k ->
-              List.iter (fun j ->
-                  incr jk;
-                  ERI.get_chem_all_i eri_ao ~j ~k ~l
-                  |> Array.iteri (fun i x -> o.{i+1,!jk} <- x)
-              ) range_ao
-          ) range_ao;
-          (* o_i_jk *)
+            let p =
+              gemm ~transa:`T ~c:p o mo_coef
+              (* p_jk_alpha  = \sum_i o_i_jk c_i_alpha *)
+            in
+            let p' = 
+              Bigarray.reshape_2 (Bigarray.genarray_of_array2 p) ao_num ao_mo_num
+              (* p_j_kalpha *)
+            in
 
-          let p =
-            gemm ~transa:`T ~c:p o mo_coef
-            (* p_jk_alpha  = \sum_i o_i_jk c_i_alpha *)
-          in
-          let p' = 
-            Bigarray.reshape_2 (Bigarray.genarray_of_array2 p) ao_num ao_mo_num
-            (* p_j_kalpha *)
-          in
-
-          let q =
-            gemm ~transa:`T ~c:q p' mo_coef
-            (* q_kalpha_beta  = \sum_j p_j_kalpha c_j_beta *)
-          in
-          let q' =
-            Bigarray.reshape_2 (Bigarray.genarray_of_array2 q) ao_num mo_num_2
-            (* q_k_alphabeta  = \sum_j p_j_kalpha c_j_beta *)
-          in
-
-          ignore @@
-            gemm ~transa:`T ~beta:1. ~c:u q' mo_coef 
-            (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *)
+            let q =
+              gemm ~transa:`T ~c:q p' mo_coef
+              (* q_kalpha_beta  = \sum_j p_j_kalpha c_j_beta *)
+            in
+            let q' =
+              Bigarray.reshape_2 (Bigarray.genarray_of_array2 q) ao_num mo_num_2
+              (* q_k_alphabeta  = \sum_j p_j_kalpha c_j_beta *)
+            in
 
+            ignore @@
+              gemm ~transa:`T ~beta:1. ~alpha:mo_coef.{l,delta} ~c:u q' mo_coef ;
+              (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *)
+          end
       ) range_ao;
       let u =
         Bigarray.reshape 
@@ -119,7 +116,7 @@ let four_index_transform ~mo_coef eri_ao =
                     ERI.set_chem eri_mo alpha beta gamma delta x
               ) (list_range ~start:1 beta) 
           ) range_mo
-      ) (list_range ~start:1 delta) 
+      ) (list_range ~start:1 delta)  
   ) range_mo;
   Printf.eprintf "\n%!";
 
@@ -172,7 +169,35 @@ let of_rhf ~frozen_core hf =
       |> Vec.of_array
   in
   let mo_coef = hf.HF.eigenvectors in
-  make ~ao_basis ~mo_type ~mo_class ~mo_occupation ~mo_coef ()
+  let result = 
+    make ~ao_basis ~mo_type ~mo_class ~mo_occupation ~mo_coef ()
+  in
+
+  let () =
+    let e = ref 0. in
+    let t = KinInt.matrix (Lazy.force result.kin_ints) in
+    let v = NucInt.matrix (Lazy.force result.eN_ints) in
+    let g = Lazy.force result.ee_ints in
+    for i = 1 to 5 do
+      e := !e +. 2. *. (t.{i,i} +. v.{i,i})
+    done;
+    Printf.printf "Energy Mono = %20.15f\n" !e;
+    let e2 = ref 0. in
+    for i = 1 to 5 do
+      for j = i+1 to 5 do
+        e2 := !e2 +. 2. *. (ERI.get_phys g i j i j -.
+                            ERI.get_phys g i j j i)
+      done;
+    done;
+    for i = 1 to 5 do
+      for j = 1 to 5 do
+        e2 := !e2 +. ERI.get_phys g i j i j 
+      done;
+    done;
+    Printf.printf "Energy bi = %20.15f\n" !e2;
+    Printf.printf "Energy = %20.15f\n" (simulation.Si.nuclear_repulsion +. !e +. !e2)
+  in
+  result
 
 
 let of_hartree_fock ~frozen_core = function
diff --git a/Utils/Util.ml b/Utils/Util.ml
index 97168fe..25c11ec 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -172,7 +172,7 @@ let list_some l =
 let stream_range ?(start=0) n = 
   Stream.from (fun i ->
     let result = i+start in
-    if result < n then
+    if result <= n then
       Some result
     else None
   )
diff --git a/Utils/Util.mli b/Utils/Util.mli
index 4314561..61c9fac 100644
--- a/Utils/Util.mli
+++ b/Utils/Util.mli
@@ -64,11 +64,11 @@ val list_some : 'a option list -> 'a list
     the [Some]. *)
 
 val list_range : ?start:int -> int -> int list
-(** Returns a list [start ; start+1 ; ... ; start+(n-1)]. Default is [start=0]. *)
+(** Returns a list [start ; start+1 ; ... ; start+n]. Default is [start=0]. *)
 
 (** {2 Useful streams} *)
 val stream_range : ?start:int -> int -> int Stream.t
-(** Returns a stream <start ; start+1 ; ... ; start+(n-1)>. Default is [start=0]. *)
+(** Returns a stream <start ; start+1 ; ... ; start+n>. Default is [start=0]. *)
 
 
 (** {2 Linear algebra } *)