From 075455fc9c1a3b3e770d20febae0c11664381e9b Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 25 Mar 2019 10:54:23 +0100
Subject: [PATCH] Fixed bug in F12 integrals

---
 Basis/F12.ml | 12 +++++++++---
 CI/F12CI.ml  | 24 ++++++++++++++++--------
 run_cas.ml   |  2 +-
 3 files changed, 26 insertions(+), 12 deletions(-)

diff --git a/Basis/F12.ml b/Basis/F12.ml
index 3b377ed..425268f 100644
--- a/Basis/F12.ml
+++ b/Basis/F12.ml
@@ -14,8 +14,14 @@ module Fis = FourIdxStorage
 
 include FourIdxStorage
 
+(*
+type gaussian_geminal =
+{
+  expo_s : float ;
+  coef_g : 
+  *)
 (** Exponent of the geminal *)
-let expo_s = 2.5
+let expo_s = 1.
 
 (** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) 
 let coef_g =
@@ -121,7 +127,7 @@ 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
+                  lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. 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))
@@ -321,7 +327,7 @@ let of_basis_parallel basis =
   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
+        lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. 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 81f3545..f920c9f 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -128,6 +128,7 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
   in
 
   Printf.printf "Building matrix\n%!";
+
   let m_H_aux, m_F_aux = 
     let rec col_vecs_list accu_H accu_F = function
     | [] ->
@@ -137,35 +138,39 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
       let h, f = 
         List.map (fun kj ->
         match hf_ij aux_basis ki kj with
-        | [ a ; b ] -> a, b
-        | _ -> assert false ) out_dets
+        | [ a ; b ] -> a, gamma *. b
+        | _ -> assert false ) in_dets
         |> List.split
       in
       col_vecs_list (h::accu_H) (f::accu_F) rest
     in
     let h, f = 
-      col_vecs_list [] [] in_dets
+      col_vecs_list [] [] out_dets
     in
     Mat.of_col_vecs_list h, 
     Mat.of_col_vecs_list f
   in
 
-  (*
+(*
   let m_H_aux = 
-    Array.map (fun ki ->
-      Array.map (fun kj ->
+    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 = 
-    Array.map (fun ki ->
-      Array.map (fun kj ->
+    List.map (fun ki ->
+      List.map (fun kj ->
         f_ij gamma aux_basis ki kj
       ) out_dets
+      |> Array.of_list
     ) in_dets
+    |> Array.of_list
     |> Mat.of_array
   in
   *)
@@ -175,6 +180,9 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
     gemm m_H_aux m_F_aux ~transb:`T 
   in
   Printf.printf "Done\n%!";
+  Printf.printf "%d %d %d %d %d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux)
+    (Mat.dim1 m_F_aux) (Mat.dim2 m_F_aux)
+    (Mat.dim1 m_HF) (Mat.dim2 m_HF);
   gemm m_HF f12_amplitudes
   |> Matrix.sparse_of_mat
 
diff --git a/run_cas.ml b/run_cas.ml
index 8e65248..89a91e2 100644
--- a/run_cas.ml
+++ b/run_cas.ml
@@ -69,7 +69,7 @@ let () =
   let space =
     DeterminantSpace.cas_of_mo_basis mos ~frozen_core:true n m
   in
-  let ci = CI.make ~frozen_core:true space in
+  let ci = CI.make space in
   Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
 
 (*