From 99f5fe0aefc42a590c6b6743dec83e307548be6d Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 25 Mar 2019 16:50:36 +0100
Subject: [PATCH] F12 matrix built as a stream

---
 Basis/F12.ml | 10 +++++++
 CI/F12CI.ml  | 85 +++++++++++++++++++++-------------------------------
 2 files changed, 44 insertions(+), 51 deletions(-)

diff --git a/Basis/F12.ml b/Basis/F12.ml
index 26a2fb6..cc2cec8 100644
--- a/Basis/F12.ml
+++ b/Basis/F12.ml
@@ -147,11 +147,13 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup
   and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple
   in
 
+  (*
   let s =
     Overlap.of_basis basis
     |> Overlap.matrix
   in
   let lambda_inv = 1. /. expo_s in
+  *)
   Array.iteri (fun i_c powers_i ->
       let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
       let xi = to_powers powers_i in
@@ -166,7 +168,10 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup
                   let xl = to_powers powers_l in
                   let key = Zkey.of_powers_twelve  xi xj xk xl in
                   let value = Zmap.find cls key in
+                  (*
                   lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. value)
+                  *)
+                  value
                   |> set_chem data i_c j_c k_c l_c
                 ) (Cs.zkey_array (Csp.shell_b shell_q))
             ) (Cs.zkey_array (Csp.shell_a shell_q))
@@ -359,15 +364,20 @@ let of_basis_parallel basis =
       Fis.create ~size:0 `Dense
   in
 
+    (*
   let s =
     Overlap.of_basis basis
     |> Overlap.matrix
   in
   let lambda_inv = 1. /. expo_s in
+  *)
   Farm.run ~ordered:false ~f input_stream
   |> Stream.iter (fun l ->
       Array.iter (fun (i_c,j_c,k_c,l_c,value) ->
+      (*
         lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. value)
+        *)
+      value
         |> set_chem eri_array i_c j_c k_c l_c) l);
 
   if Parallel.master then
diff --git a/CI/F12CI.ml b/CI/F12CI.ml
index b542be0..c42377a 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -112,65 +112,47 @@ let dressing_vector gamma ~frozen_core aux_basis f12_amplitudes ci =
 
 *)
 
-  let out_dets =
-    DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis
-    |> DeterminantSpace.determinants_array
-    |> Array.to_list
-    |> List.filter (fun i -> not (is_internal ci.CI.det_space i))
-  in
-
-  let in_dets =
-    DeterminantSpace.determinants_array ci.CI.det_space 
-    |> Array.to_list
-  in
-
   Printf.printf "Building matrix\n%!";
 
+
   let m_H_aux, m_F_aux = 
-    let rec col_vecs_list accu_H accu_F = function
-    | [] ->
-      List.rev_map Vec.of_list accu_H, 
-      List.rev_map Vec.of_list accu_F 
-    | ki :: rest ->
-      let h, f = 
-        List.map (fun kj ->
-        match hf_ij aux_basis ki kj with
-        | [ a ; b ] -> a, b
-        | _ -> assert false ) in_dets
-        |> List.split
-      in
-      col_vecs_list (h::accu_H) (f::accu_F) rest
+    let out_dets_stream =
+      DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis
+      |> DeterminantSpace.determinant_stream
+    in
+
+    let in_dets =
+      DeterminantSpace.determinants_array ci.CI.det_space 
+      |> Array.to_list
+    in
+
+    let rec col_vecs_list accu_H accu_F = 
+      try
+        let ki = Stream.next out_dets_stream in
+        if is_internal ci.CI.det_space ki then
+          raise Exit
+        else
+          let h, f = 
+            List.map (fun kj ->
+            match hf_ij aux_basis ki kj with
+            | [ a ; b ] -> a, b
+            | _ -> assert false ) in_dets
+            |> List.split
+          in
+          col_vecs_list (h::accu_H) (f::accu_F) 
+      with
+     | Exit -> col_vecs_list accu_H accu_F
+     | Stream.Failure ->
+        List.rev_map Vec.of_list accu_H, 
+        List.rev_map Vec.of_list accu_F 
     in
     let h, f = 
-      col_vecs_list [] [] out_dets
+      col_vecs_list [] [] 
     in
     Mat.of_col_vecs_list h, 
     Mat.of_col_vecs_list f
   in
 
-(*
-  let m_H_aux = 
-    List.map (fun ki ->
-      List.map (fun kj ->
-        h_ij aux_basis ki kj
-      ) out_dets
-      |> Array.of_list
-    ) in_dets
-    |> Array.of_list
-    |> Mat.of_array
-  in
-
-  let m_F_aux = 
-    List.map (fun ki ->
-      List.map (fun kj ->
-        f_ij aux_basis ki kj
-      ) out_dets
-      |> Array.of_list
-    ) in_dets
-    |> Array.of_list
-    |> Mat.of_array
-  in
-  *)
 
   Printf.printf "Matrix product\n%!";
   let m_HF = 
@@ -189,7 +171,7 @@ let dressing_vector gamma ~frozen_core aux_basis f12_amplitudes ci =
 
 let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
 
-  let gamma = 0.5 in
+  let gamma = -0.5 /. F12.expo_s  in
 
   let mo_num = MOBasis.size mo_basis in
 
@@ -316,7 +298,8 @@ TODO SINGLE STATE HERE
           (Mat.to_col_vecs eigenvectors).(0) )
       in
       if Parallel.master then
-        Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift);
+        Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift 
+          +. Simulation.nuclear_repulsion simulation);
 
       if conv > threshold then
         iteration eigenvectors