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

Filtered shell quartets

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-11 02:30:19 +01:00
parent 92596d872b
commit e72a88b106
2 changed files with 78 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Basis/TwoElectronRR.ml
View File

@ -129,17 +129,17 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let angMom_ax, angMom_ay, angMom_az = angMom_a
and angMom_cx, angMom_cy, angMom_cz = angMom_c in
match angMom_c with
| (_,0,0) -> (* 321_984 *)
| (_,0,0) ->
(angMom_ax-1, angMom_ay, angMom_az),
(angMom_cx-1, angMom_cy, angMom_cz),
(angMom_cx-2, angMom_cy, angMom_cz),
angMom_ax,angMom_cx-1, 0
| (_,_,0) -> (* 612_002 *)
| (_,_,0) ->
(angMom_ax, angMom_ay-1, angMom_az),
(angMom_cx, angMom_cy-1, angMom_cz),
(angMom_cx, angMom_cy-2, angMom_cz),
angMom_ay,angMom_cy-1, 1
| _ -> (* 1_067_324 *)
| _ ->
(angMom_ax, angMom_ay, angMom_az-1),
(angMom_cx, angMom_cy, angMom_cz-1),
(angMom_cx, angMom_cy, angMom_cz-2),
@ -152,14 +152,13 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let f2 =
expo_inv_q *. (Coordinate.coord center_pq xyz)
in
let result = empty in
let result =
if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then result else
if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then empty else
let v1 =
vrr angMom_a cm totAngMom_a (totAngMom_c-1)
in
Array.init maxsze (fun m -> result.(m) +.
f1 *. v1.(m) -. (if m = maxm then 0. else f2 *. v1.(m+1)) )
Array.init maxsze (fun m -> f1 *. v1.(m) -.
(if m = maxm then 0. else f2 *. v1.(m+1)) )
in
let result =
if cmxyz < 1 then result else

93
Basis/TwoElectronRRVectorized.ml
View File

@ -497,6 +497,61 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
Array.make (Array.length class_indices) 0.;
in
(** Screening on the product of coefficients *)
let coef_max_p =
Array.fold_left (fun accu x ->
if (abs_float x) > accu then (abs_float x) else accu)
0. shell_p.ContractedShellPair.coef
and coef_max_q =
Array.fold_left (fun accu x ->
if (abs_float x) > accu then (abs_float x) else accu)
0. shell_q.ContractedShellPair.coef
in
let rec build_list cutoff vec accu = function
| -1 -> Array.of_list accu
| k -> build_list cutoff vec (
if (abs_float vec.(k) > cutoff) then (k::accu)
else accu ) (k-1)
in
let p_list =
let vec = shell_p.ContractedShellPair.coef in
build_list (cutoff /. coef_max_q) vec [] (Array.length vec - 1)
and q_list =
let vec = shell_q.ContractedShellPair.coef in
build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1)
in
(*
if (Array.length p_list) < (Array.length shell_p.ContractedShellPair.coef) then
begin
Printf.printf "Reduced p from %d to %d:\n" (Array.length
shell_p.ContractedShellPair.coef) (Array.length p_list);
Array.iter (fun k -> Printf.printf "%d " k) p_list; print_newline ();
Printf.printf "\n%!"
end;
if (Array.length q_list) < (Array.length shell_q.ContractedShellPair.coef) then
begin
Printf.printf "Reduced q from %d to %d:\n" (Array.length
shell_q.ContractedShellPair.coef) (Array.length q_list);
Array.iter (fun k -> Printf.printf "%d " k) q_list; print_newline ();
Printf.printf "\n%!"
end;
*)
let np, nq =
Array.length p_list,
Array.length q_list
in
let filter_p vec = Array.init np (fun k -> vec.(p_list.(k)))
and filter_q vec = Array.init nq (fun k -> vec.(q_list.(k)))
in
let sp = filter_p sp
and sq = filter_q sq
in
(* Compute all integrals in the shell for each pair of significant shell pairs *)
begin
@ -504,9 +559,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
totAngMom shell_c, totAngMom shell_d) with
| Angular_momentum.(S,S,S,S) ->
contracted_class.(0) <-
let np, nq =
Array.length sp, Array.length sq
in
begin
try
let expo_inv_p =
Vec.init np (fun ab -> sp.(ab-1).ShellPair.expo_inv)
and expo_inv_q =
@ -516,13 +570,12 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let coef =
let result = Mat.make0 nq np in
Lacaml.D.ger
(Vec.of_array shell_q.ContractedShellPair.coef)
(Vec.of_array shell_p.ContractedShellPair.coef)
(Vec.of_array @@ filter_q shell_q.ContractedShellPair.coef)
(Vec.of_array @@ filter_p shell_p.ContractedShellPair.coef)
result;
result
in
let zm_array = Mat.init_cols np nq (fun i j ->
(** Screening on the product of coefficients *)
try
if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then
raise NullQuartet;
@ -544,15 +597,23 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
with NullQuartet -> 0.
) in
Mat.gemm_trace zm_array coef
with (Invalid_argument _) -> 0.
end
| _ ->
let coef =
let cp = filter_p shell_p.ContractedShellPair.coef
and cq = filter_q shell_q.ContractedShellPair.coef
in
Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) )
in
let expo_inv_p =
Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp
and expo_inv_q =
Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq
in
let np, nq =
Array.length expo_inv_p, Array.length expo_inv_q
in
let expo_b =
Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp
and expo_d =
@ -560,14 +621,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
in
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
let coef =
Array.init np (fun l ->
Array.init nq (fun k ->
shell_q.ContractedShellPair.coef.(k) *.
shell_p.ContractedShellPair.coef.(l)
) )
in
let center_pq =
Array.init 3 (fun xyz ->
Array.init np (fun ab ->
@ -610,12 +663,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
) sq
(*
|> Array.to_list
|> List.filter (fun (zero_m_array, d,
center_pq,coef_prod) -> abs_float coef_prod >= 1.e-4 *. cutoff)
|> Array.of_list
*)
in
(* Transpose result *)
for m=0 to maxm do