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

Loop over m

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-23 02:02:53 +01:00
parent fd1b05c4ec
commit 32c7c4e6c4
2 changed files with 63 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
Basis/TwoElectronRR.ml
View File

@ -61,16 +61,12 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) let f1 = expo_inv_p *. (Coordinate.get xyz center_pq)
and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab) and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab)
in in
let result = Array.create_float maxsze in let result = Array.create_float (maxsze-angMom_a.tot) in
if amxyz = 0 then if amxyz = 0 then
begin begin
let v1 = let v1 = vrr0 am in
vrr0 am Array.iteri (fun m _ ->
in result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)) result
for m=0 to maxm-1 do
result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)
done;
result.(maxm) <- -. f2 *. v1.(maxm)
end end
else else
begin begin
@ -78,11 +74,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let v3 = vrr0 amm in let v3 = vrr0 amm in
let v1 = vrr0 am in let v1 = vrr0 am in
let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
for m=0 to maxm-1 do Array.iteri (fun m _ ->
result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)
+. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) ) result
done;
result.(maxm) <- f3 *. v3.(maxm)
end; end;
result result
in Zmap.add map_1d key result; in Zmap.add map_1d key result;
@ -111,7 +105,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
and f2 = and f2 =
expo_inv_q *. (Coordinate.get xyz center_pq) expo_inv_q *. (Coordinate.get xyz center_pq)
in in
let result = Array.make maxsze 0. in let result = Array.make (maxsze - angMom_a.tot - angMom_c.tot) 0. in
if axyz > 0 then if axyz > 0 then
begin begin
let am = Powers.decr xyz angMom_a in let am = Powers.decr xyz angMom_a in
@ -122,9 +116,8 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let v5 = let v5 =
vrr am cm vrr am cm
in in
for m=0 to maxm-1 do Array.iteri (fun m _ ->
result.(m) <- result.(m) -. f5 *. v5.(m+1) result.(m) <- result.(m) -. f5 *. v5.(m+1)) result
done
end; end;
if cmxyz > 0 then if cmxyz > 0 then
begin begin
@ -138,11 +131,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let cmm = Powers.decr xyz cm in let cmm = Powers.decr xyz cm in
vrr angMom_a cmm vrr angMom_a cmm
in in
for m=0 to maxm-1 do Array.iteri (fun m _ ->
result.(m) <- result.(m) +. result.(m) <- result.(m) +.
f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) ) result
done;
result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm)
end end
end; end;
if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then
@ -150,10 +141,8 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let v1 = let v1 =
vrr angMom_a cm vrr angMom_a cm
in in
for m=0 to maxm-1 do Array.iteri (fun m _ ->
result.(m) <- result.(m) +. f1 *. v1.(m) -. f2 *. v1.(m+1) ; result.(m) <- result.(m) +. f1 *. v1.(m) -. f2 *. v1.(m+1) ) result
done;
result.(maxm) <- result.(maxm) +. f1 *. v1.(maxm) ;
end; end;
result result
in Zmap.add map_2d key result; in Zmap.add map_2d key result;

64
Basis/TwoElectronRRVectorized.ml
View File

@ -26,6 +26,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 maxm = Array.length zero_m_array - 1 in
let get_xyz angMom = let get_xyz angMom =
match angMom with match angMom with
| { y=0 ; z=0 ; _ } -> X | { y=0 ; z=0 ; _ } -> X
@ -34,74 +36,84 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
in in
(** Vertical recurrence relations *) (** Vertical recurrence relations *)
let rec vrr0_v m angMom_a = let rec vrr0_v angMom_a =
match angMom_a.tot with match angMom_a.tot with
| 0 -> zero_m_array.(m) | 0 -> zero_m_array
| _ -> | _ ->
let key = Zkey.of_powers_three angMom_a let key = Zkey.of_powers_three angMom_a
in in
try Zmap.find map_1d.(m) key with try Zmap.find map_1d key with
| Not_found -> | Not_found ->
let result = let result =
let xyz = get_xyz angMom_a in let xyz = get_xyz angMom_a in
let am = Powers.decr xyz angMom_a in let am = Powers.decr xyz angMom_a in
let cab = Coordinate.get xyz center_ab in let cab = Coordinate.get xyz center_ab in
let p0 = vrr0_v (m+1) am in let result = Array.init (maxm+1-angMom_a.tot) (fun _ -> Array.make_matrix np nq 0.) in
let v_am= vrr0_v am in
let result = Array.make_matrix np nq 0. in
begin begin
if abs_float cab >= cutoff then if abs_float cab >= cutoff then
let v0 = vrr0_v m am in Array.iteri (fun m result_m ->
Array.iteri (fun l result_l -> let v0 = v_am.(m) in
Array.iteri (fun l result_ml ->
let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab
and v0_l = v0.(l) and v0_l = v0.(l)
in in
Array.iteri (fun k v0_lk -> Array.iteri (fun k v0_lk ->
result_l.(k) <- v0_lk *. f0) v0_l ) result result_ml.(k) <- v0_lk *. f0) v0_l
) result_m
) result
end; end;
let amxyz = Powers.get xyz am in let amxyz = Powers.get xyz am in
if amxyz < 1 then if amxyz < 1 then
Array.iteri (fun l result_l -> Array.iteri (fun l expo_inv_p_l ->
let expo_inv_p_l = expo_inv_p.(l) let center_pq_xyz_l = (center_pq xyz).(l) in
and center_pq_xyz_l = (center_pq xyz).(l) Array.iteri (fun m result_m ->
and result_l = result.(l) let result_ml = result_m.(l) in
and p0_l = p0.(l) let p0 = v_am.(m+1) in
let p0_l = p0.(l)
in in
Array.iteri (fun k p0_lk -> Array.iteri (fun k p0_lk ->
result_l.(k) <- result_l.(k) result_ml.(k) <- result_ml.(k)
+. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk
) p0_l ) result ) p0_l
) result
) expo_inv_p
else else
begin begin
let amm = Powers.decr xyz am in let amm = Powers.decr xyz am in
let v1 = vrr0_v m amm in
let v2 = vrr0_v (m+1) amm in
let amxyz = float_of_int amxyz in let amxyz = float_of_int amxyz in
Array.iteri (fun l result_l -> let v_amm = vrr0_v amm in
Array.iteri (fun l expo_inv_p_l ->
let f = amxyz *. expo_inv_p.(l) *. 0.5 let f = amxyz *. expo_inv_p.(l) *. 0.5
and expo_inv_p_l = expo_inv_p.(l)
and center_pq_xyz_l = (center_pq xyz).(l) and center_pq_xyz_l = (center_pq xyz).(l)
and v1_l = v1.(l) in
and v2_l = v2.(l) Array.iteri (fun m result_m ->
and result_l = result.(l) let v1 = v_amm.(m) in
let v1_l = v1.(l) in
let result_ml = result_m.(l) in
let v2 = v_amm.(m+1) in
let p0 = v_am.(m+1) in
let v2_l = v2.(l)
in in
Array.iteri (fun k p0_lk -> Array.iteri (fun k p0_lk ->
result_l.(k) <- result_l.(k) +. result_ml.(k) <- result_ml.(k) +.
expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +.
f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l)
) p0.(l) ) p0.(l)
) result ) result
) expo_inv_p
end; end;
result result
in in
Zmap.add map_1d.(m) key result; Zmap.add map_1d key result;
result result
and vrr_v m angMom_a angMom_c = and vrr_v m angMom_a angMom_c =
match (angMom_a.tot, angMom_c.tot) with match (angMom_a.tot, angMom_c.tot) with
| (i,0) -> Some (vrr0_v m angMom_a) | (i,0) -> Some (vrr0_v angMom_a).(m)
| (_,_) -> | (_,_) ->
let key = Zkey.of_powers_six angMom_a angMom_c in let key = Zkey.of_powers_six angMom_a angMom_c in
@ -783,7 +795,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
|> Array.concat |> Array.concat
in in
let map_1d = Array.init (maxm+1) (fun _ -> Zmap.create (4*maxm)) let map_1d = Zmap.create (4*maxm)
and map_2d = Array.init (maxm+1) (fun _ -> Zmap.create (Array.length class_indices)) and map_2d = Array.init (maxm+1) (fun _ -> Zmap.create (Array.length class_indices))
in in
(* Compute the integral class from the primitive shell quartet *) (* Compute the integral class from the primitive shell quartet *)