diff --git a/Basis/Orthonormalization.ml b/Basis/Orthonormalization.ml
index a923e95..cf03f33 100644
--- a/Basis/Orthonormalization.ml
+++ b/Basis/Orthonormalization.ml
@@ -1,10 +1,7 @@
 open Util
 open Lacaml.D
 
-type t = 
-| Lowdin    of Mat.t
-| Canonical of Mat.t
-| Svd       of Mat.t
+type t = Mat.t
 
 module Am  = AngularMomentum
 module Bs  = Basis
@@ -12,38 +9,34 @@ module Cs  = ContractedShell
 
 
   
-let make_canonical_spherical ~thresh ~overlap basis =
-  let ao_num = Bs.size basis in
-  let cart_sphe = Mat.make ao_num ao_num 0. 
-  and i = ref 0 
-  and n = ref 0 in
-  Array.iter (fun shell ->
-     let submatrix =
-        SphericalToCartesian.matrix (Cs.ang_mom shell)
-      in
-      ignore @@ lacpy ~b:cart_sphe ~br:(!i+1) ~bc:(!n+1) submatrix;
-      i := !i + Mat.dim1 submatrix;
-      n := !n + Mat.dim2 submatrix;
-     ) (Bs.contracted_shells basis);
-  let s = gemm ~transa:`T  ~m:!n cart_sphe overlap in
-  let overlap = gemm s ~n:!n cart_sphe in
-  let s = canonical_ortho ~thresh ~overlap (Mat.identity !n) in
-  gemm cart_sphe ~k:!n s
+let make_canonical ~thresh ~basis ~cartesian ~overlap =
 
-
-
-let make_canonical ~cartesian ~thresh ~basis ~overlap =
-  let result = 
-    if cartesian then
-      canonical_ortho ~thresh ~overlap (Mat.identity @@ Mat.dim1 overlap)
-    else
-      match basis with
-      | None -> invalid_arg
-                "Basis.t is required when cartesian=false in make_canonical"
-      | Some basis -> 
-                make_canonical_spherical ~thresh ~overlap basis
+  let make_canonical_spherical basis =
+    let ao_num = Bs.size basis in
+    let cart_sphe = Mat.make ao_num ao_num 0. 
+    and i = ref 0 
+    and n = ref 0 in
+    Array.iter (fun shell ->
+      let submatrix =
+          SphericalToCartesian.matrix (Cs.ang_mom shell)
+        in
+        ignore @@ lacpy ~b:cart_sphe ~br:(!i+1) ~bc:(!n+1) submatrix;
+        i := !i + Mat.dim1 submatrix;
+        n := !n + Mat.dim2 submatrix;
+      ) (Bs.contracted_shells basis);
+    let s = gemm ~transa:`T  ~m:!n cart_sphe overlap in
+    let overlap = gemm s ~n:!n cart_sphe in
+    let s = canonical_ortho ~thresh ~overlap (Mat.identity !n) in
+    gemm cart_sphe ~k:!n s
   in
-  Canonical result
+
+  if cartesian then
+    canonical_ortho ~thresh ~overlap (Mat.identity @@ Mat.dim1 overlap)
+  else
+    match basis with
+    | None -> invalid_arg
+              "Basis.t is required when cartesian=false in make_canonical"
+    | Some basis -> make_canonical_spherical basis
 
 
 
@@ -59,16 +52,12 @@ let make_lowdin ~thresh ~overlap =
   let u_vec' =
     Mat.init_cols (Mat.dim1 u_vec) (Mat.dim2 u_vec) (fun i j -> u_vec.{i,j} *. u_val.{j})
   in
-  let result = 
-    gemm u_vec' ~transb:`T u_vec
-  in
-
-  Lowdin result
-     
+  gemm u_vec' ~transb:`T u_vec
 
 
-let make ~cartesian ?(thresh=1.e-12) ?basis overlap =
+
+let make ?(thresh=1.e-12) ?basis ~cartesian overlap =
   (*
   make_lowdin ~thresh ~overlap
   *)
-  make_canonical ~cartesian ~basis ~thresh ~overlap
+  make_canonical ~thresh ~basis ~cartesian ~overlap
diff --git a/Basis/Orthonormalization.mli b/Basis/Orthonormalization.mli
new file mode 100644
index 0000000..f184853
--- /dev/null
+++ b/Basis/Orthonormalization.mli
@@ -0,0 +1,12 @@
+(** Orthonormalization of the basis. *)
+
+type t = Lacaml.D.Mat.t
+
+val make: ?thresh:float -> ?basis:Basis.t -> cartesian:bool -> Overlap.t -> t
+(** Returns a matrix or orthonormal vectors in the basis. The vectors are
+    linearly independent up to a threshold [thresh]. If [cartesian] is
+    [false], the [basis] argument needs to be given and the space spanned by
+    the vectors is the same as the space spanned by the basis in spherical
+    coordinates (5d,7f,...). 
+*)
+
diff --git a/HartreeFock/RHF.ml b/HartreeFock/RHF.ml
index dfa5d67..481a11f 100644
--- a/HartreeFock/RHF.ml
+++ b/HartreeFock/RHF.ml
@@ -17,10 +17,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
 
   (* Orthogonalization matrix *)
   let m_X = 
-    match Lazy.force simulation.overlap_ortho with
-    | Lowdin    x -> x
-    | Canonical x -> x
-    | Svd       x -> x
+   Lazy.force simulation.overlap_ortho 
   in
 
 
diff --git a/Utils/Util.ml b/Utils/Util.ml
index 00b29a2..a6c186b 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -169,13 +169,13 @@ let array_product a =
    Array.fold_left ( *. ) 0. a
 
 
-let diagonalize_symm h = 
-  let v = lacpy h in
-  let w = Vec.create (Mat.dim1 h) in
+let diagonalize_symm m_H = 
+  let m_V = lacpy m_H in
+  let m_W = Vec.create (Mat.dim1 m_H) in
   let result = 
-    syevd ~vectors:true ~w v
+    syevd ~vectors:true ~w:m_W m_V
   in
-  v, result
+  m_V, result
 
 let xt_o_x ~o ~x =
   gemm o x
@@ -184,16 +184,19 @@ let xt_o_x ~o ~x =
 
 let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c =
   let d, u, _ = gesvd (lacpy overlap) in
-  let d0 = Vec.sqrt d in
-  let n = Vec.fold (fun accu x -> if x > thresh then accu + 1 else accu) 0 d in
-  let d = Vec.map (fun x ->
-    if x >= thresh then 1. /. x
-    else 0. ) ~y:d d0
+  let d_sqrt = Vec.sqrt d in
+  let n = (* Number of non-negligible singular vectors *)
+    Vec.fold (fun accu x -> if x > thresh then accu + 1 else accu) 0 d
   in
-  if n < Vec.dim d0 then 
+  let d_inv_sq = (* D^{-1/2} *)
+    Vec.map (fun x ->
+      if x >= thresh then 1. /. x
+      else 0. ) ~y:d d_sqrt
+  in
+  if n < Vec.dim d_sqrt then 
     Printf.printf "Removed linear dependencies below %f\n" (1. /. d.{n})
   ;
-  Mat.scal_cols u d ;
+  Mat.scal_cols u d_inv_sq ;
   gemm c u
 
 
diff --git a/Utils/Util.mli b/Utils/Util.mli
index a732351..fdd1b92 100644
--- a/Utils/Util.mli
+++ b/Utils/Util.mli
@@ -67,11 +67,19 @@ val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec
 (** Diagonalize a symmetric matrix. Returns the eigenvectors and the eigenvalues. *)
 
 val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
-(** Computes X{^T}.O.X *)
+(** Computes {% $\mathbf{X^\dag\, O\, X}$ %} *)
 
 val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat
-(** Canonical orthogonalization. [overlap] is the overlap matrix, and the last argument
-    contains the vectors to orthogonalize. *)
+(** 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}
+%}
+    
+    *)
 
 val debug_matrix: string -> Lacaml.D.mat -> unit
 (** Prints a matrix in stdout for debug *)