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

Reduced RAM in F12

This commit is contained in:
Anthony Scemama 2019-03-29 19:01:51 +01:00
parent 33387b6c24
commit 1896d9c45f
2 changed files with 27 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
CI/DeterminantSpace.ml
View File

@ -270,11 +270,6 @@ let arbitrary_of_mo_basis mo_basis f =
let n_det_beta = Array.length det_beta in
let n_det_alfa = Array.length det_alfa in
let ndet = n_det_alfa * n_det_beta in
if Parallel.master then
Format.printf "Number of determinants : %d %d %d\n%!"
n_det_alfa n_det_beta ndet;
let det = Array.make n_det_alfa
(Array.init n_det_beta (fun i -> i))
in

161
CI/F12CI.ml
View File

@ -86,76 +86,41 @@ let is_internal det_space =
| 1 | 2 -> false
| _ -> true
(*
if not (Bitstring.logand aux_mask alfa |> Bitstring.is_zero ) then
false
else
let beta =
Determinant.beta a
|> Spindeterminant.bitstring
in
Bitstring.logand aux_mask beta |> Bitstring.is_zero
*)
let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
(*
let i_o1_alfa = h_ij aux_basis in
let alfa_o2_i = f_ij aux_basis in
let w_alfa _ _ = 1. in
let mo_class = CI.mo_class ci in
let list_holes = List.concat
[ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
and list_particles2 = MOClass.auxiliary_mos mo_class
and list_particles1 = List.concat
[ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ; MOClass.auxiliary_mos mo_class ]
in
(* Single state here *)
let result =
CI.second_order_sum ci list_holes list_particles1 list_holes list_particles2
(is_internal ci.det_space) i_o1_alfa alfa_o2_i w_alfa f12_amplitudes ~unique:false
|> Vec.of_list
in
Matrix.sparse_of_vector_array [| Vector.sparse_of_vec result |]
*)
Printf.printf "Building matrix\n%!";
if Parallel.master then
Printf.printf "Building matrix\n%!";
(* Determinants of the FCI space as a list *)
let in_dets =
DeterminantSpace.determinants_array ci.CI.det_space
|> Array.to_list
DeterminantSpace.determinant_stream ci.CI.det_space
|> Util.stream_to_list
in
(* Stream that generates only singly and doubly excited determinants
wrt FCI space *)
let out_dets_stream =
(* Stream that generates all determinants of FCI space *)
let s =
DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis
|> DeterminantSpace.determinant_stream
(* Stream that generates all determinants of FCI space *)
let s =
DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis
|> DeterminantSpace.determinant_stream
in
(* Select only singly and 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
result ()
else
Some ki
in
result ()
with Stream.Failure -> None
)
in
(* Select only singly and 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
result ()
else
Some ki
in
result ()
with Stream.Failure -> None
)
in
let make_h_and_f alpha_list =
@ -187,9 +152,8 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
in
Printf.printf "Matrix product\n%!";
let m_HF =
let batch_size = 1 + 10_000_000 / (Mat.dim1 f12_amplitudes) in
let batch_size = 1 + 1_000_000 / (Mat.dim1 f12_amplitudes) in
let input_stream =
Stream.from (fun i ->
let rec make_batch accu = function
@ -213,7 +177,7 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
in
let iteration input =
Printf.printf "gemm\n%!";
Printf.printf ".%!";
let m_H_aux, m_F_aux = make_h_and_f input in
let m_HF =
gemm m_H_aux m_F_aux ~transb:`T
@ -225,44 +189,10 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
|> Farm.run ~ordered:false ~f:iteration
|> Stream.iter (fun hf ->
ignore @@ Mat.add result hf ~c:result );
result
Printf.printf "\n";
Parallel.broadcast (lazy result)
in
(*
let m_HF =
let in_dets =
DeterminantSpace.determinants_array ci.CI.det_space
in
let fci_space = DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis in
Array.mapi (fun i ki ->
Printf.printf "%d / %d\r%!" i (Array.length in_dets);
Array.map (fun kj ->
match Determinant.degrees ki kj with
| (0,0) | (0,1) | (0,2) | (0,3) | (0,4)
| (1,0) | (2,0) | (3,0) | (4,0)
| (1,1) | (2,1) | (3,1)
| (1,2) | (1,3)
| (2,2) ->
( let x = ref 0. in
DeterminantSpace.determinant_stream fci_space
|> Stream.iter (fun k_alfa ->
if not (is_internal ci.CI.det_space ki) then
let f = f_ij aux_basis k_alfa kj in
if f <> 0. then
let h = h_ij aux_basis ki k_alfa in
x := !x +. f *. h
);
!x )
| _ -> 0.
) in_dets
) in_dets
|> Mat.of_array
in
*)
Printf.printf "Done\n%!";
Matrix.dense_of_mat m_HF
@ -316,43 +246,6 @@ Printf.printf "det space\n%!";
Parallel.broadcast (lazy x)
in
(*
let f12_amplitudes =
fun c ->
let result = lacpy c in
Mat.scal (0.5) result ;
result
(* While in a sequential region, initiate the parallel
4-idx transformation to avoid nested parallel jobs
*)
ignore @@ MOBasis.f12_ints aux_basis;
let two_gamma_inv = -2. *. f12.F12.expo_s in
let f = fun ki kj ->
if ki <> kj then
(f_ij aux_basis ki kj)
else
two_gamma_inv +. (f_ij aux_basis ki kj)
in
let m_F =
CI.create_matrix_spin f det_space
|> Lazy.force
in
fun ci_coef ->
let result =
Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef)
|> Matrix.to_mat
in
let norm = gemm ~transa:`T result result in
Printf.printf "Norm of |F> : %f\n%!" norm.{1,1};
if abs_float norm.{1,1} > 1. then
failwith "Norm of |F> > 1";
result
in
*)
let e_shift =
let det =
DeterminantSpace.determinant_stream det_space