LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 12:23:31 +01:00
Code Releases Activity

Searching bug in CAS+EN2

Browse Source
This commit is contained in:
Anthony Scemama 2019-03-23 01:06:38 +01:00
parent 2735688efe
commit 494b0a2e3c
2 changed files with 31 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
CI/CI.ml
View File

@ -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

13
CI/F12CI.ml
View File

@ -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