Davidson F12CI OK

CI/F12CI.ml

@@ -241,47 +241,7 @@ 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
-      Util.list_range 1 (Mat.dim1 psi)
-      |> List.iter (fun i ->
-        let delta = (Matrix.get 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
-      in
-      let conv =
-        1.0 -. abs_float ( dot
-          (Mat.to_col_vecs psi).(0)
-          (Mat.to_col_vecs eigenvectors).(0) )
-      in
-      Printf.printf "Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift);
-      if conv > threshold then
-        iteration eigenvectors
-      else
-        let eigenvalues =
-          Vec.map (fun x -> x +. e_shift) eigenvalues
-        in
-        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 +.
@@ -289,19 +249,38 @@ TODO SINGLE STATE HERE
       in
       let delta = Matrix.to_mat delta in
       delta.{1,1} <- delta.{1,1} -. delta_00;
+
+(*------
+TODO SINGLE STATE HERE
+*)
+      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 delta.{1,1} *. psi.{1,1} else 0.) )
+        Vec.init (Matrix.dim1 m_H) (fun i ->
+          if i = 1 then
+            Matrix.get m_H i i +. delta.{1,1} *. f 
+          else
+            Matrix.get m_H i i
+          )
       in
-      let delta = Matrix.dense_of_mat delta in
       let matrix_prod c = 
-        Matrix.add
+        let w = 
-          (Matrix.mm ~transa:`T m_H c)
+          Matrix.mm ~transa:`T m_H c
-          delta
+          |> Matrix.to_mat
+        in
+        let c11 = Matrix.get c 1 1 in
+        Util.list_range 1 (Mat.dim1 w)
+        |> List.iter (fun i ->
+          let dci = 
+            delta.{i,1} *. f ;
+          in
+          w.{i,1} <- w.{i,1} +. dci *.  c11;
+          if (i <> 1) then
+            w.{1,1} <- w.{1,1} +. dci *. (Matrix.get c i 1);
+          );
+        Matrix.dense_of_mat w
       in
       let eigenvectors, eigenvalues = 
-        Parallel.broadcast (lazy (
+          Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
-            Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
-          ))
       in
       let conv =
         1.0 -. abs_float ( dot
@@ -318,7 +297,6 @@ TODO SINGLE STATE HERE
         eigenvectors, eigenvalues
     in
     iteration ci_coef
-*)
 
   )
   in