ROHF OK

CI/CI.ml

@@ -5,16 +5,16 @@ module Ds = Determinant_space
 type t =
   {
     det_space   : Ds.t ;
-    h_matrix    : Mat.t lazy_t ;
+    m_H         : Mat.t lazy_t ;
-    s2_matrix   : Mat.t lazy_t ;
+    m_S2        : Mat.t lazy_t ;
     eigensystem : (Mat.t * Vec.t) lazy_t;
   }
 
 let det_space    t = t.det_space
 
-let h_matrix     t = Lazy.force t.h_matrix
+let m_H          t = Lazy.force t.m_H
 
-let s2_matrix    t = Lazy.force t.s2_matrix
+let m_S2         t = Lazy.force t.m_S2
 
 let eigensystem  t = Lazy.force t.eigensystem 
 
@@ -54,7 +54,7 @@ let make det_space =
   let det  = Ds.determinants det_space in
   let mo_basis = Ds.mo_basis det_space in
   (*
-  let h_matrix = lazy (
+  let m_H = lazy (
     Util.list_range 0 (ndet-1)
     |> List.map (fun i -> let ki = det.(i) in
                   Array.init ndet (fun j -> let kj = det.(j) in
@@ -66,7 +66,7 @@ let make det_space =
   in
   *)
   (*
-  let h_matrix = lazy (
+  let m_H = lazy (
     let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in
     List.init ntasks (fun i -> 
       let m = 
@@ -93,7 +93,7 @@ let make det_space =
     |> Mat.of_col_vecs_list
   ) in
   *)
-  let h_matrix = lazy (
+  let m_H = lazy (
     let h =
       if Parallel.master then
         Array.make_matrix ndet ndet 0.
@@ -120,7 +120,7 @@ let make det_space =
       List.iter (fun (i,j,x) -> h.{i+1,j+1} <- x) l);
     Parallel.broadcast (lazy h)
   ) in
-  let s2_matrix = lazy (
+  let m_S2 = lazy (
     Array.init ndet (fun i -> let ki = det.(i) in
                       Array.init ndet (fun j -> let kj = det.(j) in
                                         CIMatrixElement.make_s2 ki kj
@@ -129,10 +129,10 @@ let make det_space =
   )
   in
   let eigensystem = lazy (
-    let h_matrix = Lazy.force h_matrix in
+    let m_H = Lazy.force m_H in
     Parallel.broadcast @@
-      lazy (Util.diagonalize_symm h_matrix)
+      lazy (Util.diagonalize_symm m_H)
   )
   in
-  { det_space ; h_matrix ; s2_matrix ; eigensystem }
+  { det_space ; m_H ; m_S2 ; eigensystem }
 

MOBasis/MOBasis.ml

@@ -45,8 +45,7 @@ let mo_matrix_of_ao_matrix ~mo_coef ao_matrix =
 
 let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix =
   let sc = gemm ao_overlap mo_coef  in
-  gemm sc @@
+  x_o_xt ~x:sc ~o:mo_matrix
-  gemm ~transb:`T mo_matrix sc
 
 
 let four_index_transform ~mo_coef eri_ao =

SCF/Fock.ml

@@ -12,11 +12,16 @@ type t =
     exchange : Mat.t ;
   }
 
+
+let fock t = t.fock
+let core t = t.core
+let coulomb t = t.coulomb
+let exchange t = t.exchange
+
+
 module Ao = AOBasis
 
-
+let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis =
-
-let make ~density ?(threshold=Constants.epsilon) ao_basis =
   let m_P  = density
   and m_T  = Ao.kin_ints ao_basis |> KinInt.matrix 
   and m_V  = Ao.eN_ints  ao_basis |> NucInt.matrix
@@ -41,47 +46,47 @@ let make ~density ?(threshold=Constants.epsilon) ao_basis =
         match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with
         | (false, false, _) -> ()
         | (true , true , true) ->
-            begin
+          begin
-              for mu = 1 to nu do
-                let integral = 
-                  ERI.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 *.
-                              ERI.get_phys m_G mu lambda nu sigma
-              done
-            end
-        | (true , true , false) ->
-           begin
-              for mu = 1 to sigma do
-                let integral = 
-                  ERI.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 *.
-                              ERI.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 *.
-                            ERI.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 *. 
+              let integral = 
-                            ERI.get_phys m_G mu lambda nu sigma
+                ERI.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 *.
+                                                    ERI.get_phys m_G mu lambda nu sigma
             done
+          end
+        | (true , true , false) ->
+          begin
+            for mu = 1 to sigma do
+              let integral = 
+                ERI.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 *. ERI.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 *. ERI.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 *. ERI.get_phys m_G mu lambda nu sigma
+          done
       done
     done;
     for mu = 1 to sigma-1 do
@@ -103,7 +108,100 @@ let make ~density ?(threshold=Constants.epsilon) ao_basis =
 
 
 
+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.kin_ints ao_basis |> KinInt.matrix 
+  and m_V  = Ao.eN_ints  ao_basis |> NucInt.matrix
+  and m_G  = Ao.ee_ints  ao_basis 
+  in
+  let nBas = Mat.dim1 m_T
+  in
+
+  let m_Hc = Mat.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 = 
+                ERI.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 *.
+                                                    ERI.get_phys m_G mu lambda nu sigma
+            done
+          end
+        | (true , true , false) ->
+          begin
+            for mu = 1 to sigma do
+              let integral = 
+                ERI.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 *. ERI.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 *. ERI.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 *. ERI.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 =  Mat.of_array m_J
+  and m_K =  Mat.of_array m_K
+  in
+  { fock = Mat.add m_Hc (Mat.add 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
@@ -117,8 +215,24 @@ let op ~f f1 f2 =
 
 
 let add = op ~f:(fun a b -> Mat.add a b)
+
 let sub = op ~f:(fun a b -> Mat.sub a b)
 
+let scale alpha f1 = 
+  let m_Hc = f1.core
+  and m_J = lacpy f1.coulomb 
+  and m_K = lacpy f1.exchange 
+  in
+  Mat.scal alpha m_J;
+  Mat.scal alpha m_K;
+  {
+    fock     = Mat.add m_Hc (Mat.add m_J m_K);
+    core     = m_Hc;
+    coulomb  = m_J;
+    exchange = m_K;
+  }
+
+
 
 let pp_fock ppf a =
   Format.fprintf ppf "@[<2>";

SCF/Fock.mli

@@ -0,0 +1,50 @@
+(** Type for the Fock operator in AO basis. *)
+
+open Lacaml.D
+
+
+type t 
+
+(** {1 Accessors} *)
+
+val fock : t -> Mat.t
+(** Fock matrix in AO basis *)
+
+val core : t -> Mat.t
+(** Core Hamiltonian : {% $\langle i | \hat{h} | j \rangle$ %} *)
+
+val coulomb : t -> Mat.t
+(** Coulomb matrix : {% $\langle i | J | j \rangle$ %} *)
+
+val exchange : t -> Mat.t
+(** Exchange matrix : {% $\langle i | K | j \rangle$ %} *)
+
+
+(** {1 Creators} *)
+
+val make_rhf : density:Mat.t -> ?threshold:float -> AOBasis.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: Mat.t -> density_other:Mat.t -> ?threshold:float ->
+               AOBasis.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. *)
+
+
+(** {1 Printers} *)
+
+val pp_fock : Format.formatter -> t -> unit
+

SCF/Guess.ml

@@ -34,8 +34,8 @@ let huckel_guess ao_basis =
   | 0 -> invalid_arg "Huckel guess needs a non-zero number of occupied MOs."
   | nocc ->
     let density = gemm ~alpha:2. ~transb:`T ~k:nocc m_X m_X in
-    let fock = Fock.make ~density ao_basis in
+    let f = Fock.make_rhf ~density ao_basis in
-    let m_F = fock.Fock.fock in
+    let m_F = Fock.fock f in
     for j=1 to ao_num do
       for i=1 to ao_num do
         if (i <> j) then

SCF/HartreeFock.ml

@@ -4,7 +4,7 @@ let make ?guess simulation =
   if Electrons.multiplicity @@ Simulation.electrons simulation  = 1 then
     RHF.make ?guess simulation
   else
-    invalid_arg "UHF or ROHF not implemented"
+    ROHF.make ?guess simulation
 
 
 let to_string = HartreeFock_type.to_string

SCF/RHF.ml

@@ -65,14 +65,14 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
         match fock_prev, threshold > 100. *. threshold_SCF with
         | Some fock_prev, true -> 
             let threshold = 1.e-8 in
-            Fock.make  ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis
+            Fock.make_rhf  ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis
             |> Fock.add fock_prev
-        | _ -> Fock.make  ~density:m_P ao_basis
+        | _ -> Fock.make_rhf  ~density:m_P ao_basis
       in
 
       let m_F, m_Hc, m_J, m_K =
         let x = fock in
-        x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange
+        Fock.(fock x, core x, coulomb x, exchange x)
       in
 
       (* Add level shift in AO basis *)

SCF/ROHF.ml

@@ -0,0 +1,290 @@
+open Util
+open Constants
+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 =
+
+  (* Number of occupied MOs *)
+  let n_alfa =
+    El.n_alfa @@ Si.electrons simulation
+  in
+
+  let nocc = n_alfa in
+
+  let n_beta = 
+    El.n_beta @@ Si.electrons simulation
+  in
+
+  let nuclear_repulsion = 
+    Si.nuclear_repulsion simulation
+  in
+
+  let ao_basis = 
+    Si.ao_basis simulation
+  in
+
+  (* Initial guess *)
+  let guess = 
+    Guess.make ~nocc ~guess ao_basis
+  in
+
+  (* Orthogonalization matrix *)
+  let m_X = 
+    Ao.ortho ao_basis
+  in
+
+
+  (* Overlap matrix *)
+  let m_S =
+    Ao.overlap ao_basis
+    |> Ov.matrix
+  in
+
+  let m_T = Ao.kin_ints ao_basis |> KinInt.matrix
+  and m_V = Ao.eN_ints  ao_basis |> NucInt.matrix
+  in
+
+  (* Level shift in MO basis *)
+  let m_LSmo =
+    Array.init (Mat.dim2 m_X) (fun i ->
+        if i > n_alfa then level_shift else 0.)
+    |> Vec.of_array
+    |> Mat.of_diag
+  in
+
+  (* SCF iterations *)
+  let rec loop nSCF iterations energy_prev m_C m_P_a_prev m_P_b_prev fock_a_prev fock_b_prev threshold diis =
+
+
+    (* Density matrix *)
+    let m_P_a =
+      gemm  ~alpha:1.  ~transb:`T  ~k:n_alfa  m_C  m_C
+    in
+
+    let m_P_b =
+      gemm  ~alpha:1.  ~transb:`T  ~k:n_beta  m_C  m_C
+    in
+
+    let m_P =
+      Mat.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:(Mat.sub m_P_a m_P_a_prev) ~density_other:(Mat.sub m_P_b 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:(Mat.sub m_P_b m_P_b_prev) ~density_other:(Mat.sub m_P_a 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, m_J_a, m_K_a =
+      let x = fock_a in
+      Fock.(fock x, core x, coulomb x, exchange x)
+    in
+
+    let m_F_b, m_Hc_b, m_J_b, m_K_b =
+      let x = fock_b in
+      Fock.(fock x, core x, coulomb x, exchange x)
+    in
+
+
+    let m_F_mo_a = 
+      xt_o_x ~o:m_F_a ~x:m_C
+    in
+
+    let m_F_mo_b = 
+      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 (Mat.dim2 m_F_mo_a) (fun i ->
+          let i = i+1 in
+          Array.init (Mat.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   ->
+                (*
+                0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) +.
+                0.5 *. (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
+                *)
+                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  ->
+                (*
+                0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) -.
+                0.5 *. (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
+                *)
+                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})
+            ) )
+      |> Mat.of_array
+    in
+
+    let m_SC =
+      gemm  m_S  m_C
+    in
+
+    let m_F = 
+      x_o_xt ~x:m_SC ~o:m_F_mo
+    in
+
+
+    (* Add level shift in AO basis *)
+    let m_F =
+      x_o_xt ~x:m_SC ~o:m_LSmo
+      |> Mat.add m_F
+    in
+
+    (* Fock matrix in orthogonal basis *)
+    let m_F_ortho =
+      xt_o_x m_F m_X
+    in
+
+    let error_fock = 
+      let fps =
+        gemm  m_F  (gemm  m_P  m_S)
+      and spf =
+        gemm  m_S  (gemm  m_P  m_F)
+      in
+      xt_o_x  (Mat.sub  fps  spf) m_X
+    in
+
+    let diis = 
+      DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
+    in
+
+    let m_F_diis =
+      let x = 
+        Bigarray.genarray_of_array1 (DIIS.next diis)
+      in
+      Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
+    in
+
+
+    (* MOs in orthogonal MO basis *)
+    let m_C', eigenvalues =
+      diagonalize_symm  m_F_diis
+    in
+
+    (* MOs in AO basis *)
+    let m_C =
+      gemm  m_X  m_C'
+    in
+
+    (* Hartree-Fock energy *)
+    let energy =
+      nuclear_repulsion +.  0.5 *. ( Mat.gemm_trace  m_P_a  (Mat.add  m_Hc_a  m_F_a) +.
+                                     Mat.gemm_trace  m_P_b  (Mat.add  m_Hc_b  m_F_b) )
+    in
+
+    (* Convergence criterion *)
+    let error =
+      error_fock 
+      |> Mat.as_vec 
+      |> amax
+      |> abs_float
+    in
+
+    let converged =
+      nSCF = max_scf  || error < threshold_SCF
+    in
+
+    let gap =
+      if nocc < Vec.dim eigenvalues then
+        eigenvalues.{nocc+1} -. eigenvalues.{nocc} 
+      else 0.
+    in
+
+    let () = 
+      match energy_prev with
+      | Some energy_prev ->
+        Printf.eprintf "%3d  %16.10f  %16.10f %11.4e  %10.4f\n%!" nSCF energy (energy -. energy_prev) error gap
+      | None -> 
+        Printf.eprintf "%3d  %16.10f  %16s %11.4e  %10.4f\n%!" nSCF energy "" error gap
+    in
+
+    if not converged then
+      loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C m_P_a m_P_b (Some fock_a) (Some fock_b) error diis
+    else
+      let iterations = 
+        List.rev ( (energy, error, gap) :: iterations )
+        |> Array.of_list
+      in
+      HartreeFock_type.(ROHF
+                          { 
+                            simulation;
+                            nocc;
+                            guess ;
+                            eigenvectors = m_C ;
+                            eigenvalues ;
+                            energy ;
+                            nuclear_repulsion;
+                            iterations ;
+                            kin_energy = Mat.gemm_trace m_P m_T;
+                            eN_energy  = Mat.gemm_trace m_P m_V;
+                            coulomb_energy  = 0.5 *. (Mat.gemm_trace m_P_a m_J_a) +. 
+                                              0.5 *. (Mat.gemm_trace m_P_b m_J_b);
+                            exchange_energy = 0.5 *. (Mat.gemm_trace m_P_a m_K_a) +.
+                                              0.5 *. (Mat.gemm_trace m_P_b m_K_b);
+                          })
+  in
+
+
+  (* Guess coefficients *)
+  let m_H = 
+    match guess with
+    | Guess.Hcore m_H -> m_H 
+    | Guess.Huckel m_H -> m_H 
+  in
+  let m_Hmo =
+    xt_o_x  m_H  m_X
+  in
+  let m_C', _ =
+    diagonalize_symm  m_Hmo
+  in
+  let m_C =
+    gemm  m_X  m_C'
+  in
+
+  let diis = DIIS.make () in
+  loop 1 [] None m_C m_C m_C None None threshold_SCF diis 
+

Utils/Util.ml

@@ -227,6 +227,10 @@ let xt_o_x ~o ~x =
   gemm o x
   |> gemm ~transa:`T x
 
+let x_o_xt ~o ~x =
+  gemm o x ~transb:`T 
+  |> gemm x
+
 
 let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c =
   let d, u, _ = gesvd (lacpy overlap) in

Utils/Util.mli

@@ -84,6 +84,9 @@ 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 {% $\mathbf{X^\dag\, O\, X}$ %} *)
 
+val x_o_xt : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
+(** Computes {% $\mathbf{X\, O\, X^\dag}$ %} *)
+
 val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat
 (** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %},
     and the last argument contains the vectors {% $\mathbf{C}$ %} to orthogonalize.

run_fci.ml

@@ -0,0 +1,68 @@
+let () =
+  let open Command_line in
+  begin
+    set_header_doc (Sys.argv.(0) ^ " - QCaml command");
+    set_description_doc "Runs a Hartree-Fock calculation";
+    set_specs
+      [ { short='b' ; long="basis" ; opt=Mandatory;
+          arg=With_arg "<string>";
+          doc="Name of the file containing the basis set"; } ;
+
+        { short='x' ; long="xyz" ; opt=Mandatory;
+          arg=With_arg "<string>";
+          doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
+
+        { short='m' ; long="multiplicity" ; opt=Optional;
+          arg=With_arg "<int>";
+          doc="Spin multiplicity (2S+1). Default is singlet"; } ;
+
+        { short='c' ; long="charge" ; opt=Optional;
+          arg=With_arg "<int>";
+          doc="Total charge of the molecule. Default is 0"; } ;
+      ]
+  end;
+
+  (* Handle options *)
+  let basis_file  = 
+    match Command_line.get "basis" with
+    | Some x -> x
+    | None -> assert false
+  in
+
+  let nuclei_file = 
+    match Command_line.get "xyz" with
+    | Some x -> x
+    | None -> assert false
+  in
+
+  let charge = 
+    match Command_line.get "charge" with
+    | Some x -> int_of_string x
+    | None   -> 0
+  in
+
+  let multiplicity =
+    match Command_line.get "multiplicity" with
+    | Some x -> int_of_string x
+    | None -> 1
+  in
+
+  let s =
+    Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file
+  in
+
+  let hf = HartreeFock.make s in
+  let mos =
+    MOBasis.of_hartree_fock hf
+  in
+  let space =
+    Determinant_space.fci_of_mo_basis ~frozen_core:false mos
+  in
+  let ci = CI.make space in
+  Format.printf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
+  (*
+  let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in
+  Util.list_range 1 (Determinant_space.size space)
+  |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i});
+  *)
+