From 4cfb4cd241359c7cf8ce7c7133f1d3430fec60ae Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 26 Aug 2019 22:34:41 +0200 Subject: [PATCH] fixed davidson in f12 --- CI/F12CI.ml | 132 ++++++++++++++++++++++++++++++++++++---------- Utils/Davidson.ml | 20 +++---- 2 files changed, 116 insertions(+), 36 deletions(-) diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 44db115..ee6c923 100644 --- a/CI/F12CI.ml +++ b/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,40 +334,112 @@ 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 = - Vec.init (Matrix.dim1 m_H) (fun i -> - if i = column_idx then - Matrix.get m_H i i +. delta.{column_idx,state} *. f - 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 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 = + +(* 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} + 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-6 ~guess:psi ~n_states:state diagonal matrix_prod + 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 (); diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml index 0af4e15..d930ede 100644 --- a/Utils/Davidson.ml +++ b/Utils/Davidson.ml @@ -5,8 +5,8 @@ type t let make ?guess ?(n_states=1) - ?(n_iter=10) - ?(threshold=1.e-6) + ?(n_iter=8) + ?(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 - else if delta < 0. then -.0.01 - else 0.01 + if abs_float delta > 1.e-2 then delta + else if delta > 0. then 1.e-2 + 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