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

Less allocations

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-10 22:19:13 +01:00
parent 943dfa9935
commit 346db83439
1 changed files with 89 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

153
Basis/TwoElectronRRVectorized.ml
View File

@ -24,9 +24,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
map_1d map_2d np nq map_1d map_2d np nq
= =
let empty = let empty = Array.make nq 0. in
Array.make nq 0. let tmp_1 = Array.make nq 0. in
in
let totAngMom_a = Angular_momentum.to_int totAngMom_a let totAngMom_a = Angular_momentum.to_int totAngMom_a
@ -54,59 +53,78 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
if amxyz < 0 then if amxyz < 0 then
None None
else else
let v1_top, p1_top = begin
if abs_float (Coordinate.coord center_ab xyz) < cutoff then let cab = Coordinate.coord center_ab xyz in
None, let v1_top, p1_top =
vrr0_v (m+1) am (totAngMom_a-1) if abs_float cab < cutoff then
else None,
vrr0_v m am (totAngMom_a-1), vrr0_v (m+1) am (totAngMom_a-1)
vrr0_v (m+1) am (totAngMom_a-1) else
in vrr0_v m am (totAngMom_a-1),
let v1_top2, p1_top2 = vrr0_v (m+1) am (totAngMom_a-1)
if amxyz < 1 then (None,None) else in
vrr0_v m amm (totAngMom_a-2), let v1_top2, p1_top2 =
vrr0_v (m+1) amm (totAngMom_a-2) if amxyz < 1 then (None,None) else
in vrr0_v m amm (totAngMom_a-2),
vrr0_v (m+1) amm (totAngMom_a-2)
in
Some ( let result = Array.make_matrix np nq 0. in
Array.init np (fun l -> for l=0 to np-1
let v1 = do
let f = let p0 =
-. expo_b.(l) *. expo_inv_p.(l) *. (Coordinate.coord center_ab xyz) match p1_top with
in | Some p1_top -> p1_top.(l)
match v1_top with | _ -> assert false
| Some v1_top -> in
v1_top.(l) |> Array.map (fun x -> f *. x) let f0 =
| None -> empty -. expo_b.(l) *. expo_inv_p.(l) *. cab
in in
let p1 = let v0 =
match p1_top with match v1_top with
| Some p1_top -> p1_top.(l) | Some v1_top -> v1_top.(l)
| _ -> assert false | None -> empty
in in
let p1 = for k=0 to nq-1
Array.init nq (fun k -> do
v1.(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p1.(k)) tmp_1.(k) <- v0.(k) *. f0
in done;
if amxyz < 1 then p1 else let v0 = tmp_1 in
let f = (float_of_int amxyz) *. expo_inv_p.(l) *. 0.5 if amxyz < 1 then
in begin
let v1 = for k=0 to nq-1
match v1_top2 with do
| Some v1_top2 -> v1_top2.(l) result.(l).(k) <- v0.(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(k)
| None -> assert false done
in end
let v2 = else
match p1_top2 with begin
| Some p1_top2 -> p1_top2.(l) let f = (float_of_int amxyz) *. expo_inv_p.(l) *. 0.5
| None -> assert false in
in let v1 =
Array.init nq (fun k -> match v1_top2 with
p1.(k) +. f *. (v1.(k) +. v2.(k) *. expo_inv_p.(l) ) ) | Some v1_top2 -> v1_top2.(l)
) | None -> assert false
) in
in Zmap.add map_1d.(m) key result; let v2 =
result match p1_top2 with
| Some p1_top2 -> p1_top2.(l)
| None -> assert false
in
for k=0 to nq-1
do
result.(l).(k) <-
v0.(k)
+. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(k)
+. f *. (v1.(k) +. v2.(k) *. expo_inv_p.(l) )
done
end
done;
Some result
end
in
Zmap.add map_1d.(m) key result;
result
and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c = and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c =
@ -140,6 +158,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
(acx-2,acy,acz), (acx-2,acy,acz),
aax,acx, 0 aax,acx, 0
in in
(*
let result = Array.make_matrix np nq 0. in
*)
let f1 = let f1 =
let f = (Coordinate.coord center_cd xyz) in let f = (Coordinate.coord center_cd xyz) in
Array.init nq (fun k -> expo_d.(k) *. expo_inv_q.(k) *. f) Array.init nq (fun k -> expo_d.(k) *. expo_inv_q.(k) *. f)
@ -484,19 +505,23 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
Array.init (maxm+1) (fun _ -> Array.init (maxm+1) (fun _ ->
Array.init np (fun _ -> Array.make nq 0. ) ) Array.init np (fun _ -> Array.make nq 0. ) )
in in
let empty = Array.make (maxm+1) 0. in
Array.iteri (fun ab shell_ab -> Array.iteri (fun ab shell_ab ->
let zero_m_array_tmp = let zero_m_array_tmp =
Array.mapi (fun cd shell_cd -> Array.mapi (fun cd shell_cd ->
let expo_pq_inv = if (abs_float coef.(ab).(cd) < cutoff) then
expo_inv_p.(ab) +. expo_inv_q.(cd) empty
in else
let norm_pq_sq = let expo_pq_inv =
center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +. expo_inv_p.(ab) +. expo_inv_q.(cd)
center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +. in
center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd) let norm_pq_sq =
in center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +.
center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +.
center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd)
in
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
) sq ) sq
(* (*
|> Array.to_list |> Array.to_list