LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

F12 without singles

This commit is contained in:
Anthony Scemama 2019-05-29 09:04:43 +02:00
parent f88409bfbd
commit 55b58268f8
1 changed files with 9 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
CI/F12CI.ml
View File

@ -61,7 +61,7 @@ let hf_ij mo_basis ki kj =
[ h_ij mo_basis ki kj ; f_ij mo_basis ki kj ]
let is_internal det_space =
let is_a_double det_space =
let mo_class = DeterminantSpace.mo_class det_space in
let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
let m l =
@ -83,8 +83,8 @@ let is_internal det_space =
and b = Bitstring.logand aux_mask beta
in
match Bitstring.popcount a + Bitstring.popcount b with
| 1 | 2 -> false
| _ -> true
| 2 -> true
| _ -> false
let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
@ -107,13 +107,12 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
|> DeterminantSpace.determinant_stream
in
(* Select only singly and doubly excited determinants
wrt FCI space *)
(* Select only doubly excited determinants wrt FCI space *)
Stream.from (fun _ ->
try
let rec result () =
let ki = Stream.next s in
if is_internal ci.CI.det_space ki then
if not (is_a_double ci.CI.det_space ki) then
result ()
else
Some ki
@ -227,15 +226,12 @@ Printf.printf "Add aux basis\n%!";
MOBasis.of_mo_basis s mo_basis
in
let () =
Printf.printf "F12 ints\n%!";
ignore @@ MOBasis.f12_ints aux_basis
in
let () =
Printf.printf "2e ints\n%!";
ignore @@ MOBasis.two_e_ints aux_basis
in
Printf.printf "det space\n%!";
let det_space =
DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num
in
@ -269,12 +265,16 @@ Printf.printf "det space\n%!";
|> Matrix.to_mat
in
Printf.printf "Cmax : %e\n" psi.{column_idx,state};
Printf.printf "Norm : %e\n" (sqrt (gemm ~transa:`T delta delta).{state,state});
let f = 1.0 /. psi.{column_idx,state} in
let delta_00 =
(* Delta_00 = {% $\sum_{j \ne x} delta_j c_j / c_x$ %} *)
f *. ( (gemm ~transa:`T delta psi).{state,state} -.
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 =