CI/F12CI.ml

@@ -130,11 +130,6 @@ 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
 
 
@@ -317,7 +312,6 @@ 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
@@ -334,112 +328,40 @@ 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 diagonal =
-      let eigenvectors, eigenvalues = 
+        Vec.init (Matrix.dim1 m_H) (fun i ->
-
+          if i = column_idx then
-(* Column dressing 
+            Matrix.get m_H i i +. delta.{column_idx,state} *. f 
-*)
+          else
-        let delta = lacpy delta in
+            Matrix.get m_H i i
-        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}
-            else
-              Matrix.get m_H i i
-            )
-        in
-
-        let matrix_prod c = 
-          let w = 
-            Matrix.mm ~transa:`T m_H c
-            |> Matrix.to_mat
-          in
-          let c = Matrix.to_mat c in
-
-          for k=1 to state do
-            for i=1 to (Mat.dim1 w) do
-              w.{i,k} <- w.{i,k} +.  delta.{i,k} *. c.{column_idx, k} ;
-              w.{column_idx,k} <- w.{column_idx,k} +.  delta.{i,k} *. c.{i,k};
-            done;
-            w.{column_idx,k} <- w.{column_idx,k} -.
-                delta.{column_idx,k} *. c.{column_idx,k};
-          done;
-
-          Matrix.dense_of_mat w
-        in
-
-
-(* Diagonal dressing *)
-(*
-        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
 
+      let matrix_prod c = 
+        let w = 
+          Matrix.mm ~transa:`T m_H c
+          |> Matrix.to_mat
+        in
+        let c11 = Matrix.get c column_idx state in
+        Util.list_range 1 (Mat.dim1 w)
+        |> List.iter (fun i ->
+          let dci = 
+            delta.{i,state} *. f ;
+          in
+          w.{i,state} <- w.{i,state} +. dci *.  c11;
+          if (i <> column_idx) then
+            w.{column_idx,state} <- w.{column_idx,state} +. dci *. (Matrix.get c i state);
+          );
+        Matrix.dense_of_mat w
+      in
 
+      let eigenvectors, eigenvalues = 
+        Parallel.broadcast (lazy (
+          Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states:state diagonal matrix_prod
+          ))
+      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=8)
+    ?(n_iter=10)
-    ?(threshold=1.e-7)
+    ?(threshold=1.e-6)
     diagonal
     matrix_prod
   =
@@ -105,14 +105,13 @@ let make
     (* Compute the residual as proposed new vectors *)
     let u_proposed = 
       Mat.init_cols n m (fun i k ->
-          let delta = lambda.{k} -. diagonal.{i} in
+          let delta = diagonal.{i} -. lambda.{k} in
           let delta = 
-            if abs_float delta > 1.e-2  then delta
+            if abs_float delta > 0.01 then delta
-            else if delta > 0. then 1.e-2
+            else if delta < 0. then -.0.01
-            else (-1.e-2)
+                 else 0.01
           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
 
@@ -131,12 +130,11 @@ 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.001 in
+    let thr = maxu *. 0.01 in
     let u_proposed = 
       Mat.map (fun x -> if abs_float x < thr then 0. else x) u_proposed
       |> Mat.to_col_vecs_list
@@ -157,6 +155,6 @@ let make
       (m_new_U |> pick_new |> Mat.of_col_vecs_list), lambda
 
   in
-  iteration [] u_new [] 1 30
+  iteration [] u_new [] 1 5