Spherical works

Basis/Orthonormalization.ml

@@ -6,16 +6,42 @@ type t =
 | Canonical of Mat.t
 | Svd       of Mat.t
 
+module Am  = AngularMomentum
+module Bs  = Basis
+module Cs  = ContractedShell
 
 
   
-let make_canonical ~cartesian ~thresh ~overlap =
+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 ~cartesian ~thresh ~basis ~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)
+      match basis with
+      | None -> invalid_arg
+                "Basis.t is required when cartesian=false in make_canonical"
+      | Some basis -> 
+                make_canonical_spherical ~thresh ~overlap basis
   in
   Canonical result
 
@@ -41,8 +67,8 @@ let make_lowdin ~thresh ~overlap =
      
 
 
-let make ~cartesian ?(thresh=1.e-12) overlap =
+let make ~cartesian ?(thresh=1.e-12) ?basis overlap =
   (*
   make_lowdin ~thresh ~overlap
   *)
-  make_canonical ~cartesian ~thresh ~overlap
+  make_canonical ~cartesian ~basis ~thresh ~overlap

Simulation.ml

@@ -35,7 +35,7 @@ let make ?cartesian:(cartesian=true)
   {
     charge ; 
     basis ; nuclei ; electrons ; overlap ;
-    overlap_ortho = lazy (Orthonormalization.make ~cartesian (Lazy.force overlap));
+    overlap_ortho = lazy (Orthonormalization.make ~cartesian ~basis (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);

Utils/SphericalToCartesian.ml

@@ -5,10 +5,17 @@ module Am = AngularMomentum
 
 let matrix = function
   | Am.S -> Mat.make 1 1 1.
-  | Am.P -> Mat.make 3 3 1.
+  | Am.P -> 
+    begin
+      let result = Mat.make 3 3 0. in
+      result.{ 3, 1} <- 1.0;
+      result.{ 1, 2} <- 1.0;
+      result.{ 2, 3} <- 1.0;
+      result
+    end
   | Am.D -> 
     begin
-      let result = Mat.make 6 5 0.;
+      let result = Mat.make 6 5 0. in
       result.{ 1, 1} <- -0.5;
       result.{ 4, 1} <- -0.5;
       result.{ 6, 1} <- 1.0;
@@ -21,7 +28,7 @@ let matrix = function
     end
   | Am.F ->
     begin
-      let result = Mat.make 10 7 0.;
+      let result = Mat.make 10 7 0. in
       result.{ 3, 1} <- -0.67082039324993690892;
       result.{ 8, 1} <- -0.67082039324993690892;
       result.{10, 1} <- 1.0;
@@ -42,7 +49,7 @@ let matrix = function
     end
   | Am.G ->
     begin
-      let result = Mat.make 15 9 0.;
+      let result = Mat.make 15 9 0. in
       result.{ 1, 1} <- 0.375;
       result.{ 4, 1} <- 0.21957751641341996535;
       result.{ 6, 1} <- -0.87831006565367986142;
@@ -75,7 +82,7 @@ let matrix = function
     end
   | Am.H ->
     begin
-      let result = Mat.make 21 11 0.;
+      let result = Mat.make 21 11 0. in
       result.{ 3, 1} <- 0.625;
       result.{ 8, 1} <- 0.36596252735569994226;
       result.{10, 1} <- -1.0910894511799619063;
@@ -126,7 +133,7 @@ let matrix = function
     end
   |Am.I ->
     begin
-      let result = Mat.make 28 13 0.;
+      let result = Mat.make 28 13 0. in
       result.{ 1, 1} <- -0.3125;
       result.{ 4, 1} <- -0.16319780245846672329;
       result.{ 6, 1} <- 0.97918681475080033975;
@@ -201,7 +208,7 @@ let matrix = function
     end
   |Am.J ->
     begin
-      let result = Mat.make 36 15 0.;
+      let result = Mat.make 36 15 0. in
       result.{ 3, 1} <- -0.60670333962134435221;
       result.{ 8, 1} <- -0.31684048566533184861;
       result.{10, 1} <- 1.4169537279434593918;
@@ -308,7 +315,7 @@ let matrix = function
     end
   |Am.K ->
     begin
-      let result = Mat.make 45 17 0.;
+      let result = Mat.make 45 17 0. in
       result.{ 1, 1} <- 0.2734375;
       result.{ 4, 1} <- 0.13566299095694674896;
       result.{ 6, 1} <- -1.0853039276555739917;
@@ -454,7 +461,7 @@ let matrix = function
     end
   |Am.L ->
     begin
-      let result = Mat.make 55 19 0.;
+      let result = Mat.make 55 19 0. in
       result.{ 3, 1} <- 0.59686501473785067702;
       result.{ 8, 1} <- 0.29612797475437320937;
       result.{10, 1} <- -1.7657660842403202261;

Utils/SphericalToCartesian.mli

@@ -1,6 +1,6 @@
 (** Conversion from spherical coordinate to cartesian corrdinates. *)
 
-val matrix : AngularMomentum.t -> float array
+val matrix : AngularMomentum.t -> Lacaml.D.Mat.t
 (** Returns a transformation matrix to rotate between the basis of atom-centered
     spherical coordinates to x,y,z coordinates.
 

Utils/Util.ml

@@ -197,6 +197,10 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c =
   gemm c u
 
 
+let debug_matrix name a =
+  Format.printf "@[<2>%s =@\n@\n@[%a@]@]@\n@\n" name pp_mat a
+
+
 
 (** {2 Printers} *)
 

Utils/Util.mli

@@ -73,6 +73,8 @@ val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> La
 (** Canonical orthogonalization. [overlap] is the overlap matrix, and the last argument
     contains the vectors to orthogonalize. *)
 
+val debug_matrix: string -> Lacaml.D.mat -> unit
+(** Prints a matrix in stdout for debug *)
 
 (** {2 Printers} *)
 val pp_float_array_size : Format.formatter -> float array -> unit