From b40a9b8db8ea3b5d1429f5435b5bd280a8b83d80 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 22 Mar 2019 22:42:20 +0100
Subject: [PATCH] Working n F12

---
 CI/F12CI.ml | 42 +++++++++++++++++++++++-------------------
 1 file changed, 23 insertions(+), 19 deletions(-)

diff --git a/CI/F12CI.ml b/CI/F12CI.ml
index b338be3..2f19cd3 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -241,23 +241,29 @@ 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);
 
+      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;
+
 (*------
 TODO SINGLE STATE HERE
 *)
+(*
       let m_H_dressed = Matrix.to_mat m_H in
       let f = 1.0 /. psi.{1,1} in
       Util.list_range 1 (Mat.dim1 psi)
       |> List.iter (fun i ->
-        let delta = (Matrix.get delta i 1) *. f in
+        let delta = delta.{i,1} *. f in
         m_H_dressed.{i,1} <- m_H_dressed.{i,1} +. delta;
         if i <> 1 then
           begin
             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
@@ -277,26 +283,25 @@ TODO SINGLE STATE HERE
         eigenvectors, eigenvalues
     in
     iteration ci_coef
+    *)
 (*
 -------    *)
-(*
       let n_states = ci.CI.n_states in
-      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)
-          delta
+        let w = 
+          Matrix.mm ~transa:`T m_H c
+          |> Matrix.to_mat
+        in
+        Util.list_range 1 (Mat.dim1 w)
+        |> List.iter (fun i ->
+          w.{i,1} <- w.{i,1} +. delta.{i,1} ;
+          if (i <> 1) then
+            w.{1,1} <- w.{1,1} +. delta.{i,1} *. f *. (Matrix.get c i 1);
+          );
+        Matrix.dense_of_mat w
       in
       let eigenvectors, eigenvalues = 
         Parallel.broadcast (lazy (
@@ -318,7 +323,6 @@ TODO SINGLE STATE HERE
         eigenvectors, eigenvalues
     in
     iteration ci_coef
-*)
 
   )
   in