From 4dec632e5a750480c051dbacee72b39434731cfb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 4 Jul 2018 19:21:45 +0200 Subject: [PATCH] Introduced difference of Fock matrices for SCF convergence --- SCF/Fock.ml | 67 ++++++++++++++++------------------------------------- SCF/RHF.ml | 27 +++++++++++---------- 2 files changed, 35 insertions(+), 59 deletions(-) diff --git a/SCF/Fock.ml b/SCF/Fock.ml index 7da46d7..ced2eac 100644 --- a/SCF/Fock.ml +++ b/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; diff --git a/SCF/RHF.ml b/SCF/RHF.ml index 04894e9..0deeca4 100644 --- a/SCF/RHF.ml +++ b/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 fock = + 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.add fock_prev + | _ -> Fock.make ~density:m_P ao_basis + in + let m_F, m_Hc, m_J, m_K = - let x = - Fock.make ~density:m_P ao_basis - in + 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