From c1a53f40e89d92d5d93e32f5c1db17c3cf375eec Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 13 Jun 2018 19:03:42 +0200
Subject: [PATCH] Added MatrixOnBasis

---
 Basis/KinInt.ml              |  1 +
 Basis/KinInt.mli             |  9 ++-------
 Basis/MatrixOnBasis.mli      | 14 ++++++++++++++
 Basis/NucInt.ml              |  5 +++++
 Basis/NucInt.mli             |  7 +++----
 Basis/Orthonormalization.ml  | 15 ++++++++++-----
 Basis/Orthonormalization.mli |  4 +++-
 Basis/Overlap.ml             |  3 +++
 Basis/Overlap.mli            |  6 ++++++
 SCF/Fock.ml                  |  6 +++---
 SCF/Guess.ml                 | 11 ++++++-----
 SCF/RHF.ml                   |  6 ++++--
 Utils/NonNegativeFloat.mli   |  1 +
 13 files changed, 61 insertions(+), 27 deletions(-)
 create mode 100644 Basis/MatrixOnBasis.mli
 create mode 100644 Basis/Overlap.mli

diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
index c5186ea..daff5f8 100644
--- a/Basis/KinInt.ml
+++ b/Basis/KinInt.ml
@@ -12,6 +12,7 @@ module Ps  = PrimitiveShell
 module Psp = PrimitiveShellPair
 
 type t = Mat.t
+external matrix : t -> Mat.t = "%identity"
 
 let cutoff = integrals_cutoff
 
diff --git a/Basis/KinInt.mli b/Basis/KinInt.mli
index dbccc0b..30432e6 100644
--- a/Basis/KinInt.mli
+++ b/Basis/KinInt.mli
@@ -1,4 +1,4 @@
-(** Electronic kinetic energy integrals, expressed as a matrix in a {Basis.t}.
+(** Electronic kinetic energy integrals, expressed as a matrix in a {!Basis.t}.
 
 {%
 $$
@@ -7,11 +7,6 @@ $$
 %}
 *)
 
-type t = Lacaml.D.Mat.t
+include module type of MatrixOnBasis
 
-val of_basis : Basis.t -> t
-(** Build from a {Basis.t}. *)
-
-val to_file : filename:string -> t -> unit
-(** Write the integrals in a file. *)
 
diff --git a/Basis/MatrixOnBasis.mli b/Basis/MatrixOnBasis.mli
new file mode 100644
index 0000000..c8955e8
--- /dev/null
+++ b/Basis/MatrixOnBasis.mli
@@ -0,0 +1,14 @@
+(** Signature for matrices built on the {!Basis.t} *)
+
+open Lacaml.D
+
+type t 
+
+val of_basis : Basis.t -> t
+(** Build from a {!Basis.t}. *)
+
+val to_file : filename:string -> t -> unit
+(** Write the integrals in a file. *)
+
+val matrix : t -> Mat.t
+(** Returns the matrix suitable for Lacaml *)
diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index dc92009..c5baff4 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -5,6 +5,7 @@ open Constants
 open Lacaml.D
 
 type t = Mat.t
+external matrix : t -> Mat.t = "%identity"
 
 module Am = AngularMomentum
 module Bs = Basis
@@ -108,3 +109,7 @@ let to_file ~filename eni_array =
   done;
   close_out oc
 
+
+let of_basis basis =
+  invalid_arg "of_basis_nuclei should be called for NucInt"
+
diff --git a/Basis/NucInt.mli b/Basis/NucInt.mli
index 28cab01..8ea1860 100644
--- a/Basis/NucInt.mli
+++ b/Basis/NucInt.mli
@@ -1,4 +1,4 @@
-(** Electron-Nucleus attractive potential integrals, expressed as a matrix in a {Basis.t}.
+(** Electron-Nucleus attractive potential integrals, expressed as a matrix in a {!Basis.t}.
 
 {% $$
 V_{ij} = \left \langle \chi_i \left| \sum_A \frac{-Z_A}{ | \mathbf{r} - \mathbf{R}_A |} \right| \chi_j \right \rangle
@@ -6,10 +6,9 @@ $$ %}
 
 *)
 
-type t = Lacaml.D.Mat.t
+include module type of MatrixOnBasis
 
 val of_basis_nuclei : basis:Basis.t -> Nuclei.t -> t
 (** Build from a {Basis.t} and the nuclei (coordinates and charges). *)
 
-val to_file : filename:string -> t -> unit
-(** Write the integrals in a file. *)
+
diff --git a/Basis/Orthonormalization.ml b/Basis/Orthonormalization.ml
index cf03f33..37b4f3f 100644
--- a/Basis/Orthonormalization.ml
+++ b/Basis/Orthonormalization.ml
@@ -6,11 +6,14 @@ type t = Mat.t
 module Am  = AngularMomentum
 module Bs  = Basis
 module Cs  = ContractedShell
+module Ov  = Overlap
 
 
   
 let make_canonical ~thresh ~basis ~cartesian ~overlap =
 
+  let overlap_matrix =  Ov.matrix overlap in
+
   let make_canonical_spherical basis =
     let ao_num = Bs.size basis in
     let cart_sphe = Mat.make ao_num ao_num 0. 
@@ -24,14 +27,14 @@ let make_canonical ~thresh ~basis ~cartesian ~overlap =
         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
+    let s = gemm ~transa:`T  ~m:!n cart_sphe overlap_matrix in
+    let overlap_matrix = gemm s ~n:!n cart_sphe in
+    let s = canonical_ortho ~thresh ~overlap:overlap_matrix (Mat.identity !n) in
     gemm cart_sphe ~k:!n s
   in
 
   if cartesian then
-    canonical_ortho ~thresh ~overlap (Mat.identity @@ Mat.dim1 overlap)
+    canonical_ortho ~thresh ~overlap:overlap_matrix (Mat.identity @@ Mat.dim1 overlap_matrix)
   else
     match basis with
     | None -> invalid_arg
@@ -42,7 +45,8 @@ let make_canonical ~thresh ~basis ~cartesian ~overlap =
 
 let make_lowdin ~thresh ~overlap =
   
-  let u_vec, u_val = diagonalize_symm overlap in
+  let overlap_matrix = Ov.matrix overlap in
+  let u_vec, u_val = diagonalize_symm overlap_matrix in
 
   Vec.iter (fun x -> if x < thresh then
     invalid_arg (__FILE__^": make_lowdin") ) u_val;
@@ -61,3 +65,4 @@ let make ?(thresh=1.e-12) ?basis ~cartesian overlap =
   make_lowdin ~thresh ~overlap
   *)
   make_canonical ~thresh ~basis ~cartesian ~overlap
+
diff --git a/Basis/Orthonormalization.mli b/Basis/Orthonormalization.mli
index f184853..f6adc77 100644
--- a/Basis/Orthonormalization.mli
+++ b/Basis/Orthonormalization.mli
@@ -1,6 +1,8 @@
 (** Orthonormalization of the basis. *)
 
-type t = Lacaml.D.Mat.t
+open Lacaml.D
+
+type t = 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
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index 3ab474b..7451cbf 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -3,6 +3,8 @@ open Constants
 open Lacaml.D
 
 type t = Mat.t
+external matrix : t -> Mat.t = "%identity"
+
 
 module Am  = AngularMomentum
 module Bs  = Basis
@@ -12,6 +14,7 @@ module Csp = ContractedShellPair
 module Po  = Powers
 module Psp = PrimitiveShellPair
 
+
 let cutoff = integrals_cutoff
 
 let to_powers x = 
diff --git a/Basis/Overlap.mli b/Basis/Overlap.mli
new file mode 100644
index 0000000..57be445
--- /dev/null
+++ b/Basis/Overlap.mli
@@ -0,0 +1,6 @@
+(** Overlap matrix of the atomic orbitals.
+
+{% $$ \langle \chi_i | \chi_j \rangle $$ %}
+*)
+
+include module type of MatrixOnBasis
diff --git a/SCF/Fock.ml b/SCF/Fock.ml
index b4c4d9b..0f51f20 100644
--- a/SCF/Fock.ml
+++ b/SCF/Fock.ml
@@ -15,9 +15,9 @@ module Ao = AOBasis
 
 let make ~density ao_basis =
   let m_P  = density
-  and m_T  = Lazy.force ao_basis.Ao.kin_ints
-  and m_V  = Lazy.force ao_basis.Ao.eN_ints
-  and m_G  = Lazy.force ao_basis.Ao.ee_ints
+  and m_T  = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix 
+  and m_V  = Lazy.force ao_basis.Ao.eN_ints  |> NucInt.matrix
+  and m_G  = Lazy.force ao_basis.Ao.ee_ints  
   in
   let nBas = Mat.dim1 m_T
   in
diff --git a/SCF/Guess.ml b/SCF/Guess.ml
index 3b5c4a4..e347ccf 100644
--- a/SCF/Guess.ml
+++ b/SCF/Guess.ml
@@ -9,10 +9,11 @@ type t = guess
 
 module Ao = AOBasis
 module El = Electrons
+module Ov = Overlap
 
 let hcore_guess ao_basis = 
-  let eN_ints  = Lazy.force ao_basis.Ao.eN_ints
-  and kin_ints = Lazy.force ao_basis.Ao.kin_ints
+  let eN_ints  = Lazy.force ao_basis.Ao.eN_ints  |> NucInt.matrix
+  and kin_ints = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix
   in
   Mat.add eN_ints kin_ints
 
@@ -20,9 +21,9 @@ let hcore_guess ao_basis =
 let huckel_guess ao_basis =
   let c = 0.5 *. 1.75 in
   let ao_num = Basis.size ao_basis.Ao.basis in
-  let eN_ints  = Lazy.force ao_basis.Ao.eN_ints
-  and kin_ints = Lazy.force ao_basis.Ao.kin_ints
-  and overlap  = Lazy.force ao_basis.Ao.overlap
+  let eN_ints  = Lazy.force ao_basis.Ao.eN_ints  |> NucInt.matrix
+  and kin_ints = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix
+  and overlap  = Lazy.force ao_basis.Ao.overlap  |> Ov.matrix
   and m_X = Lazy.force ao_basis.Ao.ortho
   in
   let diag = Array.init (ao_num+1) (fun i -> if i=0 then 0. else
diff --git a/SCF/RHF.ml b/SCF/RHF.ml
index 5200afb..04894e9 100644
--- a/SCF/RHF.ml
+++ b/SCF/RHF.ml
@@ -4,6 +4,7 @@ open Lacaml.D
 module Si = Simulation
 module El = Electrons
 module Ao = AOBasis
+module Ov = Overlap
 
 let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1)
   ?threshold_SCF:(threshold_SCF=1.e-5) simulation =
@@ -35,10 +36,11 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
   (* Overlap matrix *)
   let m_S =
     Lazy.force  ao_basis.Ao.overlap
+    |> Ov.matrix
   in
 
-  let m_T = Lazy.force ao_basis.Ao.kin_ints
-  and m_V = Lazy.force ao_basis.Ao.eN_ints
+  let m_T = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix
+  and m_V = Lazy.force ao_basis.Ao.eN_ints  |> NucInt.matrix
   in
 
   (* Level shift in MO basis *)
diff --git a/Utils/NonNegativeFloat.mli b/Utils/NonNegativeFloat.mli
index e68fd79..1bd9c06 100644
--- a/Utils/NonNegativeFloat.mli
+++ b/Utils/NonNegativeFloat.mli
@@ -1,6 +1,7 @@
 (** Floats >= 0. *)
 type t = private float
 val of_float : float -> t
+val unsafe_of_float : float -> t
 val to_float : t -> float
 val to_string : t -> string
 val of_string : string -> t