fixed davidson in f12

CI/F12CI.ml

@@ -130,6 +130,11 @@ let p12 det_space =
                   (Spindeterminant.of_bitstring  @@
                     Bitstring.(logor a (logand not_aux_mask beta))
                   ) )
+      | 1, 0 
+      | 0, 1 -> Some (Determinant.negate_phase k)
+(*
+      | 0, 1 -> Some (Determinant.vac 1)
+*)
       | _ -> None
 
 
@@ -312,6 +317,7 @@ Printf.printf "Add aux basis\n%!";
       
     let rec iteration ~state psi = 
       let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
+
       let delta = 
          (* delta_i = {% $\sum_j c_j H_{ij}$ %} *)
          dressing_vector ~frozen_core aux_basis psi ci
@@ -328,12 +334,25 @@ Printf.printf "Add aux basis\n%!";
           delta.{column_idx,state} *. psi.{column_idx,state} )
       in
       Printf.printf "Delta_00 : %e %e\n" delta.{column_idx,state} delta_00;
+
       delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00;
 
+
+      let eigenvectors, eigenvalues = 
+
+(* Column dressing 
+*)
+        let delta = lacpy delta in
+        Mat.scal f delta;
+        for k=1 to state-1 do
+          for i=1 to Mat.dim1 delta do
+                delta.{i,k} <- delta.{i,state}
+          done;
+        done;
         let diagonal =
           Vec.init (Matrix.dim1 m_H) (fun i ->
             if i = column_idx then
-            Matrix.get m_H i i +. delta.{column_idx,state} *. f 
+              Matrix.get m_H i i +. delta.{column_idx,state}
             else
               Matrix.get m_H i i
             )
@@ -344,24 +363,83 @@ Printf.printf "Add aux basis\n%!";
             Matrix.mm ~transa:`T m_H c
             |> Matrix.to_mat
           in
-        let c11 = Matrix.get c column_idx state in
+          let c = Matrix.to_mat c in
-        Util.list_range 1 (Mat.dim1 w)
+
-        |> List.iter (fun i ->
+          for k=1 to state do
-          let dci = 
+            for i=1 to (Mat.dim1 w) do
-            delta.{i,state} *. f ;
+              w.{i,k} <- w.{i,k} +.  delta.{i,k} *. c.{column_idx, k} ;
-          in
+              w.{column_idx,k} <- w.{column_idx,k} +.  delta.{i,k} *. c.{i,k};
-          w.{i,state} <- w.{i,state} +. dci *.  c11;
+            done;
-          if (i <> column_idx) then
+            w.{column_idx,k} <- w.{column_idx,k} -.
-            w.{column_idx,state} <- w.{column_idx,state} +. dci *. (Matrix.get c i state);
+                delta.{column_idx,k} *. c.{column_idx,k};
-          );
+          done;
+
           Matrix.dense_of_mat w
         in
 
-      let eigenvectors, eigenvalues = 
+
-        Parallel.broadcast (lazy (
+(* Diagonal dressing *)
-          Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states:state diagonal matrix_prod
+(*
-          ))
+        let diagonal =
+          Vec.init (Matrix.dim1 m_H) (fun i ->
+            Matrix.get m_H i i +.
+              if (abs_float psi.{i,state} > 1.e-8) then
+                delta.{i,state} /. psi.{i,state}
+              else  0.
+            )
         in
+
+        let matrix_prod c = 
+          let w = 
+            Matrix.mm ~transa:`T m_H c
+            |> Matrix.to_mat
+          in
+          for i=1 to (Mat.dim1 w) do
+            w.{i,state} <- w.{i,state} +. delta.{i,state}
+          done;
+          Matrix.dense_of_mat w
+        in
+*)
+
+
+        Parallel.broadcast (lazy (
+          Davidson.make ~threshold:1.e-9 ~guess:psi ~n_states:state diagonal matrix_prod
+          ))
+
+(*
+        let m_H = Matrix.to_mat m_H |> lacpy in
+*)
+
+(* DIAGONAL TEST 
+        for i=1 to Mat.dim1 m_H do
+          if (abs_float psi.{i,state} > 1.e-8) then
+            m_H.{i,i} <- m_H.{i,i} +. delta.{i,state} /. psi.{i,state};
+        done;
+*)
+
+
+(* COLUMN TEST
+        for i=1 to Mat.dim1 m_H do
+          let d = delta.{i,state} /. psi.{column_idx,state} in
+          m_H.{i,column_idx} <- m_H.{i,column_idx} +. d;
+          if (i <> column_idx) then
+            begin
+              m_H.{column_idx,i} <- m_H.{column_idx,i} +. d;
+              m_H.{column_idx,column_idx} <- m_H.{column_idx,column_idx} -.
+                d *. psi.{i,state}
+            end
+        done;
+*)
+
+
+
+(*
+        let m_v = syev m_H in
+        m_H, m_v
+*)
+      in
+
+
       Vec.iter (fun energy -> Printf.printf "%g\t" energy) eigenvalues;
       print_newline ();
 

Utils/Davidson.ml

@@ -5,8 +5,8 @@ type t
 let make
     ?guess
     ?(n_states=1)
-    ?(n_iter=10)
+    ?(n_iter=8)
-    ?(threshold=1.e-6)
+    ?(threshold=1.e-7)
     diagonal
     matrix_prod
   =
@@ -105,13 +105,14 @@ let make
     (* Compute the residual as proposed new vectors *)
     let u_proposed = 
       Mat.init_cols n m (fun i k ->
-          let delta = diagonal.{i} -. lambda.{k} in
+          let delta = lambda.{k} -. diagonal.{i} in
           let delta = 
-            if abs_float delta > 0.01 then delta
+            if abs_float delta > 1.e-2  then delta
-            else if delta < 0. then -.0.01
+            else if delta > 0. then 1.e-2
-                 else 0.01
+            else (-1.e-2)
           in
-        (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. delta )
+            (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. delta
+          )
       |> Mat.to_col_vecs_list
     in
 
@@ -130,11 +131,12 @@ let make
       end;
 
     (* Make new vectors sparse *)
+
     let u_proposed = 
       Mat.of_col_vecs_list u_proposed
     in
     let maxu = lange u_proposed ~norm:`M in
-    let thr = maxu *. 0.01 in
+    let thr = maxu *. 0.001 in
     let u_proposed = 
       Mat.map (fun x -> if abs_float x < thr then 0. else x) u_proposed
       |> Mat.to_col_vecs_list
@@ -155,6 +157,6 @@ let make
       (m_new_U |> pick_new |> Mat.of_col_vecs_list), lambda
 
   in
-  iteration [] u_new [] 1 5
+  iteration [] u_new [] 1 30