Accelerated SCF

SCF/Fock.ml

@@ -29,31 +29,8 @@ let make ~density ao_basis =
   and m_J = Array.make_matrix nBas nBas 0.
   and m_K = Array.make_matrix nBas nBas 0.
   in
+
 (*
-  let permutations  i j k l  =
-      [  [ i ; j ; k ; l ] ; 
-         [ k ; j ; i ; l ] ; 
-         [ i ; l ; k ; j ] ; 
-         [ k ; l ; i ; j ] ; 
-         [ j ; i ; l ; k ] ; 
-         [ j ; k ; l ; i ] ; 
-         [ l ; i ; j ; k ] ; 
-         [ l ; k ; j ; i ]
-         ] 
-  in
-  ERI.to_stream m_G
-  |> Stream.iter (fun { ERI.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } ->
-     permutations i_r1 j_r2 k_r1 l_r2 
-     |> List.iter ( fun ijkl ->
-        match ijkl with
-        | mu :: lambda :: nu :: sigma :: [] ->
-          let p = m_P.{lambda,sigma} in
-          if abs_float p > epsilon then
-            let m_Jnu = m_J.(nu-1) in
-            m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. p *. value 
-        | _ -> assert false
-    ));
-  *)
   for sigma = 1 to nBas do
     for nu = 1 to nBas do
       let m_Jnu = m_J.(nu-1) in
@@ -96,6 +73,73 @@ let make ~density ao_basis =
       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
+      let m_Jnu = m_J.(nu-1) in
+      for lambda = 1 to nBas do
+        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
+        | (false, false, _) -> ()
+        | (true , true , true) ->
+            begin
+              for mu = 1 to nu do
+                let integral = 
+                  ERI.get_phys m_G mu lambda nu sigma
+                in
+                if (integral <> 0.) then begin
+                  m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  pJ *. integral;
+                  m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral
+                end
+              done;
+              for mu = nu+1 to sigma do
+                m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *.
+                              ERI.get_phys m_G mu lambda nu sigma
+              done
+            end
+        | (true , true , false) ->
+           begin
+              for mu = 1 to sigma do
+                let integral = 
+                  ERI.get_phys m_G mu lambda nu sigma
+                in
+                if (integral <> 0.) then begin
+                  m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  pJ *. integral;
+                  m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral
+                end
+              done;
+              for mu = sigma+1 to nu do
+                m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *.
+                              ERI.get_phys m_G mu lambda nu sigma
+              done
+            end
+        | (false, true , _) ->
+            for mu = 1 to sigma do
+              m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *.
+                            ERI.get_phys m_G mu lambda nu sigma
+            done
+        | (true , false, _) ->
+            for mu = 1 to nu do
+              m_Jnu.(mu-1) <- m_Jnu.(mu-1) +.  pJ *. 
+                            ERI.get_phys m_G mu lambda nu sigma
+            done
+      done
+    done;
+    for mu = 1 to sigma-1 do
+      m_K.(mu-1).(sigma-1) <- m_Ksigma.(mu-1);
+    done
+  done;
+  for nu = 1 to nBas do
+    let m_Jnu = m_J.(nu-1) in
+    for mu = 1 to nu-1 do
+      m_J.(mu-1).(nu-1) <- m_Jnu.(mu-1)
+    done
+  done;
+
   let m_J =  Mat.of_array m_J
   and m_K =  Mat.of_array m_K
   in

Utils/FourIdxStorage.ml

@@ -173,11 +173,6 @@ let to_stream d =
   and k = ref 1
   and l = ref 1
   in
-  let f i k = 
-    let p, r = 
-      if i <= k then i, k else k, i
-    in p+ (r*r-r)/2
-  in
   let rec f_dense _ =
     i := !i+1;
     if !i > !k then begin