Documentation

Basis/Orthonormalization.ml

@@ -1,10 +1,7 @@
 open Util
 open Lacaml.D
 
-type t = 
+type t = Mat.t
-| Lowdin    of Mat.t
-| Canonical of Mat.t
-| Svd       of Mat.t
 
 module Am  = AngularMomentum
 module Bs  = Basis
@@ -12,7 +9,9 @@ module Cs  = ContractedShell
 
 
   
-let make_canonical_spherical ~thresh ~overlap basis =
+let make_canonical ~thresh ~basis ~cartesian ~overlap =
+
+  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 
@@ -29,21 +28,15 @@ let make_canonical_spherical ~thresh ~overlap basis =
     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
 
-
-
-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 -> 
+    | Some basis -> make_canonical_spherical basis
-                make_canonical_spherical ~thresh ~overlap basis
-  in
-  Canonical result
 
 
 
@@ -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
 
 
 
-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

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,...). 
+*)
+

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
+   Lazy.force simulation.overlap_ortho 
-    | Lowdin    x -> x
-    | Canonical x -> x
-    | Svd       x -> x
   in
 
 

Utils/Util.ml

@@ -169,13 +169,13 @@ let array_product a =
    Array.fold_left ( *. ) 0. a
 
 
-let diagonalize_symm h = 
+let diagonalize_symm m_H = 
-  let v = lacpy h in
+  let m_V = lacpy m_H in
-  let w = Vec.create (Mat.dim1 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 d_sqrt = Vec.sqrt d in
-  let n = Vec.fold (fun accu x -> if x > thresh then accu + 1 else accu) 0 d in
+  let n = (* Number of non-negligible singular vectors *)
-  let d = Vec.map (fun x ->
+    Vec.fold (fun accu x -> if x > thresh then accu + 1 else accu) 0 d
-    if x >= thresh then 1. /. x
-    else 0. ) ~y:d d0
   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
 
 

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
+(** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %},
-    contains the vectors to orthogonalize. *)
+    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 *)