Added Hartree_Fock

ao/lib/ao_dim.ml

@@ -0,0 +1 @@
+type t

ao/lib/ao_dim.mli

@@ -0,0 +1 @@
+type t

ao/lib/basis.ml

@@ -9,6 +9,8 @@ type t =
     cartesian : bool
   }
 
+type ao = Ao_dim.t
+
 let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false)
     ~nuclei filename =
   match kind with

ao/lib/basis.mli

@@ -10,6 +10,7 @@ type basis =
   | Gaussian of Basis_gaussian.t
 
 type t 
+type ao = Ao_dim.t
 
 (** {1 Accessors} *)
 
@@ -19,37 +20,37 @@ val size : t -> int
 val ao_basis : t -> basis
 (** One-electron basis set *)
 
-val overlap : t -> ('a,'a) Matrix.t
+val overlap : t -> (ao,ao) Matrix.t
 (** Overlap matrix *)
 
-val multipole : t -> ('a,'a) Matrix.t array
+val multipole : t -> (ao,ao) Matrix.t array
 (** Multipole matrices *)
 
-val ortho : t -> ('a,'a) Matrix.t
+val ortho : t -> (ao,'a) Matrix.t
 (** Orthonormalization matrix of the overlap *)
 
-val eN_ints : t -> ('a,'a) Matrix.t
+val eN_ints : t -> (ao,ao) Matrix.t
 (** Electron-nucleus potential integrals *)
 
-val kin_ints : t -> ('a,'a) Matrix.t
+val kin_ints : t -> (ao,ao) Matrix.t
 (** Kinetic energy integrals *)
 
-val ee_ints : t -> 'a Four_idx_storage.t
+val ee_ints : t -> ao Four_idx_storage.t
 (** Electron-electron potential integrals *)
 
-val ee_lr_ints : t -> 'a Four_idx_storage.t
+val ee_lr_ints : t -> ao Four_idx_storage.t
 (** Electron-electron long-range potential integrals *)
 
-val f12_ints : t -> 'a Four_idx_storage.t
+val f12_ints : t -> ao Four_idx_storage.t
 (** Electron-electron potential integrals *)
 
-val f12_over_r12_ints : t -> 'a Four_idx_storage.t
+val f12_over_r12_ints : t -> ao Four_idx_storage.t
 (** Electron-electron potential integrals *)
 
 val cartesian : t -> bool
 (** If true, use cartesian Gaussians (6d, 10f, ...) *)
 
-val values : t -> Coordinate.t -> 'a Vector.t
+val values : t -> Coordinate.t -> ao Vector.t
 (** Values of the AOs evaluated at a given point *)
 
 

ao/test/ao_basis.ml

@@ -33,21 +33,21 @@ let make_tests name t =
 
   let test_overlap () =
     let reference =
-      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref")
+      Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref")
     in
     check_matrix "Overlap" (overlap t) reference
   in
 
   let test_eN_ints () =
     let reference =
-      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref")
+      Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref")
     in
     check_matrix "eN_ints" (eN_ints t) reference
   in
 
   let test_kin_ints () =
     let reference =
-      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^  "_kin.ref")
+      Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^  "_kin.ref")
     in
     check_matrix "kin_ints" (kin_ints t) reference
   in

ao/test/ao_basis_gaussian.ml

@@ -35,21 +35,21 @@ let make_tests name t =
 
   let test_overlap () =
     let reference =
-      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref")
+      Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref")
     in
     check_matrix "Overlap" (overlap t) reference
   in
 
   let test_eN_ints () =
     let reference =
-      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref")
+      Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref")
     in
     check_matrix "eN_ints" (eN_ints t) reference
   in
 
   let test_kin_ints () =
     let reference =
-      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^  "_kin.ref")
+      Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^  "_kin.ref")
     in
     check_matrix "kin_ints" (kin_ints t) reference
   in
@@ -69,7 +69,6 @@ let make_tests name t =
     in
     check_eri (ee_lr_ints t) reference
   in
-
   [
     "Overlap",     `Quick,  test_overlap;
     "eN_ints",     `Quick,  test_eN_ints;

linear_algebra/lib/four_idx_storage.ml

@@ -426,13 +426,13 @@ let four_index_transform_dense_sparse ds coef source =
               Matrix.gemm_inplace ~transa:`T ~c:p o coef;
               
               (* p_j_kalpha *)
-              let p' = Matrix.reshape p ao_num ao_mo_num in
+              let p' = Matrix.reshape_inplace ao_num ao_mo_num p in
 
               (* q_kalpha_beta  = \sum_j p_j_kalpha c_j_beta *)
               Matrix.gemm_inplace ~transa:`T ~c:q p' coef;
               
               (* q_k_alphabeta  = \sum_j p_j_kalpha c_j_beta *)
-              let q' = Matrix.reshape q ao_num mo_num_2 in
+              let q' = Matrix.reshape_inplace ao_num mo_num_2 q in
 
               (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *)
               Matrix.gemm_inplace ~transa:`T ~beta:1. ~alpha:coefx.{l,delta} ~c:u q' coef ;

linear_algebra/lib/linear_algebra.ml

@@ -0,0 +1,10 @@
+let (%.) = Vector.(%.)
+let (%:) = Matrix.(%:)
+
+module Conventions = Conventions
+module Diis = Diis
+module Four_idx_storage = Four_idx_storage
+module Matrix = Matrix
+module Orthonormalization = Orthonormalization
+module Spherical_to_cartesian = Spherical_to_cartesian
+module Vector = Vector

linear_algebra/lib/matrix.ml

@@ -48,7 +48,7 @@ external of_bigarray_inplace : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigar
 let to_bigarray t = lacpy t
 let of_bigarray t = lacpy t
 
-let reshape a m n =
+let reshape_inplace m n a =
   assert ( (dim1 a) * (dim2 a) = m*n);
   let b = 
     to_bigarray a 
@@ -56,6 +56,14 @@ let reshape a m n =
   in
   Bigarray.reshape_2 b m n
 
+let reshape_vec_inplace m n v =
+  assert ( Vector.dim v = m*n);
+  let b = 
+    Vector.to_bigarray_inplace v 
+    |> Bigarray.genarray_of_array1
+  in
+  Bigarray.reshape_2 b m n
+
 let col_inplace t j =
   Mat.col t j
   |> Vector.of_bigarray_inplace
@@ -200,6 +208,19 @@ let of_col_vecs_list a =
   |> List.rev
   |> Mat.of_col_vecs_list 
 
+let of_array a =
+  Bigarray.Array2.of_array Bigarray.Float64 Bigarray.fortran_layout a
+    
+let to_array a =
+  let result = Array.make_matrix (Mat.dim1 a) (Mat.dim2 a) 0. in
+  for i=1 to Mat.dim1 a do
+    let result_i = result.(i-1) in
+    for j=1 to Mat.dim2 a do
+      result_i.(j-1) <- a.{i,j}
+    done
+  done;
+  result
+    
 let normalize_mat_inplace t = 
   let norm = Mat.as_vec t |> nrm2 in
   Mat.scal norm t
@@ -227,6 +248,9 @@ let x_o_xt ~o ~x =
   gemm o x ~transb:`T
   |> gemm x
 
+let amax t =
+  Mat.as_vec t |> amax 
+    
 let pp ppf m =
   let rows = Mat.dim1 m
   and cols = Mat.dim2 m
@@ -263,7 +287,7 @@ let outer_product ?(alpha=1.0) u v =
   outer_product_inplace m ~alpha u v;
   m
 
-let matrix_of_file filename =
+let of_file filename =
   let ic = Scanf.Scanning.open_in filename in
   let rec read_line accu =
     let result =
@@ -290,9 +314,9 @@ let matrix_of_file filename =
   result
 
 
-let sym_matrix_of_file filename =
+let sym_of_file filename =
   let result =
-    matrix_of_file filename
+    of_file filename
   in
   for j=1 to Mat.dim1 result do
     for i=1 to j do

linear_algebra/lib/matrix.mli

@@ -17,9 +17,12 @@ val make0 : int -> int -> ('a,'b) t
 val create : int -> int -> ('a,'b) t
 (** Creates an uninitialized matrix. *)
 
-val reshape : ('a,'b) t -> int -> int -> ('c,'d) t
+val reshape_inplace : int -> int -> ('a,'b) t -> ('c,'d) t
 (** Changes the dimensions of the matrix *)
 
+val reshape_vec_inplace : int -> int -> ('a*'b) Vector.t -> ('a,'b) t
+(** Reshapres a vector into a matrix *)
+
 val init_cols : int -> int -> (int -> int -> float) -> ('a,'b) t
 (** Creates an uninitialized matrix. *)
 
@@ -59,6 +62,9 @@ val mul : ('a,'b) t -> ('a,'b) t -> ('a,'b) t
 val div : ('a,'b) t -> ('a,'b) t -> ('a,'b) t
 (** Divides two matrices element-wise *)
   
+val amax : ('a,'b) t -> float
+(** Maximum of the absolute values of the elements of the matrix. *)
+  
 val add_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit
 (** [add_inplace c a b] : performs [c = a+b] in-place. *)
   
@@ -99,6 +105,12 @@ val of_col_vecs : 'a Vector.t array -> ('a,'b) t
 val of_col_vecs_list : 'a Vector.t list -> ('a,'b) t
 (** Converts a list of vectors  into a matrix *)
 
+val to_array : ('a,'b) t -> float array array
+(** Converts the matrix into an array of arrays *)
+
+val of_array : float array array -> ('a,'b) t
+(** Converts an array of arrays into a matrix *)
+
 val copy: ?m:int -> ?n:int -> ?br:int -> ?bc:int -> ?ar:int -> ?ac:int -> ('a,'b) t -> ('a,'b) t
 (** Copies all or part of a two-dimensional matrix A to a new matrix B *)
 
@@ -228,7 +240,7 @@ val gemm_nn: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
 (** Performs gemm with [~transa=`N] and [~transb=`N]. *)
 
 val gemm_nt: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
-  ?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('c,'b) t -> ('c,'a) t 
+  ?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('c,'b) t -> ('a,'c) t 
 (** Performs gemm with [~transa=`N] and [~transb=`T]. *)
 
 val gemm_tn: ?m:int -> ?n:int -> ?k:int -> ?beta:float -> 
@@ -236,7 +248,7 @@ val gemm_tn: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
 (** Performs gemm with [~transa=`T] and [~transb=`N]. *)
 
 val gemm_tt: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
-  ?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('c,'b) t -> ('c,'b) t 
+  ?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('c,'b) t -> ('a,'c) t 
 (** Performs gemm with [~transa=`T] and [~transb=`T]. *)
 
 
@@ -285,13 +297,13 @@ val x_o_xt : o:('b,'b) t -> x:('a,'b) t -> ('a,'a) t
 val debug_matrix: string -> ('a,'b) t -> unit                                                     
 (** Prints a matrix in stdout for debug *)
 
-val matrix_of_file : string -> ('a,'b) t
+val of_file : string -> ('a,'b) t
 (** Reads a matrix from a file with format "%d %d %f" corresponding to
     [i, j, A.{i,j}]. *)
 
 val relabel : ('a,'b) t -> ('c,'d) t
 
-val sym_matrix_of_file : string -> ('a,'b) t
+val sym_of_file : string -> ('a,'b) t
 (** Reads a symmetric matrix from a file with format "%d %d %f" corresponding to
     [i, j, A.{i,j}]. *)
 

linear_algebra/lib/vector.ml

@@ -32,6 +32,7 @@ let sub t1 t2 = Vec.sub t1 t2
 let mul t1 t2 = Vec.mul t1 t2
 let div t1 t2 = Vec.div t1 t2
 let dot t1 t2 = dot t1 t2
+let amax    t = amax t
 
 let create   n = Vec.create   n
 let make0    n = Vec.make0    n

linear_algebra/lib/vector.mli

@@ -67,6 +67,9 @@ val asin  : 'a t -> 'a t
 val atan  : 'a t -> 'a t
 (** [atan  t = map (f x -> atan x) t *)
   
+val amax : 'a t -> float
+(** Maximum of the absolute values of the elements of the vector. *)
+
 val normalize : 'a t -> 'a t
 (** Returns the vector normalized *)
 

linear_algebra/test/vector.ml

@@ -41,8 +41,8 @@ let test_all () =
   check (float 1.e-14) "norm" (sqrt (dot u2 u2)) (Vector.norm v2);
   check (float 1.e-14) "sum" (Vec.sum u1) (Vector.sum v1);
   check (float 1.e-14) "sum" (Vec.sum u2) (Vector.sum v2);
-  check (float 1.e-14) "at" (u1.{n/2}) (Vector.at v1 (n/2));
+  check (float 1.e-14) "at" (u1.{n/2}) (v1%.(n/2));
-  check (float 1.e-14) "at" (u2.{n/2}) (Vector.at v2 (n/2));
+  check (float 1.e-14) "at" (u2.{n/2}) (v2%.(n/2));
   check (bool) "of_list" true (v1 = Vector.of_list @@ Array.to_list a1);
   check (bool) "of_list" true (v2 = Vector.of_list @@ Array.to_list a2);
   check (bool) "to_list" true (Vector.to_list v1 = Array.to_list a1);

mo/lib/basis.ml

@@ -1,12 +1,13 @@
 open Linear_algebra
-open Common.Util
-open Common.Constants
 
 (** One-electron orthogonal basis set, corresponding to Molecular Orbitals. *)
 
-module HF = HartreeFock
+module HF = Hartree_fock
 module Si = Simulation
 
+type ao = Ao.Ao_dim.t
+type mo = Mo_dim.t
+
 type mo_type =
   | RHF | ROHF | UHF | CASSCF | Projected
   | Natural of string
@@ -16,12 +17,12 @@ type t =
   {
     simulation      : Simulation.t;      (* Simulation which produced the MOs *)
     mo_type         : mo_type;           (* Kind of MOs (RHF, CASSCF, Localized...) *)
-    mo_occupation   : t Vector.t;        (* Occupation numbers *)
+    mo_occupation   : mo Vector.t;       (* Occupation numbers *)
-    mo_coef         : (Ao.Basis.t,t) Matrix.t;       (* Matrix of the MO coefficients in the AO basis *)
+    mo_coef         : (ao,mo) Matrix.t;  (* Matrix of the MO coefficients in the AO basis *)
-    eN_ints         : (t,t) Matrix.t lazy_t;         (* Electron-nucleus potential integrals *)
+    eN_ints         : (mo,mo) Matrix.t lazy_t;      (* Electron-nucleus potential integrals *)
-    ee_ints         : t Four_idx_storage.t lazy_t; (* Electron-electron potential integrals *)
+    ee_ints         : mo Four_idx_storage.t lazy_t; (* Electron-electron potential integrals *)
-    kin_ints        : (t,t) Matrix.t lazy_t;      (* Kinetic energy integrals *)
+    kin_ints        : (mo,mo) Matrix.t lazy_t;      (* Kinetic energy integrals *)
-    one_e_ints      : (t,t) Matrix.t lazy_t;      (* One-electron integrals *)
+    one_e_ints      : (mo,mo) Matrix.t lazy_t;      (* One-electron integrals *)
     (* TODO 
     f12_ints        : F12.t lazy_t;      (* F12 integrals *)
     *)
@@ -50,16 +51,16 @@ let mo_energies t =
   let m_C = mo_coef t in
   let f =
     let m_N = Matrix.of_diag @@ mo_occupation t in
-    let m_P = Matrix.x_o_xt m_N m_C in
+    let m_P = Matrix.x_o_xt ~o:m_N ~x:m_C in
     match t.mo_type with
     | RHF  -> Fock.make_rhf  ~density:m_P (ao_basis t)
     | Projected
-    | ROHF -> (Matrix.scal 0.5 m_P;
+    | ROHF -> (Matrix.scale_inplace 0.5 m_P;
                Fock.make_uhf  ~density_same:m_P ~density_other:m_P (ao_basis t))
     | _ -> failwith "Not implemented"
   in
   let m_F0 = Fock.fock f in
-  Matrix.xt_o_x m_F0 m_C
+  Matrix.xt_o_x ~o:m_F0 ~x:m_C
   |> Matrix.diag
 
 
@@ -86,7 +87,7 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef ()  =
   )
   and ee_ints = lazy (
     Ao.Basis.ee_ints ao_basis
-    |> Eri.four_index_transform mo_coef
+    |> Four_idx_storage.four_index_transform mo_coef
   )
 (*
   and f12_ints = lazy (
@@ -106,7 +107,7 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef ()  =
 let values t point =
   let c = mo_coef t in
   let a = Ao.Basis.values (Simulation.ao_basis t.simulation) point in
-  Matrix.gemv ~trans:`T c a
+  Matrix.gemv_t c a
 
 let of_hartree_fock hf =
   let mo_coef       = HF.eigenvectors hf in
@@ -114,24 +115,24 @@ let of_hartree_fock hf =
   let mo_occupation = HF.occupation   hf in
   let mo_type =
     match HF.kind hf with
-    | HartreeFock.RHF  -> RHF
+    | HF.RHF  -> RHF
-    | HartreeFock.ROHF -> ROHF
+    | HF.ROHF -> ROHF
-    | HartreeFock.UHF  -> UHF
+    | HF.UHF  -> UHF
   in
   make ~simulation ~mo_type ~mo_occupation ~mo_coef ()
 
-
+(*
 let of_mo_basis simulation other =
 
   let mo_coef =
     let basis = Simulation.ao_basis simulation in
     let basis_other = ao_basis other in
     let m_S =
-      Overlap.(matrix @@ of_basis_pair
+      Ao.Overlap.(of_basis_pair
-        (AOBasis.basis basis)
+        (Ao.Basis.ao_basis basis)
-        (AOBasis.basis basis_other) )
+        (Ao.Basis.ao_basis basis_other) )
     in
-    let m_X = AOBasis.ortho basis in
+    let m_X = Ao.Basis.ortho basis in
     (* Project other vectors in the current basis *)
     let m_C =
       gemm m_S @@ mo_coef other
@@ -156,15 +157,15 @@ let of_mo_basis simulation other =
       else 0.)
   in
   make ~simulation ~mo_type:Projected ~mo_occupation ~mo_coef ()
+*)
 
 
 
 
 
 let pp ?(start=1) ?(finish=0) ppf t =
-    let open Lacaml.Io in
+    let rows = Matrix.dim1 t.mo_coef
-    let rows = Mat.dim1 t.mo_coef
+    and cols = Matrix.dim2 t.mo_coef
-    and cols = Mat.dim2 t.mo_coef
     in
     let finish =
       match finish with
@@ -182,7 +183,7 @@ let pp ?(start=1) ?(finish=0) ppf t =
           Array.iteri (fun i x ->
               if (i+1 >= first) && (i+1 <= first+4 ) then
                 Format.fprintf ppf "%12f@ " x)
-            (Vec.to_array @@ mo_energies t);
+            (Vector.to_array @@ mo_energies t);
 
           Format.fprintf ppf "@]@;";
           Format.fprintf ppf "@[%a@]"
@@ -193,7 +194,8 @@ let pp ?(start=1) ?(finish=0) ppf t =
                 (Array.init (min 5 (cols-first+1)) (fun i -> Printf.sprintf "-- %d --" (i + first) ))
               ~print_right:false
               ~print_foot:false
-              () ) (lacpy ~ac:first ~n:(min 5 (cols-first+1)) (t.mo_coef)) ;
+              () ) (Matrix.to_bigarray_inplace t.mo_coef
+                    |> Lacaml.D.lacpy ~ac:first ~n:(min 5 (cols-first+1)) ) ;
           Format.fprintf ppf "@]@;@;@]";
           (aux [@tailcall]) (first+5)
         end

mo/lib/basis.mli

@@ -13,6 +13,8 @@ type mo_type =
   | Localized of string
       
 type t
+type mo = Mo_dim.t
+type ao = Ao.Ao_dim.t
 
 (** {1 Accessors} *)
 
@@ -25,25 +27,25 @@ val mo_type : t -> mo_type
 val ao_basis : t -> Ao.Basis.t
 (** Matrix of the MO coefficients in the AO basis *)
 
-val mo_occupation : t -> t Vector.t
+val mo_occupation : t -> mo Vector.t
 (** Occupation numbers *)
 
-val mo_coef : t -> (Ao.Basis.t, t) Matrix.t
+val mo_coef : t -> (ao, mo) Matrix.t
 (** Molecular orbitcal coefficients *)
 
-val eN_ints : t -> (t,t) Matrix.t
+val eN_ints : t -> (mo,mo) Matrix.t
 (** Electron-nucleus potential integrals *)
 
-val ee_ints : t -> t Four_idx_storage.t
+val ee_ints : t -> mo Four_idx_storage.t
 (** Electron-electron repulsion integrals *)
 
-val kin_ints : t -> (t,t) Matrix.t
+val kin_ints : t -> (mo,mo) Matrix.t
 (** Kinetic energy integrals *)
 
-val one_e_ints : t -> (t,t) Matrix.t
+val one_e_ints : t -> (mo,mo) Matrix.t
 (** One-electron integrals {% $\hat{T} + V$ %} *)
 
-val two_e_ints : t -> t Four_idx_storage.t
+val two_e_ints : t -> mo Four_idx_storage.t
 (** Electron-electron repulsion integrals *)
 
 (* TODO
@@ -54,39 +56,38 @@ val f12_ints : t -> F12.t
 val size : t -> int
 (** Number of molecular orbitals in the basis *)
 
-val mo_energies : t -> t Vector.t
+val mo_energies : t -> mo Vector.t
 (** Fock MO energies *)
 
-val values : t -> Coordinate.t -> t Vector.t
+val values : t -> Coordinate.t -> mo Vector.t
 (** Values of the MOs evaluated at a given coordinate. *)
 
 (** {1 Creators} *)
 
 val make : simulation:Simulation.t ->
   mo_type:mo_type ->
-  mo_occupation:t Vector.t ->
+  mo_occupation:mo Vector.t ->
-  mo_coef:(Ao.Basis.t,t) Matrix.t ->
+  mo_coef:(ao,mo) Matrix.t ->
   unit -> t
 (** Function to build a data structure representing the molecular orbitals. *)
 
-val of_hartree_fock : HartreeFock.t -> t
+val of_hartree_fock : Hartree_fock.t -> t
 (** Build MOs from a Restricted Hartree-Fock calculation. *)
 
+(*
 val of_mo_basis : Simulation.t -> t -> t
 (** Project the MOs of the other basis on the current one. *)
+*)
 
 
 val mo_matrix_of_ao_matrix :
-  mo_coef:(Ao.Basis.t,t) Matrix.t ->
+  mo_coef:(ao,mo) Matrix.t ->
-  (Ao.Basis.t,Ao.Basis.t) Matrix.t ->
+  (ao,ao) Matrix.t -> (mo,mo) Matrix.t
-  (t,t) Matrix.t
 (** Build a matrix in MO basis from a matrix in AO basis. *)
 
 val ao_matrix_of_mo_matrix :
-  mo_coef:(Ao.Basis.t,t) Matrix.t ->
+  mo_coef:(ao,mo) Matrix.t -> ao_overlap:(ao,ao) Matrix.t ->
-  ao_overlap:(Ao.Basis.t,Ao.Basis.t) Matrix.t ->
+  (mo,mo) Matrix.t -> (ao,ao) Matrix.t
-  (t,t) Matrix.t ->
-  (Ao.Basis.t,Ao.Basis.t) Matrix.t
 (** Build a matrix in AO basis from a matrix in MO basis. *)
 
 (** {1 Printers} *)

mo/lib/class.ml

@@ -1,4 +1,5 @@
 open Particles
+open Common
     
 type mo_class =                                                                              
   | Core       of int  (* Always doubly occupied *)

mo/lib/fock.ml

@@ -0,0 +1,240 @@
+open Linear_algebra
+open Common
+
+type ao = Ao.Ao_dim.t
+
+type t =
+  {
+    fock     : (ao,ao) Matrix.t ;
+    core     : (ao,ao) Matrix.t ;
+    coulomb  : (ao,ao) Matrix.t ;
+    exchange : (ao,ao) Matrix.t ;
+  }
+
+
+let fock t = t.fock
+let core t = t.core
+let coulomb t = t.coulomb
+let exchange t = t.exchange
+
+
+let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis =
+  let m_P  = density
+  and m_T  = Ao.Basis.kin_ints ao_basis
+  and m_V  = Ao.Basis.eN_ints  ao_basis
+  and m_G  = Ao.Basis.ee_ints  ao_basis 
+  in
+  let nBas = Matrix.dim1 m_T
+  in
+
+  let m_Hc = Matrix.add m_T m_V 
+  and m_J = Array.make_matrix nBas nBas 0.
+  and m_K = Array.make_matrix nBas nBas 0.
+  in
+
+  for sigma = 1 to nBas do
+    let m_Ksigma = m_K.(sigma-1) in
+    for nu = 1 to nBas do
+      let m_Jnu = m_J.(nu-1) in
+      for lambda = 1 to nBas do
+        let pJ = m_P%:(lambda,sigma)
+        and pK = 0.5 *. (m_P%:(lambda,nu))
+        in
+        match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with
+        | (false, false, _) -> ()
+        | (true , true , true) ->
+            begin
+              for mu = 1 to nu do
+                let integral = 
+                  Four_idx_storage.get_phys m_G mu lambda nu sigma
+                in
+                if (integral <> 0.) then begin
+                  m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  pJ *. integral;
+                  m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
+                end
+              done;
+              for mu = nu+1 to sigma do
+                m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +.
+                                   pK *. Four_idx_storage.get_phys m_G mu lambda nu sigma
+              done
+            end
+        | (true , true , false) ->
+            begin
+              for mu = 1 to sigma do
+                let integral = 
+                  Four_idx_storage.get_phys m_G mu lambda nu sigma
+                in
+                if (integral <> 0.) then begin
+                  m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  pJ *. integral;
+                  m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
+                end
+              done;
+              for mu = sigma+1 to nu do
+                m_Jnu.(mu-1) <-
+                  m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma
+              done
+            end
+        | (false, true , _) ->
+            for mu = 1 to sigma do
+              m_Ksigma.(mu-1) <-
+                m_Ksigma.(mu-1) +. pK *. Four_idx_storage.get_phys m_G mu lambda nu sigma
+            done
+        | (true , false, _) ->
+            for mu = 1 to nu do
+              m_Jnu.(mu-1) <-
+                m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma
+            done
+      done
+    done;
+    for mu = 1 to sigma-1 do
+      m_K.(mu-1).(sigma-1) <- m_Ksigma.(mu-1);
+    done
+  done;
+  for nu = 1 to nBas do
+    let m_Jnu = m_J.(nu-1) in
+    for mu = 1 to nu-1 do
+      m_J.(mu-1).(nu-1) <- m_Jnu.(mu-1)
+    done
+  done;
+
+  let m_J = Matrix.of_array m_J
+  and m_K = Matrix.of_array m_K
+  in
+  { fock = Matrix.add m_Hc (Matrix.sub m_J m_K) ;
+    core = m_Hc ; coulomb = m_J ; exchange = m_K }
+
+
+
+let make_uhf ~density_same ~density_other ?(threshold=Constants.epsilon) ao_basis =
+  let m_P_a  = density_same
+  and m_P_b  = density_other
+  and m_T  = Ao.Basis.kin_ints ao_basis
+  and m_V  = Ao.Basis.eN_ints  ao_basis
+  and m_G  = Ao.Basis.ee_ints  ao_basis 
+  in
+  let nBas = Matrix.dim1 m_T
+  in
+
+  let m_Hc = Matrix.add m_T m_V 
+  and m_J = Array.make_matrix nBas nBas 0.
+  and m_K = Array.make_matrix nBas nBas 0.
+  in
+
+  for sigma = 1 to nBas do
+    let m_Ksigma = m_K.(sigma-1) in
+    for nu = 1 to nBas do
+      let m_Jnu = m_J.(nu-1) in
+      for lambda = 1 to nBas do
+        let pJ = m_P_a%:(lambda,sigma) +. m_P_b%:(lambda,sigma)
+        and pK = m_P_a%:(lambda,nu)
+        in
+        match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with
+        | (false, false, _) -> ()
+        | (true , true , true) ->
+            begin
+              for mu = 1 to nu do
+                let integral = 
+                  Four_idx_storage.get_phys m_G mu lambda nu sigma
+                in
+                if (integral <> 0.) then begin
+                  m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  pJ *. integral;
+                  m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
+                end
+              done;
+              for mu = nu+1 to sigma do
+                m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +.
+                                   pK *.Four_idx_storage.get_phys m_G mu lambda nu sigma
+              done
+            end
+        | (true , true , false) ->
+            begin
+              for mu = 1 to sigma do
+                let integral = 
+                  Four_idx_storage.get_phys m_G mu lambda nu sigma
+                in
+                if (integral <> 0.) then begin
+                  m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  pJ *. integral;
+                  m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
+                end
+              done;
+              for mu = sigma+1 to nu do
+                m_Jnu.(mu-1) <-
+                  m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma
+              done
+            end
+        | (false, true , _) ->
+            for mu = 1 to sigma do
+              m_Ksigma.(mu-1) <-
+                m_Ksigma.(mu-1) +. pK *. Four_idx_storage.get_phys m_G mu lambda nu sigma
+            done
+        | (true , false, _) ->
+            for mu = 1 to nu do
+              m_Jnu.(mu-1) <-
+                m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma
+            done
+      done
+    done;
+    for mu = 1 to sigma-1 do
+      m_K.(mu-1).(sigma-1) <- m_Ksigma.(mu-1);
+    done
+  done;
+  for nu = 1 to nBas do
+    let m_Jnu = m_J.(nu-1) in
+    for mu = 1 to nu-1 do
+      m_J.(mu-1).(nu-1) <- m_Jnu.(mu-1)
+    done
+  done;
+  
+  let m_J = Matrix.of_array m_J
+  and m_K = Matrix.of_array m_K
+  in
+  { fock = Matrix.add m_Hc (Matrix.sub m_J m_K) ;
+    core = m_Hc ; coulomb = m_J ; exchange = m_K }
+
+
+
+
+
+
+
+let op ~f f1 f2 =
+  assert (f1.core = f2.core);
+  let m_Hc = f1.core
+  and m_J = f f1.coulomb f2.coulomb
+  and m_K = f f1.exchange f2.exchange
+  in
+  {
+    fock     = Matrix.add m_Hc (Matrix.sub m_J m_K);
+    core     = m_Hc;
+    coulomb  = m_J;
+    exchange = m_K;
+  }
+
+
+let add = op ~f:(fun a b -> Matrix.add a b)
+
+let sub = op ~f:(fun a b -> Matrix.sub a b)
+
+let scale alpha f1 = 
+  let m_Hc = f1.core
+  and m_J = Matrix.copy f1.coulomb 
+  and m_K = Matrix.copy f1.exchange 
+  in
+  Matrix.scale_inplace alpha m_J;
+  Matrix.scale_inplace alpha m_K;
+  {
+    fock     = Matrix.add m_Hc (Matrix.sub m_J m_K);
+    core     = m_Hc;
+    coulomb  = m_J;
+    exchange = m_K;
+  }
+
+
+
+let pp ppf a =
+  Format.fprintf ppf "@[<2>";
+  Format.fprintf ppf "@[ Fock matrix:@[<2>@[%a@]@.]@]" Matrix.pp a.fock;
+  Format.fprintf ppf "@[ Core Hamiltonian:@[<2>@[%a@]@.]@]" Matrix.pp a.core;
+  Format.fprintf ppf "@[ Coulomb matrix:@[<2>@[%a@]@.]@]" Matrix.pp a.coulomb;
+  Format.fprintf ppf "@[ Exchange matrix:@[<2>@[%a@]@.]@]" Matrix.pp a.exchange;
+  Format.fprintf ppf "@]"

mo/lib/fock.mli

@@ -0,0 +1,54 @@
+(** Type for the Fock operator in AO basis. *)
+
+open Linear_algebra
+
+type t 
+type ao = Ao.Ao_dim.t
+
+(** {1 Accessors} *)
+
+val fock : t -> (ao,ao) Matrix.t
+(** Fock matrix in AO basis *)
+
+val core : t -> (ao,ao) Matrix.t
+(** Core Hamiltonian : {% $\langle i | \hat{h} | j \rangle$ %} *)
+
+val coulomb : t -> (ao,ao) Matrix.t
+(** Coulomb matrix : {% $\langle i | J | j \rangle$ %} *)
+
+val exchange : t -> (ao,ao) Matrix.t
+(** Exchange matrix : {% $\langle i | K | j \rangle$ %} *)
+
+
+(** {1 Creators} *)
+
+val make_rhf : density:(ao,ao) Matrix.t -> ?threshold:float -> Ao.Basis.t -> t
+(** Create a Fock operator in the RHF formalism. Expected density is
+    {% $2 \mathbf{C\, C}^\dagger$ %}. [threshold] is a threshold on the
+    integrals. *)
+
+val make_uhf :
+  density_same:(ao,ao) Matrix.t ->
+  density_other:(ao,ao) Matrix.t ->
+  ?threshold:float ->
+  Ao.Basis.t -> t
+(** Create a Fock operator in the UHF formalism. Expected density is
+    {% $\mathbf{C\, C}^\dagger$ %}. When building the {% $\alpha$ %} Fock
+    operator, [density_same] is the {% $\alpha$ %} density and [density_other]
+    is the {% $\beta$ %} density. [threshold] is a threshold on the integrals. *)
+
+
+(** {1 Operations} *)
+
+val add : t -> t -> t
+(** Add two Fock operators sharing the same core Hamiltonian. *)
+
+val sub : t -> t -> t
+(** Subtract two Fock operators sharing the same core Hamiltonian. *)
+
+val scale : float -> t -> t
+(** Scale the Fock matrix by a constant *)
+
+(** {1 Printers} *)
+
+val pp : Format.formatter -> t -> unit

mo/lib/guess.ml

@@ -0,0 +1,56 @@
+open Linear_algebra
+
+type ao = Ao.Ao_dim.t
+type mo = Mo_dim.t
+
+type guess =
+| Hcore  of (ao,ao) Matrix.t
+| Huckel of (ao,ao) Matrix.t
+| Matrix of (ao,mo) Matrix.t
+
+type t = guess
+
+module El = Particles.Electrons
+
+let hcore_guess ao_basis = 
+  let eN_ints  = Ao.Basis.eN_ints  ao_basis
+  and kin_ints = Ao.Basis.kin_ints ao_basis
+  in
+  Matrix.add eN_ints kin_ints
+
+
+let huckel_guess ao_basis =
+  let c = 0.5 *. 1.75 in
+  let eN_ints  = Ao.Basis.eN_ints  ao_basis
+  and kin_ints = Ao.Basis.kin_ints ao_basis
+  in
+  let m_F = 
+    Matrix.add eN_ints kin_ints
+  in
+  let ao_num   = Ao.Basis.size     ao_basis
+  and overlap  = Ao.Basis.overlap  ao_basis
+  in
+  let diag = Vector.init ao_num (fun i -> m_F%:(i,i) ) in
+
+  function 
+  | 0 -> invalid_arg "Huckel guess needs a non-zero number of occupied MOs."
+  | _nocc -> 
+    Matrix.init_cols ao_num ao_num (fun i j ->
+        if (i<>j) then
+          if (diag%.(i) +. diag%.(j)) < 0. then
+            c *. (overlap%:(i,j)) *. (diag%.(i) +. diag%.(j)) +. m_F%:(i,j) (*TODO Pseudo *)
+          else
+            m_F%:(i,j) (*TODO Pseudo *)
+        else
+          diag%.(i)
+    )
+
+
+let make ?(nocc=0) ~guess ao_basis = 
+  match guess with
+  | `Hcore  -> Hcore  (hcore_guess  ao_basis)
+  | `Huckel -> Huckel (huckel_guess ao_basis nocc)
+  | `Matrix m -> Matrix m
+
+
+

mo/lib/guess.mli

@@ -0,0 +1,20 @@
+open Linear_algebra
+
+(** Guess for Hartree-Fock calculations. *)
+
+type ao = Ao.Ao_dim.t
+type mo = Mo_dim.t
+
+type guess =
+| Hcore  of (ao,ao) Matrix.t (* Core Hamiltonian Matrix *)
+| Huckel of (ao,ao) Matrix.t (* Huckel Hamiltonian Matrix *)
+| Matrix of (ao,mo) Matrix.t (* Guess Eigenvectors *)
+
+type t = guess
+
+
+val make :
+  ?nocc:int ->
+  guess:[ `Hcore | `Huckel | `Matrix of (ao,mo) Matrix.t ] ->
+  Ao.Basis.t -> t
+

mo/lib/hartree_fock.ml

@@ -0,0 +1,721 @@
+open Linear_algebra
+open Particles
+open Common
+
+type ao = Ao.Ao_dim.t
+type mo = Mo_dim.t
+
+type hartree_fock_data =
+  {
+    iteration       :  int              ;
+    coefficients    :  (ao, mo) Matrix.t   option  ;
+    eigenvalues     :  mo Vector.t option   ;
+    error           :  float       option   ;
+    diis            :  (mo*mo) Diis.t   option     ;
+    energy          :  float       option   ;
+    density         :  (ao,ao) Matrix.t   option   ;
+    density_a       :  (ao,ao) Matrix.t   option   ;
+    density_b       :  (ao,ao) Matrix.t   option   ;
+    fock            :  Fock.t  option   ;
+    fock_a          :  Fock.t  option   ;
+    fock_b          :  Fock.t  option   ;
+  }
+
+and hartree_fock_kind  =
+| RHF  (** Restricted Hartree-Fock *)
+| ROHF (** Restricted Open-shell Hartree-Fock *)
+| UHF  (** Unrestricted Hartree-Fock *)
+
+and t =
+  {
+    kind              : hartree_fock_kind;
+    simulation        : Simulation.t;
+    guess             : Guess.t;
+    data              : hartree_fock_data option lazy_t array;
+    nocc              : int   ;
+  }
+
+
+let empty = 
+  { 
+    iteration       =  0                ; 
+    coefficients    =  None             ; 
+    eigenvalues     =  None             ; 
+    error           =  None             ; 
+    diis            =  None             ; 
+    energy          =  None             ; 
+    density         =  None             ; 
+    density_a       =  None             ; 
+    density_b       =  None             ; 
+    fock            =  None             ; 
+    fock_a          =  None             ; 
+    fock_b          =  None             ; 
+  } 
+
+
+
+
+
+
+module Si = Simulation
+module El = Electrons
+
+
+let  kind               t  =  t.kind
+let  simulation         t  =  t.simulation
+let  guess              t  =  t.guess
+let  nocc               t  =  t.nocc
+
+
+
+let n_iterations t =
+  Array.fold_left (fun accu x ->
+    match Lazy.force x with
+    | Some _ -> accu + 1
+    | None -> accu
+    ) 0 t.data
+
+
+let last_iteration t =
+  Util.of_some @@ Lazy.force (t.data.(n_iterations t - 1))
+
+let eigenvectors t =
+  let data = last_iteration t in
+  Util.of_some data.coefficients
+
+let eigenvalues t =
+  let data = last_iteration t in
+  Util.of_some data.eigenvalues
+
+let density t =
+  let data = last_iteration t in
+  match kind t with
+  | RHF  -> Util.of_some data.density
+  | ROHF -> Matrix.add (Util.of_some data.density_a) (Util.of_some data.density_b)
+  | _    -> failwith "Not implemented"
+
+let occupation t =
+  let n_alfa, n_beta =
+    El.n_alfa @@ Simulation.electrons @@ simulation t,
+    El.n_beta @@ Simulation.electrons @@ simulation t
+  in
+  match kind t with
+  | RHF  -> Vector.init (Matrix.dim2 @@ eigenvectors t) (fun i ->
+                if i <= nocc t then 2.0 else 0.0)
+  | ROHF ->  Vector.init (Matrix.dim2 @@ eigenvectors t) (fun i ->
+                 if i <= n_beta then 2.0 else
+                 if i <= n_alfa then 1.0 else
+                 0.0)
+  | _    -> failwith "Not implemented"
+
+
+let energy t = 
+  let data = last_iteration t in
+  Util.of_some data.energy
+
+
+let nuclear_repulsion t =
+  Si.nuclear_repulsion (simulation t)
+
+
+let ao_basis t =
+  Si.ao_basis (simulation t)
+
+
+let kin_energy t =
+  let m_T    =
+    ao_basis t
+    |> Ao.Basis.kin_ints 
+  in
+  let m_P = density t in
+  Matrix.gemm_trace m_P m_T
+
+
+let eN_energy t =
+  let m_V =
+    ao_basis t
+    |> Ao.Basis.eN_ints 
+  in
+  let m_P = density t in
+  Matrix.gemm_trace m_P m_V
+
+
+let coulomb_energy t =
+  let data =
+    last_iteration t
+  in
+  match kind t with
+  | RHF ->  let  m_P   =  Util.of_some    data.density  in
+            let  fock  =  Util.of_some    data.fock     in
+            let  m_J   =  Fock.coulomb    fock          in
+            0.5 *. Matrix.gemm_trace m_P m_J
+
+  | ROHF -> let  m_P_a  =  Util.of_some    data.density_a in
+            let  m_P_b  =  Util.of_some    data.density_b in
+            let  fock_a =  Util.of_some    data.fock_a    in
+            let  fock_b =  Util.of_some    data.fock_b    in
+            let  m_J_a  =  Fock.coulomb    fock_a         in
+            let  m_J_b  =  Fock.coulomb    fock_b         in
+            0.5 *. ( (Matrix.gemm_trace m_P_a m_J_a) +. (Matrix.gemm_trace m_P_b m_J_b) )
+
+  | _    -> failwith "Not implemented"
+
+
+let exchange_energy t = 
+  let data =
+    last_iteration t
+  in
+  match kind t with
+  | RHF ->  let  m_P   =  Util.of_some  data.density  in
+            let  fock  =  Util.of_some  data.fock     in
+            let  m_K   =  Fock.exchange fock          in
+            -. 0.5 *. Matrix.gemm_trace m_P m_K
+
+  | ROHF -> let  m_P_a  =  Util.of_some   data.density_a in
+            let  m_P_b  =  Util.of_some   data.density_b in
+            let  fock_a =  Util.of_some   data.fock_a    in
+            let  fock_b =  Util.of_some   data.fock_b    in
+            let  m_K_a  =  Fock.exchange  fock_a         in
+            let  m_K_b  =  Fock.exchange  fock_b         in
+            -. 0.5 *. ( (Matrix.gemm_trace m_P_a m_K_a) +. (Matrix.gemm_trace m_P_b m_K_b) )
+
+  | _    -> failwith "Not implemented"
+
+
+
+
+
+
+let make
+    ?kind
+    ?guess:(guess=`Huckel)
+    ?max_scf:(max_scf=64)
+    ?level_shift:(level_shift=0.2)
+    ?threshold_SCF:(threshold_SCF=1.e-8)
+    simulation = 
+
+
+  (* Number of occupied MOs *)
+  let n_alfa, n_beta =
+    El.n_alfa @@ Si.electrons simulation,
+    El.n_beta @@ Si.electrons simulation
+  in
+
+  let nocc = n_alfa in
+
+  let kind = 
+    match kind with
+    | Some kind -> kind
+    | None -> if (n_alfa = n_beta) then RHF else ROHF
+  in
+
+  let nuclear_repulsion = 
+    Si.nuclear_repulsion simulation
+  in
+
+  let ao_basis = 
+    Si.ao_basis simulation
+  in
+
+
+  (* Orthogonalization matrix *)
+  let m_X = 
+    Ao.Basis.ortho ao_basis
+  in
+
+  (* Overlap matrix *)
+  let m_S =
+    Ao.Basis.overlap ao_basis
+  in
+
+  (* Level shift in MO basis *)
+  let m_LSmo =
+    Array.init (Matrix.dim2 m_X) (fun i ->
+        if i > nocc then level_shift else 0.)
+    |> Vector.of_array
+    |> Matrix.of_diag
+  in
+
+  (* Guess coefficients *)
+  let guess =
+    Guess.make ~nocc ~guess ao_basis
+  in
+
+  let m_C : (ao,mo) Matrix.t = 
+    let c_of_h m_H = 
+      let m_Hmo   = Matrix.xt_o_x  ~o:m_H  ~x:m_X in
+      let m_C', _ = Matrix.diagonalize_symm  m_Hmo in
+      Matrix.gemm_nn  m_X  m_C'
+    in
+    match guess with
+      | Guess.Hcore  m_H -> c_of_h m_H
+      | Guess.Huckel m_H -> c_of_h m_H
+      | Guess.Matrix m_C -> m_C 
+  in
+
+  (* A single SCF iteration *)
+  let scf_iteration_rhf data =
+
+    let  nSCF         =  data.iteration + 1
+    and  m_C          =  Util.of_some data.coefficients
+    and  m_P_prev     =  data.density
+    and  fock_prev    =  data.fock
+    and  diis         =  
+      match data.diis with
+      | Some diis -> diis
+      | None      -> Diis.make ()
+    and  threshold    =  
+      match data.error with
+      | Some error -> error
+      | None       -> threshold_SCF *. 2.
+    in
+
+    (* Density matrix over nocc occupied MOs *)
+    let m_P =
+      Matrix.gemm_nt  ~alpha:2.  ~k:nocc  m_C  m_C
+    in
+
+    (* Fock matrix in AO basis *)
+    let fock = 
+      match fock_prev, m_P_prev, threshold > 100. *. threshold_SCF with
+      | Some fock_prev, Some m_P_prev, true -> 
+        let threshold = 1.e-8 in
+        Fock.make_rhf  ~density:(Matrix.sub m_P m_P_prev) ~threshold ao_basis
+        |> Fock.add fock_prev
+      | _ -> Fock.make_rhf  ~density:m_P ao_basis
+    in
+
+    let m_F0, m_Hc, _, _ =
+      let x = fock in
+      Fock.(fock x, core x, coulomb x, exchange x)
+    in
+
+    (* Add level shift in AO basis *)
+    let m_F =
+      let m_SC =
+        Matrix.gemm_nn m_S m_C
+      in
+      Matrix.gemm_nn m_SC (Matrix.gemm_nt m_LSmo m_SC)
+      |> Matrix.add m_F0
+    in
+
+
+    (* Fock matrix in orthogonal basis *)
+    let m_F_ortho =
+      Matrix.xt_o_x ~o:m_F ~x:m_X
+    in
+
+    let error_fock = 
+      let fps =
+        Matrix.gemm_nn  m_F  (Matrix.gemm_nn  m_P  m_S)
+      and spf =
+        Matrix.gemm_nn  m_S  (Matrix.gemm_nn  m_P  m_F)
+      in
+      Matrix.xt_o_x  ~o:(Matrix.sub  fps  spf)  ~x:m_X
+    in
+
+    let diis, m_F_diis = 
+      let diis =
+        Diis.append
+          ~p:(Matrix.as_vec_inplace m_F_ortho)
+          ~e:(Matrix.as_vec_inplace error_fock) diis
+      in
+
+      try
+        let m_F_diis =
+          Diis.next diis
+          |> Matrix.reshape_vec_inplace (Matrix.dim1 m_F_ortho) (Matrix.dim2 m_F_ortho)
+        in
+        diis, m_F_diis
+
+      with Failure _ -> (* Failure in DIIS.next *)
+        Diis.make (), m_F_ortho
+    in
+    let diis = 
+      Diis.append
+        ~p:(Matrix.as_vec_inplace m_F_ortho)
+        ~e:(Matrix.as_vec_inplace error_fock) diis
+    in
+
+
+    (* MOs in orthogonal MO basis *)
+    let m_C', _ =
+      Matrix.diagonalize_symm  m_F_diis
+    in
+
+    (* Re-compute eigenvalues to remove level-shift *)
+    let eigenvalues = 
+      let m_F_ortho =
+        Matrix.xt_o_x ~o:m_F0 ~x:m_X
+      in
+      Matrix.xt_o_x ~o:m_F_ortho ~x:m_C'
+      |> Matrix.diag 
+    in
+
+    (* MOs in AO basis *)
+    let m_C =
+      Matrix.gemm_nn  m_X  m_C'
+      |> Conventions.rephase
+    in
+
+    (* Hartree-Fock energy *)
+    let energy =
+      nuclear_repulsion +. 0.5 *.
+                           Matrix.gemm_trace  m_P  (Matrix.add  m_Hc  m_F)
+    in
+
+    (* Convergence criterion *)
+    let error =
+      error_fock 
+      |> Matrix.amax
+      |> abs_float
+    in
+
+    { empty with
+      iteration      = nSCF              ;
+      eigenvalues    = Some eigenvalues  ;
+      coefficients   = Some m_C          ;
+      error          = Some error        ;
+      diis           = Some diis         ;
+      energy         = Some energy       ;
+      density        = Some m_P          ;
+      fock           = Some fock         ;
+    }
+
+  in
+
+  let scf_iteration_rohf data = 
+
+    let  nSCF         =  data.iteration + 1
+    and  m_C          =  Util.of_some data.coefficients
+    and  m_P_a_prev   =  data.density_a
+    and  m_P_b_prev   =  data.density_b
+    and  fock_a_prev  =  data.fock_a
+    and  fock_b_prev  =  data.fock_b
+    and  diis         =  
+      match data.diis with
+      | Some diis -> diis
+      | None      -> Diis.make ()
+    and  threshold    =  
+      match data.error with
+      | Some error -> error
+      | None       -> threshold_SCF *. 2.
+    in
+    
+    (* Density matrix *)
+    let m_P_a =
+      Matrix.gemm_nt  ~alpha:1.  ~k:n_alfa  m_C  m_C
+    in
+
+    let m_P_b =
+      Matrix.gemm_nt  ~alpha:1.  ~k:n_beta  m_C  m_C
+    in
+
+    let m_P =
+      Matrix.add m_P_a m_P_b
+    in
+
+    (* Fock matrix in AO basis *)
+    let fock_a = 
+      match fock_a_prev, threshold > 100. *. threshold_SCF with
+      | Some fock_a_prev, true -> 
+        let threshold = 1.e-8 in
+        Fock.make_uhf  ~density_same:(Matrix.sub m_P_a @@ Util.of_some m_P_a_prev) ~density_other:(Matrix.sub m_P_b @@ Util.of_some m_P_b_prev) ~threshold ao_basis
+        |> Fock.add fock_a_prev
+      | _ -> Fock.make_uhf  ~density_same:m_P_a  ~density_other:m_P_b ao_basis
+    in
+
+    let fock_b = 
+      match fock_b_prev, threshold > 100. *. threshold_SCF with
+      | Some fock_b_prev, true -> 
+        let threshold = 1.e-8 in
+        Fock.make_uhf  ~density_same:(Matrix.sub m_P_b @@ Util.of_some m_P_b_prev) ~density_other:(Matrix.sub m_P_a @@ Util.of_some m_P_a_prev) ~threshold ao_basis
+        |> Fock.add fock_b_prev
+      | _ -> Fock.make_uhf  ~density_same:m_P_b  ~density_other:m_P_a ao_basis
+    in
+
+    let m_F_a, m_Hc_a, _, _ =
+      let x = fock_a in
+      Fock.(fock x, core x, coulomb x, exchange x)
+    in
+
+    let m_F_b, m_Hc_b, _, _ =
+      let x = fock_b in
+      Fock.(fock x, core x, coulomb x, exchange x)
+    in
+
+
+    let m_F_mo_a = 
+      Matrix.xt_o_x ~o:m_F_a ~x:m_C
+    in
+
+    let m_F_mo_b = 
+      Matrix.xt_o_x ~o:m_F_b ~x:m_C
+    in
+
+    let m_F_mo  =
+      let space k =
+        if k <= n_beta then
+          `Core
+        else if k <= n_alfa then
+          `Active
+        else
+          `Virtual
+      in
+      Array.init (Matrix.dim2 m_F_mo_a) (fun i ->
+          let i = i+1 in
+          Array.init (Matrix.dim1 m_F_mo_a) (fun j ->
+              let j = j+1 in
+              match (space i), (space j) with
+              |  `Core     ,  `Core     ->  
+                0.5 *. (m_F_mo_a%:(i,j) +. m_F_mo_b%:(i,j)) -.
+                (m_F_mo_b%:(i,j) -. m_F_mo_a%:(i,j))
+
+              |  `Active   ,  `Core
+              |  `Core     ,  `Active   ->
+                m_F_mo_b%:(i,j)
+
+              |  `Core     ,  `Virtual
+              |  `Virtual  ,  `Core     
+              |  `Active   ,  `Active   ->
+                0.5 *. (m_F_mo_a%:(i,j) +. m_F_mo_b%:(i,j))
+
+              |  `Virtual  ,  `Active   
+              |  `Active   ,  `Virtual  ->
+                m_F_mo_a%:(i,j)
+
+              |  `Virtual  ,  `Virtual  ->
+                0.5 *. (m_F_mo_a%:(i,j) +. m_F_mo_b%:(i,j)) +.
+                (m_F_mo_b%:(i,j) -. m_F_mo_a%:(i,j))
+            ) )
+      |> Matrix.of_array
+    in
+
+    let m_SC =
+      Matrix.gemm_nn  m_S  m_C
+    in
+
+    let m_F0 = 
+      Matrix.x_o_xt ~x:m_SC ~o:m_F_mo
+    in
+
+
+    (* Add level shift in AO basis *)
+    let m_F =
+      Matrix.x_o_xt ~x:m_SC ~o:m_LSmo
+      |> Matrix.add m_F0
+    in
+
+    (* Fock matrix in orthogonal basis *)
+    let m_F_ortho =
+      Matrix.xt_o_x ~o:m_F ~x:m_X
+    in
+
+    let error_fock = 
+      let fps =
+        Matrix.gemm_nn  m_F  (Matrix.gemm_nn  m_P  m_S)
+      and spf =
+        Matrix.gemm_nn  m_S  (Matrix.gemm_nn  m_P  m_F)
+      in
+      Matrix.xt_o_x ~o:(Matrix.sub  fps  spf) ~x:m_X
+    in
+
+    let diis, m_F_diis = 
+      let diis =
+        Diis.append
+          ~p:(Matrix.as_vec_inplace m_F_ortho)
+          ~e:(Matrix.as_vec_inplace error_fock) diis
+      in
+
+      try
+        let m_F_diis =
+          Diis.next diis
+          |> Matrix.reshape_vec_inplace (Matrix.dim1 m_F_ortho) (Matrix.dim2 m_F_ortho)
+        in
+        diis, m_F_diis
+
+      with Failure _ -> (* Failure in DIIS.next *)
+        Diis.make (), m_F_ortho
+    in
+
+
+    (* MOs in orthogonal MO basis *)
+    let m_C', _ =
+      Matrix.diagonalize_symm  m_F_diis
+    in
+
+    (* Re-compute eigenvalues to remove level-shift *)
+    let eigenvalues = 
+      let m_F_ortho =
+        Matrix.xt_o_x ~o:m_F0 ~x:m_X
+      in
+      Matrix.xt_o_x ~o:m_F_ortho ~x:m_C'
+      |> Matrix.diag 
+    in
+
+    (* MOs in AO basis *)
+    let m_C =
+      Matrix.gemm_nn  m_X  m_C'
+      |> Conventions.rephase
+    in
+
+    (* Hartree-Fock energy *)
+    let energy =
+      nuclear_repulsion +.  0.5 *. ( Matrix.gemm_trace  m_P_a  (Matrix.add  m_Hc_a  m_F_a) +.
+                                     Matrix.gemm_trace  m_P_b  (Matrix.add  m_Hc_b  m_F_b) )
+    in
+
+    (* Convergence criterion *)
+    let error =
+      error_fock 
+      |> Matrix.amax
+      |> abs_float
+    in
+    { empty with
+      iteration      = nSCF              ;
+      eigenvalues    = Some eigenvalues  ;
+      coefficients   = Some m_C          ;
+      error          = Some error        ;
+      diis           = Some diis         ;
+      energy         = Some energy       ;
+      density_a      = Some m_P_a        ;
+      density_b      = Some m_P_b        ;
+      fock_a         = Some fock_a       ;
+      fock_b         = Some fock_b       ;
+    }
+
+  in
+
+
+  let scf_iteration data =
+    match kind with
+    | RHF  -> scf_iteration_rhf data
+    | ROHF -> scf_iteration_rohf data
+    | _    -> failwith "Not implemented"
+  in
+
+  let array_data =
+      
+    let storage =
+      Array.make max_scf None
+    in
+
+    let rec iteration = function
+      | 0 -> Some (scf_iteration { empty with coefficients = Some m_C })
+      | n -> begin
+          match storage.(n) with
+          | Some result -> Some result
+          | None ->
+            begin
+              let data = iteration (n-1) in
+              match data with
+              | None -> None
+              | Some data -> 
+                begin
+                  (* Check convergence *)
+                  let converged, _error =
+                    match data.error with
+                    | None -> false, 0.
+                    | Some error -> (data.iteration = max_scf || error < threshold_SCF), error
+                  in
+                  if converged then
+                    None
+                  else
+                    begin
+                      storage.(n-1) <- Some { empty with
+                                              iteration   = data.iteration;
+                                              energy      = data.energy ;
+                                              eigenvalues = data.eigenvalues ;
+                                              error       = data.error ;
+                                            };
+                      storage.(n) <- Some (scf_iteration data);
+                      storage.(n);
+                    end
+                end
+            end
+        end
+    in
+    Array.init max_scf (fun i -> lazy (iteration i))
+  in
+  { 
+      kind;
+      simulation;
+      guess ;
+      data = array_data;
+      nocc;
+  }
+
+
+
+
+
+
+
+
+
+
+let linewidth = 60
+
+let pp_iterations ppf t =
+  let line = (String.make linewidth '-') in
+  Format.fprintf ppf "@[%4s%s@]@." "" line;
+  Format.fprintf ppf "@[%4s@[%5s@]@,@[%16s@]@,@[%16s@]@,@[%11s@]@]@." 
+                      ""    "#" "HF energy  " "Convergence" "HOMO-LUMO";
+  Format.fprintf ppf "@[%4s%s@]@." "" line;
+  let nocc = nocc t in
+  Array.iter (fun data ->
+    let data = Lazy.force data in
+    match data with
+    | None -> ()
+    | Some data ->
+      let e = Util.of_some data.eigenvalues in
+      let gap = e%.(nocc+1) -. e%.(nocc) in
+      begin
+        Format.fprintf ppf "@[%4s@[%5d@]@,@[%16.8f@]@,@[%16.4e@]@,@[%11.4f@]@]@." ""
+          (data.iteration) (Util.of_some data.energy) (Util.of_some data.error) gap;
+      end
+    ) t.data;
+  Format.fprintf ppf "@[%4s%s@]@." "" line
+
+
+let pp_summary ppf t =
+  let print text value = 
+    Format.fprintf ppf "@[@[%30s@]@,@[%16.10f@]@]@;" text value;
+  and line () = 
+    Format.fprintf ppf "@[  %s  @]@;" (String.make (linewidth-4) '-');
+  in
+  let ha_to_ev = Constants.ha_to_ev in
+  let e = eigenvalues t in
+
+  Format.fprintf ppf "@[%s@]@;" (String.make linewidth '=')
+  ; Format.fprintf ppf "@[<v>"
+  ; print "One-electron energy" (kin_energy t +. eN_energy t) 
+  ; print "Kinetic"   (kin_energy t) 
+  ; print "Potential" (eN_energy  t) 
+  ; line () 
+  ; print "Two-electron energy" (coulomb_energy t +. exchange_energy t) 
+  ; print "Coulomb"  (coulomb_energy  t) 
+  ; print "Exchange" (exchange_energy t) 
+  ; line ()
+  ; print "HF HOMO" (ha_to_ev *. (e%.(nocc t)))
+  ; print "HF LUMO" (ha_to_ev *. (e%.(nocc t + 1)))
+  ; print "HF LUMO-LUMO" (ha_to_ev *. (e%.(nocc t + 1) -. e%.(nocc t)))
+  ; line ()
+  ; print "Electronic   energy" (energy t -. nuclear_repulsion t) 
+  ; print "Nuclear   repulsion" (nuclear_repulsion t) 
+  ; print "Hartree-Fock energy" (energy t) 
+  ; Format.fprintf ppf "@]"
+  ; Format.fprintf ppf "@[%s@]@;" (String.make linewidth '=')
+
+
+
+
+let pp ppf t =
+  Format.fprintf ppf "@.@[%s@]@." (String.make 70 '=');
+  Format.fprintf ppf "@[%34s %-34s@]@." (match t.kind with
+  | UHF  -> "Unrestricted"
+  | RHF  -> "Restricted"
+  | ROHF -> "Restricted Open-shell") "Hartree-Fock";
+  Format.fprintf ppf "@[%s@]@.@." (String.make 70 '=');
+  Format.fprintf ppf "@[%a@]@." pp_iterations t;
+  Format.fprintf ppf "@[<v 4>@;%a@]@." pp_summary    t

mo/lib/hartree_fock.mli

@@ -0,0 +1,75 @@
+open Linear_algebra
+
+(** Data structure representing the output of a Hartree-Fock caculation *)
+
+type hartree_fock_data
+
+type hartree_fock_kind =
+| RHF       (** Restricted Hartree-Fock *)
+| ROHF      (** Restricted Open-shell Hartree-Fock *)
+| UHF       (** Unrestricted Hartree-Fock *)
+
+
+type t
+type mo = Mo_dim.t
+type ao = Ao.Ao_dim.t
+
+val kind : t -> hartree_fock_kind
+(** Kind of simulation  : RHF, ROHF or UHF. *)
+
+val simulation : t -> Simulation.t
+(** Simulation which was used for HF calculation *)
+
+val guess : t -> Guess.t
+(** Initial guess *)
+
+val eigenvectors : t -> (ao, mo) Matrix.t
+(** Final eigenvectors *)
+
+val eigenvalues : t -> mo Vector.t
+(** Final eigenvalues  *)
+
+val occupation : t -> mo Vector.t
+(** Diagonal of the density matrix *)
+
+val energy : t -> float
+(** Final energy *)
+
+val nuclear_repulsion : t -> float
+(** Nucleus-Nucleus potential energy *)
+
+val kin_energy : t -> float
+(** Kinetic energy *)
+
+val eN_energy : t -> float
+(** Electron-nucleus potential energy *)
+                                                                                              
+val coulomb_energy : t -> float 
+(** Electron-Electron potential energy *)
+
+val exchange_energy : t -> float 
+(** Exchange energy *)
+
+val nocc : t -> int
+(** Number of occupied MOs *)
+
+val empty: hartree_fock_data
+(** Empty data *)
+
+
+val make :
+  ?kind:hartree_fock_kind ->
+  ?guess:[ `Hcore | `Huckel | `Matrix of (ao,mo) Matrix.t ] ->
+  ?max_scf:int ->
+  ?level_shift:float -> ?threshold_SCF:float ->
+  Simulation.t -> t
+
+
+(** {1 Printers} *)
+
+val pp : Format.formatter -> t -> unit
+
+val pp_iterations : Format.formatter -> t -> unit
+
+val pp_summary : Format.formatter -> t -> unit
+

mo/lib/mo_dim.ml

@@ -0,0 +1 @@
+type t

mo/lib/mo_dim.mli

@@ -0,0 +1 @@
+type t

mo/test/dune

@@ -0,0 +1,12 @@
+(library
+ (name test_mo_basis)
+ (libraries
+  alcotest
+  qcaml.mo
+ )
+ (synopsis "Tests for the MO basis"))
+
+
+
+
+

mo/test/guess.ml

@@ -0,0 +1,50 @@
+open Alcotest
+open Particles
+open Operators
+open Mo.Guess
+open Linear_algebra
+    
+let wd = Common.Qcaml.root ^ Filename.dir_sep ^ "test"
+
+let test_case ao_basis =
+
+  let test_hcore () =
+    match make ~guess:`Hcore ao_basis with
+    | Hcore matrix -> 
+      let a = Matrix.to_array matrix in
+      let reference =
+        Matrix.add 
+          (Ao.Basis.eN_ints  ao_basis)
+          (Ao.Basis.kin_ints ao_basis)
+        |> Matrix.to_array
+      in
+      Array.iteri (fun i x ->
+        let message =
+          Printf.sprintf "Guess line %d" (i)
+        in
+        check (array (float 1.e-15)) message a.(i) x) reference
+    | _ -> assert false
+  
+  in
+  [
+    "HCore", `Quick, test_hcore;
+  ]
+
+let tests = 
+  List.concat [ 
+    let nuclei =
+      wd ^ Filename.dir_sep ^ "water.xyz"
+    |> Nuclei.of_xyz_file 
+    in
+    let basis_filename =
+      wd ^ Filename.dir_sep ^ "cc-pvdz"
+    in
+    let rs = 0.5 in
+    let operators = [ Operator.of_range_separation rs ] in
+    let ao_basis = 
+      Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
+        ~operators ~cartesian:false ~nuclei basis_filename
+    in
+    test_case ao_basis ;
+  ]
+    

test/dune

@@ -8,6 +8,7 @@
   test_gaussian_basis
   test_gaussian_integrals
   test_ao_basis
+  test_mo_basis
 ))
 
 (alias

test/run_tests.ml

@@ -11,6 +11,7 @@ let test_suites: unit Alcotest.test list = [
   "Gaussian_basis.General_basis", Test_gaussian_basis.General_basis.tests;
   "Ao_basis.Ao_basis_gaussian", Test_ao_basis.Ao_basis_gaussian.tests;
   "Ao_basis.Ao_basis", Test_ao_basis.Ao_basis.tests;
+  "Mo_basis.Guess", Test_mo_basis.Guess.tests;
 ]
 
 let () = Alcotest.run "QCaml" test_suites