From 494b0a2e3cad3e704e691c1ad3f077d283a9a328 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 23 Mar 2019 01:06:38 +0100 Subject: [PATCH] Searching bug in CAS+EN2 --- CI/CI.ml | 55 ++++++++++++++++++++++++++--------------------------- CI/F12CI.ml | 13 ++++--------- 2 files changed, 31 insertions(+), 37 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index ea31578..f8fa0a7 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -414,14 +414,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } list_holes2 list_particles2 is_internal i_o1_alfa alfa_o2_i w_alfa psi0 = - let list_holes1 = Array.of_list list_holes1 - and list_holes2 = Array.of_list list_holes2 - and list_particles1 = Array.of_list list_particles1 - and list_particles2 = Array.of_list list_particles2 - in + let list_holes1 = Array.of_list list_holes1 in + let list_particles1 = Array.of_list list_particles1 in + let list_holes2 = Array.of_list list_holes2 in + let list_particles2 = Array.of_list list_particles2 in let psi0 = - let stream = Ds.determinant_stream det_space in @@ -448,7 +446,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } || aux (j-1) in aux (i-1) - ) else + ) + else is_internal in @@ -513,18 +512,17 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let double = Array.fold_left (fun accu particle' -> - if particle' > particle || particle' = hole then + if particle' >= particle || particle' = hole then accu else accu +. Array.fold_left (fun accu hole' -> - if hole' = particle' || hole' = particle || hole' < hole then + if hole' = particle' || hole' = particle || hole' <= hole then accu else let alfa = - Determinant.double_excitation - spin hole particle - spin hole' particle' det_i + Determinant.single_excitation + spin hole' particle' alfa in if Determinant.is_none alfa || already_generated alfa then @@ -550,25 +548,26 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } in if Determinant.is_none alfa then accu else - let double = + let double_ab = Array.fold_left (fun accu particle' -> - accu +. - Array.fold_left (fun accu hole' -> - if hole' = particle' then accu else - let alfa = - Determinant.double_excitation - Spin.Alfa hole particle - Spin.Beta hole' particle' det_i - in - if Determinant.is_none alfa || - already_generated alfa then - accu - else - accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa + accu +. + Array.fold_left (fun accu hole' -> + if hole' = particle' then accu else + let alfa = + Determinant.double_excitation + Spin.Alfa hole particle + Spin.Beta hole' particle' det_i + in + if Determinant.is_none alfa || + already_generated alfa then + accu + else + accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa ) 0. list_holes1 ) 0. list_particles1 in - accu +. double + + accu +. double_ab ) 0. list_holes2 ) 0. list_particles2 in @@ -757,7 +756,7 @@ let pt2_en_reference ci = in Array.map (fun ki -> Array.mapi (fun j kj -> - (h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j)) + (h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j)) ) out_dets ) in_dets |> Mat.of_array diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 051ccf2..f81f17e 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -197,7 +197,6 @@ debug "Four-idx transform of f12 intergals"; else 1. +. gamma *. (f_ij aux_basis ki kj) in -debug "Computing F matrix"; let m_F = CI.create_matrix_spin f det_space |> Lazy.force @@ -228,18 +227,10 @@ Util.debug_matrix "H" (Matrix.to_mat m_H); *) let rec iteration ?(state=1) psi = -(* debug "Iteration"; -Util.debug_matrix "T" (f12_amplitudes psi); -*) let delta = dressing_vector aux_basis (f12_amplitudes psi) ci in - (* -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 = @@ -262,6 +253,7 @@ TODO SINGLE STATE HERE Matrix.get m_H i i ) in + let matrix_prod c = let w = Matrix.mm ~transa:`T m_H c @@ -279,15 +271,18 @@ TODO SINGLE STATE HERE ); Matrix.dense_of_mat w in + let eigenvectors, eigenvalues = Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod 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