LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-07 06:33:39 +01:00
Code Releases Activity

Accelerated scalar RR

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-11 23:41:18 +01:00
parent 8c35202b29
commit 264181ebd1
3 changed files with 65 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Basis/ContractedShellPair.ml
View File

@ -16,7 +16,8 @@ type t =
center_ab : Coordinate.t; (* A-B *) center_ab : Coordinate.t; (* A-B *)
norm_sq : float; (* |A-B|^2 *) norm_sq : float; (* |A-B|^2 *)
norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *) norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *)
totAngMomInt : int (* Total angular Momentum *) totAngMomInt : int; (* Total angular Momentum *)
monocentric : bool;
} }
@ -119,6 +120,7 @@ let create ?cutoff p_a p_b =
shell_pairs ; center_ab=shell_pairs.(0).center_ab; shell_pairs ; center_ab=shell_pairs.(0).center_ab;
norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq; norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq;
totAngMomInt = shell_pairs.(0).ShellPair.totAngMomInt; totAngMomInt = shell_pairs.(0).ShellPair.totAngMomInt;
monocentric = shell_pairs.(0).ShellPair.monocentric;
} }

103
Basis/TwoElectronRR.ml
View File

@ -86,10 +86,6 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)
done; done;
result.(maxm) <- -. f2 *. v1.(maxm) result.(maxm) <- -. f2 *. v1.(maxm)
(*
Array.init maxsze (fun m ->
if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
*)
end end
else else
begin begin
@ -107,14 +103,6 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
result.(maxm) <- f3 *. v3.(maxm) result.(maxm) <- f3 *. v3.(maxm)
end; end;
result result
(*
Array.init maxsze (fun m ->
(if m = maxm then 0. else
(f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
+. f3 *. (v3.(m) +. if m = maxm then 0. else
expo_inv_p *. v3.(m+1))
)
*)
in Zmap.add map_1d key result; in Zmap.add map_1d key result;
result result
@ -164,47 +152,55 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
(angMom_cx, angMom_cy, angMom_cz-2), (angMom_cx, angMom_cy, angMom_cz-2),
angMom_az,angMom_cz-1, 2 angMom_az,angMom_cz-1, 2
in in
if cmxyz < 0 then empty else if cmxyz < 0 then empty
else
let f1 = let f1 =
-. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) -. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz)
in and f2 =
let f2 =
expo_inv_q *. (Coordinate.coord center_pq xyz) expo_inv_q *. (Coordinate.coord center_pq xyz)
in in
let result = let result = Array.make maxsze 0. in
if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then empty else if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then
begin
let v1 = let v1 =
vrr angMom_a cm totAngMom_a (totAngMom_c-1) vrr angMom_a cm totAngMom_a (totAngMom_c-1)
in in
Array.init maxsze (fun m -> f1 *. v1.(m) -. for m=0 to maxm-1 do
(if m = maxm then 0. else f2 *. v1.(m+1)) ) result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ;
in done;
let result = result.(maxm) <- f1 *. v1.(maxm) ;
if cmxyz < 1 then result else end;
if cmxyz > 0 then
begin
let f3 = let f3 =
(float_of_int cmxyz) *. expo_inv_q *. 0.5 (float_of_int cmxyz) *. expo_inv_q *. 0.5
in in
if (abs_float f3 < cutoff) && (abs_float (f3 *. expo_inv_q) < cutoff) then result else if (abs_float f3 > cutoff) ||
( (abs_float (f3 *. expo_inv_q) > cutoff) then
begin
let v3 = let v3 =
vrr angMom_a cmm totAngMom_a (totAngMom_c-2) vrr angMom_a cmm totAngMom_a (totAngMom_c-2)
in in
Array.init maxsze (fun m -> result.(m) +. for m=0 to maxm-1 do
f3 *. (v3.(m) +. (if m=maxm then 0. else expo_inv_q *. v3.(m+1)) )) result.(m) <- result.(m) +.
) f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1))
in done;
let result = result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm)
if (axyz < 1) || (cmxyz < 0) then result else end
end;
if (axyz > 0) && (cmxyz >= 0) then
begin
let f5 = let f5 =
(float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5
in in
if (abs_float f5 < cutoff) then result else if (abs_float f5 > cutoff) then
let v5 = let v5 =
vrr am cm (totAngMom_a-1) (totAngMom_c-1) vrr am cm (totAngMom_a-1) (totAngMom_c-1)
in in
Array.init (maxsze) (fun m -> for m=0 to maxm-1 do
result.(m) -. (if m = maxm then 0. else f5 *. v5.(m+1))) result.(m) <- result.(m) -. f5 *. v5.(m+1)
in done
end;
result result
in Zmap.add map_2d key result; in Zmap.add map_2d key result;
result result
@ -352,6 +348,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
Array.make (Array.length class_indices) 0.; Array.make (Array.length class_indices) 0.;
in in
let monocentric =
shell_p.ContractedShellPair.monocentric &&
shell_q.ContractedShellPair.monocentric
in
(* Compute all integrals in the shell for each pair of significant shell pairs *) (* Compute all integrals in the shell for each pair of significant shell pairs *)
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
@ -365,7 +366,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
in in
(** Screening on the product of coefficients *) (** Screening on the product of coefficients *)
try try
if (abs_float coef_prod) < 1.e-4*.cutoff then if (abs_float coef_prod) < 1.e-3*.cutoff then
raise NullQuartet; raise NullQuartet;
let expo_pq_inv = let expo_pq_inv =
@ -397,18 +398,19 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let map_1d' = Zmap.create (4*maxm) in let map_1d' = Zmap.create (4*maxm) in
let map_2d' = Zmap.create (Array.length class_indices) in let map_2d' = Zmap.create (Array.length class_indices) in
let norm_coef_scale = let norm_coef_scale =
Array.map (fun v1 -> Array.to_list norm_coef_scale_p
Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q |> List.map (fun v1 ->
) norm_coef_scale_p Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q)
|> Array.to_list
|> Array.concat |> Array.concat
in in
(* let () =
let monocentric = if debug then (
shell_p.(ab).ShellPair.monocentric && if monocentric then
shell_q.(cd).ShellPair.monocentric Printf.printf "Mono-centric\n"
else
Printf.printf "Multi-centric\n"
)
in in
*)
(* Compute the integral class from the primitive shell quartet *) (* Compute the integral class from the primitive shell quartet *)
class_indices class_indices
|> Array.iteri (fun i key -> |> Array.iteri (fun i key ->
@ -418,16 +420,19 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
| _ -> assert false | _ -> assert false
in in
try try
(*
if monocentric then if monocentric then
begin begin
if ( ((1 land (a.(0) + a.(3) + a.(6) + a.( 9)))=1) || let ax,ay,az = angMomA
((1 land (a.(1) + a.(4) + a.(7) + a.(10)))=1) || and bx,by,bz = angMomB
((1 land (a.(2) + a.(5) + a.(8) + a.(11)))=1) and cx,cy,cz = angMomC
and dx,dy,dz = angMomD
in
if ( ((1 land ax+bx+cx+dx)=1) ||
((1 land ay+by+cy+dy)=1) ||
((1 land az+bz+cz+dz)=1)
) then ) then
raise NullQuartet raise NullQuartet
end; end;
*)
(* Schwartz screening *) (* Schwartz screening *)
(* (*
let schwartz_p = let schwartz_p =

29
Basis/TwoElectronRRVectorized.ml
View File

@ -522,24 +522,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1) build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1)
in 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 = let np, nq =
Array.length p_list, Array.length p_list,
Array.length q_list Array.length q_list
@ -675,11 +657,12 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
in in
let norm = let norm =
let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in let norm_coef_scale_q =
Array.map (fun v1 -> shell_q.ContractedShellPair.norm_coef_scale
Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q in
) norm_coef_scale_p Array.to_list norm_coef_scale_p
|> Array.to_list |> List.map (fun v1 ->
Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q)
|> Array.concat |> Array.concat
in in