From 15390759f984a9d132b652494b539f9ff11744fa Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 May 2018 18:50:34 +0200 Subject: [PATCH] Accelerated RHF --- SCF/Fock.ml | 19 ++++++++++++------- Utils/DIIS.ml | 5 +++-- Utils/FourIdxStorage.ml | 19 ++++++++++--------- 3 files changed, 25 insertions(+), 18 deletions(-) diff --git a/SCF/Fock.ml b/SCF/Fock.ml index 62749b3..20a5d8e 100644 --- a/SCF/Fock.ml +++ b/SCF/Fock.ml @@ -22,39 +22,44 @@ let make ~density simulation = in let m_Hc = Mat.add m_T m_V - and m_J = Mat.make0 nBas nBas - and m_K = Mat.make0 nBas nBas + 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 for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in for lambda = 1 to nBas do let p = m_P.{lambda,sigma} in if abs_float p > epsilon then for mu = 1 to nu do - m_J.{mu,nu} <- m_J.{mu,nu} +. p *. + 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 + 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_K.{mu,nu} <- m_K.{mu,nu} -. p *. + m_Knu.(mu-1) <- m_Knu.(mu-1) -. p *. ERI.get_phys m_G mu lambda sigma nu done done done done; for nu = 1 to nBas do - for mu = 1 to nu do - m_J.{nu,mu} <- m_J.{mu,nu}; - m_K.{nu,mu} <- m_K.{mu,nu}; + for mu = 1 to nu-1 do + m_J.(mu-1).(nu-1) <- m_J.(nu-1).(mu-1); + m_K.(mu-1).(nu-1) <- m_K.(nu-1).(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 } diff --git a/Utils/DIIS.ml b/Utils/DIIS.ml index 147a956..e3b5acc 100644 --- a/Utils/DIIS.ml +++ b/Utils/DIIS.ml @@ -44,8 +44,9 @@ let next diis = let a = Mat.make (m+1) (m+1) 1. in a.{m+1,m+1} <- 0.; ignore @@ lacpy ~b:a (gemm ~transa:`T ~m ~n:m e e); - if sycon (lacpy a) > 1.e-14 then a - else aux (m-1) + if m > 1 && sycon (lacpy a) > 1.e-14 then + aux (m-1) + else a in aux diis.m in diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index a90e1f9..97e7394 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -26,7 +26,8 @@ let key_of_indices ~r1 ~r2 = let get_four_index ~r1 ~r2 t = match t.four_index with - | Dense t -> let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in t.{i,j,k,l} + | Dense t -> let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in + Bigarray.Genarray.get t [| i; j; k; l|] | Sparse t -> let key = key_of_indices ~r1 ~r2 in try Hashtbl.find t key with Not_found -> 0. @@ -35,14 +36,14 @@ let get_four_index ~r1 ~r2 t = let set_four_index ~r1 ~r2 ~value t = match t.four_index with | Dense t -> let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in - t.{i,j,k,l} <- value; - t.{k,j,i,l} <- value; - t.{i,l,k,j} <- value; - t.{k,l,i,j} <- value; - t.{j,i,l,k} <- value; - t.{j,k,l,i} <- value; - t.{l,i,j,k} <- value; - t.{l,k,j,i} <- value; + Bigarray.Genarray.set t [| i; j; k; l|] value; + Bigarray.Genarray.set t [| k; j; i; l|] value; + Bigarray.Genarray.set t [| i; l; k; j|] value; + Bigarray.Genarray.set t [| k; l; i; j|] value; + Bigarray.Genarray.set t [| j; i; l; k|] value; + Bigarray.Genarray.set t [| j; k; l; i|] value; + Bigarray.Genarray.set t [| l; i; j; k|] value; + Bigarray.Genarray.set t [| l; k; j; i|] value; | Sparse t -> let key = key_of_indices ~r1 ~r2 in Hashtbl.replace t key value