From 69faafb6f797a928228a81c6c45d181896b4a83e Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 26 Mar 2019 18:02:22 +0100
Subject: [PATCH] FCIF12 OK for 2 electrons

---
 Basis/F12.ml       | 40 +++++++++++++++++++++++-----------------
 CI/F12CI.ml        | 17 +++++++++++++++--
 Utils/Bitstring.ml |  2 +-
 3 files changed, 39 insertions(+), 20 deletions(-)

diff --git a/Basis/F12.ml b/Basis/F12.ml
index 22ffaa2..4f6e6e5 100644
--- a/Basis/F12.ml
+++ b/Basis/F12.ml
@@ -31,12 +31,24 @@ let make_gaussian_corr_factor expo_s coef_g expo_sg =
   }
 
 
-(* exp (-expo_s r) *)
+(* -1/expo_s *. exp (-expo_s r) *)
 let gaussian_geminal expo_s =
   let coef_g =
-    [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |];
+    [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |]
+    |> Array.map (fun x -> -. x /. expo_s)
   and expo_sg =
-    [| 0.2209 ; 1.004 ;  3.622 ; 12.16 ;   45.87 ;  254.4 |];
+    [| 0.2209 ; 1.004 ;  3.622 ; 12.16 ;   45.87 ;  254.4 |]
+  in
+  make_gaussian_corr_factor expo_s coef_g expo_sg
+
+
+
+(* exp (-expo_s r) *)
+let simple_gaussian_geminal expo_s =
+  let coef_g =
+    [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |]
+  and expo_sg =
+    [| 0.2209 ; 1.004 ;  3.622 ; 12.16 ;   45.87 ;  254.4 |]
   in
   make_gaussian_corr_factor expo_s coef_g expo_sg
 
@@ -52,7 +64,7 @@ let gaussian_geminal_times_r12 expo_s =
 
 
 (* exp (-expo_s r) *)
-let gaussian_geminal' expo_s =
+let simple_gaussian_geminal' expo_s =
   let coef_g =
     [|
     -3.4793465193721626604883567779324948787689208984375 ;
@@ -135,7 +147,7 @@ let filter_contracted_shell_pairs f12 ?(cutoff=integrals_cutoff) shell_pairs =
 *)
 
 
-let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = 
+let store_class basis ?(cutoff=integrals_cutoff) data f12 contracted_shell_pair_couple cls = 
   let to_powers x = 
     let open Zkey in
     match to_powers x with
@@ -148,11 +160,7 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup
   in
 
   (*
-  let s =
-    Overlap.of_basis basis
-    |> Overlap.matrix
-  in
-  let lambda_inv = 1. /. f12.expo_s in
+  let lambda_inv = -. 1. /. f12.expo_s in
   *)
   Array.iteri (fun i_c powers_i ->
       let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
@@ -170,6 +178,7 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup
                   let value = Zmap.find cls key in
                   (*
                   lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. value)
+                  lambda_inv *. value
                   *)
                   value
                   |> set_chem data i_c j_c k_c l_c
@@ -240,7 +249,7 @@ let of_basis_serial f12 basis =
               let cls =
                 class_of_contracted_shell_pair_couple f12 cspc
               in
-              store_class basis ~cutoff eri_array cspc cls
+              store_class basis ~cutoff eri_array f12 cspc cls
             | None -> ()
           ) shell_pairs
       with Exit -> ()
@@ -364,18 +373,15 @@ let of_basis_parallel f12 basis =
       Fis.create ~size:0 `Dense
   in
 
-    (*
-  let s =
-    Overlap.of_basis basis
-    |> Overlap.matrix
-  in
-  let lambda_inv = 1. /. f12.expo_s in
+  (*
+  let lambda_inv = -. 1. /. f12.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)
+        lambda_inv *. value
         *)
         value
         |> set_chem eri_array i_c j_c k_c l_c) l);
diff --git a/CI/F12CI.ml b/CI/F12CI.ml
index 772cb49..b4f6704 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -244,6 +244,11 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
 
 
   let f12_amplitudes = 
+    fun c -> 
+      let result = lacpy c in
+      Mat.scal (0.5) result ; 
+      result
+  (*
     (* While in a sequential region, initiate the parallel
        4-idx transformation to avoid nested parallel jobs
     *)
@@ -262,8 +267,16 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
       |> Lazy.force
     in
     fun ci_coef -> 
-    Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef)
-    |> Matrix.to_mat
+      let result = 
+        Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef)
+        |> Matrix.to_mat
+      in
+      let norm = gemm ~transa:`T result result in
+      Printf.printf "Norm of |F> : %f\n%!" norm.{1,1};
+      if abs_float norm.{1,1}  > 1. then
+        failwith "Norm of |F> > 1";
+      result
+      *)
   in
 
   let e_shift = 
diff --git a/Utils/Bitstring.ml b/Utils/Bitstring.ml
index 3e6382c..ff86504 100644
--- a/Utils/Bitstring.ml
+++ b/Utils/Bitstring.ml
@@ -189,7 +189,7 @@ let permtutations m n =
       let t'' = shift_right (minus_one (logand (lognot t) t')) (trailing_zeros u + 1) in
       aux (k-1) (logor t' t'') (u :: rest)
   in
-  aux (Util.binom n m) (minus_one (shift_left_one m m)) []
+  aux (Util.binom n m) (minus_one (shift_left_one n m)) []
 
 
 (*-----------------------------------------------------------------------------------*)