Added MatrixOnBasis

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
 

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

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 *)

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"
+

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

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
+

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

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 = 

Basis/Overlap.mli

@@ -0,0 +1,6 @@
+(** Overlap matrix of the atomic orbitals.
+
+{% $$ \langle \chi_i | \chi_j \rangle $$ %}
+*)
+
+include module type of MatrixOnBasis

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

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

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 *)

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