Cleaned RHF

SCF/RHF.ml

@@ -7,6 +7,31 @@ module El = Electrons
 module Ao = AOBasis
 module Ov = Overlap
 
+
+type hartree_fock_data =
+  {
+    iteration       :  int              ;
+    coefficients    :  Mat.t   option   ;
+    eigenvalues     :  Vec.t   option   ;
+    error           :  float   option   ;
+    diis            :  DIIS.t  option   ;
+    energy          :  float   option   ;
+    density         :  Mat.t   option   ;
+    fock            :  Fock.t  option   ;
+  }
+
+let empty =
+  {
+    iteration       =  0                ;
+    coefficients    =  None             ;
+    eigenvalues     =  None             ;
+    error           =  None             ;
+    diis            =  None             ;
+    energy          =  None             ;
+    density         =  None             ;
+    fock            =  None             ;
+  }
+
 let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation =
   (* Number of occupied MOs *)
   let nocc =
@@ -50,8 +75,23 @@ let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation =
     |> Mat.of_diag
   in
 
-  (* SCF iterations *)
+
-  let rec loop nSCF iterations energy_prev m_C m_P_prev fock_prev threshold diis =
+  (* A single SCF iteration *)
+  let scf_iteration data =
+
+    let  nSCF         =  data.iteration + 1
+    and  m_C          =  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 =
@@ -60,8 +100,8 @@ let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation =
 
     (* Fock matrix in AO basis *)
     let fock = 
-        match fock_prev, threshold > 100. *. threshold_SCF with
+      match fock_prev, m_P_prev, threshold > 100. *. threshold_SCF with
-        | Some fock_prev, true -> 
+      | Some fock_prev, Some m_P_prev, true -> 
         let threshold = 1.e-8 in
         Fock.make_rhf  ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis
         |> Fock.add fock_prev
@@ -132,51 +172,60 @@ let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation =
       |> amax
       |> abs_float
     in
+    {
+      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         ;
+    }
 
-      let converged =
-        nSCF = max_scf  || error < threshold_SCF
   in
 
-      let gap =
+
-        if nocc < Vec.dim eigenvalues then
+  let rec make_iterations_list data =
-          eigenvalues.{nocc+1} -. eigenvalues.{nocc} 
+
-        else 0.
+    let energy_prev = data.energy in
+
+    (** Perform SCF iteration *)
+    let data = scf_iteration data in
+
+    (** Check convergence *)
+    let converged, error =
+      match data.error with
+      | None -> false, 0.
+      | Some error -> (data.iteration = max_scf || error < threshold_SCF), error
     in
 
+    (** Print values *)
+    let nSCF = data.iteration in
+
+    let energy = of_some data.energy 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
+        Printf.eprintf "%3d  %16.10f  %16.10f %11.4e\n%!" nSCF energy (energy -. energy_prev) error
       | None -> 
-          Printf.eprintf "%3d  %16.10f  %16s %11.4e  %10.4f\n%!" nSCF energy "" error gap
+        Printf.eprintf "%3d  %16.10f  %16s %11.4e\n%!" nSCF energy "" error
     in
 
-      if not converged then
+    if converged then
-          loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C m_P (Some fock) error diis
+      [ data ]
     else
-          let iterations = 
+      { empty  with 
-            List.rev ( (energy, error, gap) :: iterations )
+        iteration   = data.iteration;
-            |> Array.of_list
+        energy      = data.energy ;
-          in
+        eigenvalues = data.eigenvalues ;
-          HartreeFock_type.(RHF
+        error       = data.error ;
-          { 
+      } :: (make_iterations_list data)
-            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 m_J;
-            exchange_energy = 0.5 *. Mat.gemm_trace m_P m_K;
-            occupation = Mat.copy_diag m_P;
-            })
   in
 
 
+
   (* Guess coefficients *)
   let m_H = 
     match guess with
@@ -193,8 +242,50 @@ let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation =
     gemm  m_X  m_C'
   in
 
-  let diis = DIIS.make () in
+  let iterations_list =
-  loop 1 [] None m_C m_C None threshold_SCF diis 
+    make_iterations_list { empty with coefficients = Some m_C }
+  in
+
+  let iterations, data = 
+    List.map (fun data -> 
+        let gap =
+          let eigenvalues = of_some data.eigenvalues in
+          if nocc < Vec.dim eigenvalues then
+            eigenvalues.{nocc+1} -. eigenvalues.{nocc} 
+          else 0.
+        and energy = of_some data.energy
+        and error  = of_some data.error
+        in
+        (energy, error, gap)
+      ) iterations_list
+    |> Array.of_list,
+    List.hd (List.rev iterations_list)
+  in
+
+
+  let energy = of_some data.energy  in
+  let m_P    = of_some data.density in
+  let fock   = of_some data.fock    in
+  let m_J    = Fock.coulomb fock    in
+  let m_K    = Fock.exchange fock   in
+
+  HartreeFock_type.(
+    RHF { 
+      simulation;
+      nocc;
+      guess ;
+      eigenvectors = of_some data.coefficients ;
+      eigenvalues  = of_some data.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 m_J;
+      exchange_energy = 0.5 *. Mat.gemm_trace m_P m_K;
+      occupation = Mat.copy_diag m_P;
+    }
+  )
 
 
 

Utils/Util.ml

@@ -171,6 +171,11 @@ let boys_function ~maxm t =
     end
 
 
+let of_some = function
+  | Some a -> a
+  | None   -> assert false
+
+
 (** {2 List functions} *)
 
 let list_some l =

Utils/Util.mli

@@ -36,6 +36,7 @@ val chop : float -> (unit -> float) -> float
   than {!Constants.epsilon}, and return [a *. f ()].
   *)
 
+val of_some : 'a option -> 'a
 
 
 (** {2 Functions related to the Boys function} *)