Canonical orthogonalization OK in cart. TODO: spherical.

Basis/ERI.ml

@@ -77,6 +77,7 @@ let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs =
         | Some cspc -> 
           let cls = class_of_contracted_shell_pair_couple cspc in
           (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
+          (* TODO \sum_k |coef_k * integral_k| *)
         | None -> (pair, -1.) 
     ) shell_pairs
   |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff)

Basis/Orthonormalization.ml

@@ -7,7 +7,21 @@ type t =
 | Svd       of Mat.t
 
 
-let make_lowdin ?(thresh=1.e-12) ~overlap =
+  
+
+let make_canonical ~cartesian ~thresh ~overlap =
+  let result = 
+    if cartesian then
+      canonical_ortho ~thresh ~overlap (Mat.identity @@ Mat.dim1 overlap)
+    else
+      (* TODO *)
+      canonical_ortho ~thresh ~overlap (Mat.identity @@ Mat.dim1 overlap)
+  in
+  Canonical result
+
+
+
+let make_lowdin ~thresh ~overlap =
   
   let u_vec, u_val = diagonalize_symm overlap in
 
@@ -26,4 +40,9 @@ let make_lowdin ?(thresh=1.e-12) ~overlap =
   Lowdin result
      
 
-let make = make_lowdin
+
+let make ~cartesian ?(thresh=1.e-12) overlap =
+  (*
+  make_lowdin ~thresh ~overlap
+  *)
+  make_canonical ~cartesian ~thresh ~overlap

Basis/Overlap.ml

@@ -127,7 +127,6 @@ let of_basis basis =
 
 
 
-
 (** Write all overlap integrals to a file *)
 let to_file ~filename overlap =
 
@@ -146,4 +145,3 @@ let to_file ~filename overlap =
 
 
 
-

HartreeFock/RHF.ml

@@ -23,6 +23,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
     | Svd       x -> x
   in
 
+
   (* Core Hamiltonian *)
   let m_Hc =
     let m_T = Lazy.force simulation.kin_ints

Simulation.ml

@@ -9,9 +9,15 @@ type t = {
   kin_ints : KinInt.t lazy_t;
   ee_ints : ERI.t lazy_t;
   nuclear_repulsion : float;
+  cartesian : bool;
 }
 
-let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis =
+let make ?cartesian:(cartesian=true)
+         ?multiplicity:(multiplicity=1)
+         ?charge:(charge=0)
+         ~nuclei
+         basis
+  =
 
   (* Tune Garbage Collector *)
   let gc = Gc.get () in
@@ -29,20 +35,21 @@ let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis =
   {
     charge ; 
     basis ; nuclei ; electrons ; overlap ;
-    overlap_ortho = lazy (Orthonormalization.make (Lazy.force overlap));
+    overlap_ortho = lazy (Orthonormalization.make ~cartesian (Lazy.force overlap));
     eN_ints  = lazy (NucInt.of_basis_nuclei basis nuclei);
     kin_ints = lazy (KinInt.of_basis basis);
     ee_ints  = lazy (ERI.of_basis basis);
     nuclear_repulsion =  Nuclei.repulsion nuclei;
+    cartesian ;
   }
 
 
-let of_filenames ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis_filename =
+let of_filenames ?cartesian:(cartesian=false) ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis_filename =
   let nuclei =
     Nuclei.of_filename nuclei
   in
   let basis =
     Basis.of_nuclei_and_basis_filename ~nuclei basis_filename
   in
-  make ~charge ~multiplicity ~nuclei basis 
+  make ~cartesian ~charge ~multiplicity ~nuclei basis 
 

Utils/SphericalToCartesian.ml

@@ -647,7 +647,7 @@ let matrix = function
   | _ -> invalid_arg "Not implemented"
 
 
-let index nx ny nz =
+let index ~nx ~ny ~nz =
   let l = nx + ny + nz in
   ((l-nx)*(l-nx+1))/2 + nz + 1
 

Utils/SphericalToCartesian.mli

@@ -1,3 +1,5 @@
+(** Conversion from spherical coordinate to cartesian corrdinates. *)
+
 val matrix : AngularMomentum.t -> float array
 (** Returns a transformation matrix to rotate between the basis of atom-centered
     spherical coordinates to x,y,z coordinates.
@@ -12,12 +14,12 @@ val matrix : AngularMomentum.t -> float array
 ]}
 *)
 
-val index : int -> int -> int -> int
+val index : nx:int -> ny:int -> nz:int -> int
 (** Unique index given to a triplet of powers. Used to identify a cartesian shell.
 
     Example:
 {[
   let nx, ny, nz = 3, 2, 1 in
-  SphericalToCartesian.index nx ny nz -> 8
+  SphericalToCartesian.index ~nx ~ny ~nz -> 8
 ]}
 *)

Utils/Util.ml

@@ -182,6 +182,20 @@ let xt_o_x ~o ~x =
   |> gemm ~transa:`T 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
+  in
+  if n < Vec.dim d0 then 
+    Printf.printf "Removed linear dependencies below %f\n" (1. /. d.{n})
+  ;
+  Mat.scal_cols u d ;
+  gemm c u
+
 
 
 (** {2 Printers} *)

Utils/Util.mli

@@ -69,6 +69,10 @@ val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec
 val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
 (** Computes X{^T}.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. *)
+
 
 (** {2 Printers} *)
 val pp_float_array_size : Format.formatter -> float array -> unit