From 29ddf41c0bfcfa6b32050675f5c4f4e4d05a606f Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 22 Mar 2019 19:15:59 +0100
Subject: [PATCH] F12 almost working

---
 CI/CI.ml    |  2 +-
 CI/F12CI.ml | 46 ++++++++++++++++++++++++++++++----------------
 run_fci.ml  |  2 +-
 3 files changed, 32 insertions(+), 18 deletions(-)

diff --git a/CI/CI.ml b/CI/CI.ml
index 1476f47..ea31578 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -737,7 +737,7 @@ let pt2_en_reference ci =
 
   let aux_basis = mo_basis in
   let ds =
-    DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis
+    DeterminantSpace.fci_of_mo_basis ~frozen_core:true aux_basis
   in
   let out_dets =
     ds
diff --git a/CI/F12CI.ml b/CI/F12CI.ml
index 52f1117..b338be3 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -136,7 +136,6 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
     |> Mat.of_array
   in
 
-Printf.printf "%d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux);
   let m_HF = 
     gemm m_H_aux m_F_aux ~transb:`T 
   in
@@ -149,7 +148,7 @@ Printf.printf "%d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux);
 
 let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
 
-  let gamma = 1. in
+  let gamma = 1.0 in
 
   let mo_num = MOBasis.size mo_basis in
 
@@ -236,9 +235,14 @@ Util.debug_matrix "T" (f12_amplitudes psi);
       let delta = 
          dressing_vector aux_basis (f12_amplitudes psi) ci
       in
-Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
+      (*
+Util.debug_matrix "psi" psi;
+*)
+Format.printf "Amplitude (1,1) : %f@." (f12_amplitudes psi).{1,1};
+Format.printf "Dressing vector(1,1) : %f@." (Matrix.get delta 1 1);
 
 (*------
+TODO SINGLE STATE HERE
 *)
       let m_H_dressed = Matrix.to_mat m_H in
       let f = 1.0 /. psi.{1,1} in
@@ -251,6 +255,9 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
             m_H_dressed.{1,i} <- m_H_dressed.{1,i} +. delta;
             m_H_dressed.{1,1} <- m_H_dressed.{1,1} -. delta *. psi.{i,1} *. f;
           end
+      (* DIAGONAL DRESSING
+      m_H_dressed.{i,i} <- m_H_dressed.{i,i} +. (Matrix.get delta i 1) /. (psi.{i,1} +. 1.e-5);
+          *)
       );
       let eigenvectors, eigenvalues =
         Util.diagonalize_symm m_H_dressed
@@ -260,7 +267,7 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
           (Mat.to_col_vecs psi).(0)
           (Mat.to_col_vecs eigenvectors).(0) )
       in
-      Printf.printf "Convergence : %f  %f\n" conv (eigenvalues.{1} +. e_shift);
+      Printf.printf "Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift);
       if conv > threshold then
         iteration eigenvectors
       else
@@ -272,13 +279,20 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
     iteration ci_coef
 (*
 -------    *)
-
 (*
-(*TODO : enlever le double comptage de la symmetrisation*)
       let n_states = ci.CI.n_states in
-      let diagonal =
-        Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i +. (if i=1 then Matrix.get delta 1 1 else 0.) )
+      let f = 1.0 /. (psi.{1,1} ) in
+      let delta_00 =
+        Util.list_range 2 (Mat.dim1 psi)
+        |> List.fold_left (fun accu i -> accu +.
+            (Matrix.get delta i 1) *. psi.{i,1} *. f) 0. 
       in
+      let delta = Matrix.to_mat delta in
+      delta.{1,1} <- delta.{1,1} -. delta_00;
+      let diagonal =
+        Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i +. (if i=1 then delta.{1,1} *. psi.{1,1} else 0.) )
+      in
+      let delta = Matrix.dense_of_mat delta in
       let matrix_prod c = 
         Matrix.add
           (Matrix.mm ~transa:`T m_H c)
@@ -286,24 +300,24 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
       in
       let eigenvectors, eigenvalues = 
         Parallel.broadcast (lazy (
-            Davidson.make ~threshold:1.e-6 ~guess:(Matrix.to_mat psi) ~n_states diagonal matrix_prod
+            Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
           ))
       in
-      let m =
-        Matrix.mm ~transa:`T psi (Matrix.dense_of_mat eigenvectors)
-        |> Matrix.to_mat
+      let conv =
+        1.0 -. abs_float ( dot
+          (Mat.to_col_vecs psi).(0)
+          (Mat.to_col_vecs eigenvectors).(0) )
       in
-      let conv = Mat.sum m -. (Vec.sum (Mat.copy_diag m)) in
-      Printf.printf "Convergence : %f  %f\n" conv eigenvalues.{1};
+      Printf.printf "Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift);
       if conv > threshold then
-        iteration (Matrix.dense_of_mat eigenvectors)
+        iteration eigenvectors
       else
         let eigenvalues =
           Vec.map (fun x -> x +. e_shift) eigenvalues
         in
         eigenvectors, eigenvalues
     in
-    iteration (Matrix.dense_of_mat ci_coef)
+    iteration ci_coef
 *)
 
   )
diff --git a/run_fci.ml b/run_fci.ml
index b95b7e4..1d1c886 100644
--- a/run_fci.ml
+++ b/run_fci.ml
@@ -56,7 +56,7 @@ let () =
   in
 
   let space =
-    DeterminantSpace.fci_of_mo_basis ~frozen_core:false mos
+    DeterminantSpace.fci_of_mo_basis  ~frozen_core:false  mos
   in
   let ci = CI.make space in
   Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);