Accelerated RHF

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_J = Array.make_matrix nBas nBas 0.
-  and m_K = Mat.make0 nBas nBas 
+  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
+    for mu = 1 to nu-1 do
-      m_J.{nu,mu} <- m_J.{mu,nu};
+      m_J.(mu-1).(nu-1) <- m_J.(nu-1).(mu-1);
-      m_K.{nu,mu} <- m_K.{mu,nu};
+      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 }
 

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
+      if m > 1 && sycon (lacpy a) > 1.e-14 then 
-      else  aux (m-1)
+        aux (m-1)
+      else a
     in
     aux diis.m
   in

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;
+                Bigarray.Genarray.set t [| i; j; k; l|] value;
-                t.{k,j,i,l} <- value;
+                Bigarray.Genarray.set t [| k; j; i; l|] value;
-                t.{i,l,k,j} <- value;
+                Bigarray.Genarray.set t [| i; l; k; j|] value;
-                t.{k,l,i,j} <- value;
+                Bigarray.Genarray.set t [| k; l; i; j|] value;
-                t.{j,i,l,k} <- value;
+                Bigarray.Genarray.set t [| j; i; l; k|] value;
-                t.{j,k,l,i} <- value;
+                Bigarray.Genarray.set t [| j; k; l; i|] value;
-                t.{l,i,j,k} <- value;
+                Bigarray.Genarray.set t [| l; i; j; k|] value;
-                t.{l,k,j,i} <- 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