Introduced difference of Fock matrices for SCF convergence

SCF/Fock.ml

@@ -16,7 +16,7 @@ module Ao = AOBasis
 
 
 
-let make ~density ao_basis =
+let make ~density ?(threshold=Constants.epsilon) ao_basis =
   let m_P  = density
   and m_T  = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix 
   and m_V  = Lazy.force ao_basis.Ao.eN_ints  |> NucInt.matrix
@@ -30,51 +30,6 @@ let make ~density ao_basis =
   and m_K = Array.make_matrix nBas nBas 0.
   in
 
-(*
-  for sigma = 1 to nBas do
-    for nu = 1 to nBas do
-      let m_Jnu = m_J.(nu-1) in
-      for lambda = 1 to sigma do
-        let p = 
-          if lambda < sigma then
-            2. *. m_P.{lambda,sigma}
-          else
-            m_P.{lambda,sigma}
-        in
-        if abs_float p > epsilon then
-          for mu = 1 to nu do
-            m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  p *. 
-                          ERI.get_phys m_G mu lambda nu sigma
-          done
-      done
-    done
-  done;
-  for nu = 1 to nBas do
-    for mu = 1 to nu-1 do
-      m_J.(mu-1).(nu-1) <- m_J.(nu-1).(mu-1);
-    done
-  done;
-
-  for nu = 1 to nBas do
-    let m_Knu = m_K.(nu-1) in
-    for sigma = 1 to nBas do
-      for lambda = 1 to nBas do
-        let p = 
-            0.5 *. m_P.{lambda,sigma}
-        in
-        if abs_float p > epsilon then
-          for mu = 1 to nu do
-            m_Knu.(mu-1) <- m_Knu.(mu-1) -. p *.
-                          ERI.get_phys m_G mu lambda sigma nu
-          done
-      done
-    done;
-    for mu = 1 to nu-1 do
-      m_K.(mu-1).(nu-1) <- m_Knu.(mu-1);
-    done
-  done;
-  *)
-
   for sigma = 1 to nBas do
     let m_Ksigma = m_K.(sigma-1) in
     for nu = 1 to nBas do
@@ -83,7 +38,7 @@ let make ~density ao_basis =
         let pJ = m_P.{lambda,sigma}
         and pK = 0.5 *. m_P.{lambda,nu}
         in
-        match (abs_float pJ > epsilon , abs_float pK > epsilon, nu < sigma) with
+        match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with
         | (false, false, _) -> ()
         | (true , true , true) ->
             begin
@@ -147,6 +102,24 @@ let make ~density ao_basis =
     core = m_Hc ; coulomb = m_J ; exchange = m_K }
 
 
+
+let op ~f f1 f2 =
+  let m_Hc = f1.core
+  and m_J = f f1.coulomb f2.coulomb
+  and m_K = f f1.exchange f2.exchange
+  in
+  {
+    fock     = Mat.add m_Hc (Mat.add m_J m_K);
+    core     = m_Hc;
+    coulomb  = m_J;
+    exchange = m_K;
+  }
+
+
+let add = op ~f:(fun a b -> Mat.add a b)
+let sub = op ~f:(fun a b -> Mat.sub a b)
+
+
 let pp_fock ppf a =
   Format.fprintf ppf "@[<2>";
   Format.fprintf ppf "@[ Fock matrix:@[<2>@[%a@]@.]@]" pp_matrix a.fock;

SCF/RHF.ml

@@ -1,4 +1,5 @@
 open Util
+open Constants
 open Lacaml.D
 
 module Si = Simulation
@@ -52,7 +53,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
   in
 
   (* SCF iterations *)
-  let rec loop nSCF iterations energy_prev m_C diis =
+  let rec loop nSCF iterations energy_prev m_C m_P_prev fock_prev threshold diis =
 
       (* Density matrix over nocc occupied MOs *)
       let m_P =
@@ -60,18 +61,20 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
       in
 
       (* Fock matrix in AO basis *)
-      let m_F, m_Hc, m_J, m_K =
+      let fock = 
-        let x = 
+        match fock_prev, threshold > 100. *. threshold_SCF with
-          Fock.make  ~density:m_P  ao_basis
+        | Some fock_prev, true -> 
+            let threshold = 1.e-8 in
+            Fock.make  ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis
+            |> Fock.add fock_prev
+        | _ -> Fock.make  ~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
       in
-(*
+
-debug_matrix "Fock" m_F;
-debug_matrix "Coulomb" m_J;
-debug_matrix "Exchange" m_K;
-debug_matrix "HCore" m_Hc;
-*)
       (* Add level shift in AO basis *)
       let m_F =
         let m_SC =
@@ -149,7 +152,7 @@ debug_matrix "HCore" m_Hc;
       in
 
       if not converged then
-          loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C diis
+          loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C m_P (Some fock) error diis
       else
           let iterations = 
             List.rev ( (energy, error, gap) :: iterations )
@@ -188,7 +191,7 @@ debug_matrix "HCore" m_Hc;
   in
 
   let diis = DIIS.make () in
-  loop 1 [] None m_C diis
+  loop 1 [] None m_C m_C None threshold_SCF diis