From 32c7c4e6c4e62ee08c7d90787bb6ad336322aa46 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Feb 2018 02:02:53 +0100 Subject: [PATCH] Loop over m --- Basis/TwoElectronRR.ml | 37 +++++--------- Basis/TwoElectronRRVectorized.ml | 88 ++++++++++++++++++-------------- 2 files changed, 63 insertions(+), 62 deletions(-) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 6aefaa2..1f3a498 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -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) and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab) in - let result = Array.create_float maxsze in + let result = Array.create_float (maxsze-angMom_a.tot) in if amxyz = 0 then begin - let v1 = - vrr0 am - in - for m=0 to maxm-1 do - result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) - done; - result.(maxm) <- -. f2 *. v1.(maxm) + let v1 = vrr0 am in + Array.iteri (fun m _ -> + result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)) result end else begin @@ -78,11 +74,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let v3 = vrr0 amm in let v1 = vrr0 am 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) - +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) - done; - result.(maxm) <- f3 *. v3.(maxm) + +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) ) result end; 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 = expo_inv_q *. (Coordinate.get xyz center_pq) 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 begin 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 = vrr am cm in - for m=0 to maxm-1 do - result.(m) <- result.(m) -. f5 *. v5.(m+1) - done + Array.iteri (fun m _ -> + result.(m) <- result.(m) -. f5 *. v5.(m+1)) result end; if cmxyz > 0 then 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 vrr angMom_a cmm in - for m=0 to maxm-1 do + Array.iteri (fun m _ -> result.(m) <- result.(m) +. - f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) - done; - result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm) + f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) ) result end end; 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 = vrr angMom_a cm in - for m=0 to maxm-1 do - result.(m) <- result.(m) +. f1 *. v1.(m) -. f2 *. v1.(m+1) ; - done; - result.(maxm) <- result.(maxm) +. f1 *. v1.(maxm) ; + Array.iteri (fun m _ -> + result.(m) <- result.(m) +. f1 *. v1.(m) -. f2 *. v1.(m+1) ) result end; result in Zmap.add map_2d key result; diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 0564806..788e76c 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -26,6 +26,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) map_1d map_2d np nq = + let maxm = Array.length zero_m_array - 1 in + let get_xyz angMom = match angMom with | { y=0 ; z=0 ; _ } -> X @@ -34,74 +36,84 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in (** Vertical recurrence relations *) - let rec vrr0_v m angMom_a = + let rec vrr0_v angMom_a = match angMom_a.tot with - | 0 -> zero_m_array.(m) + | 0 -> zero_m_array | _ -> let key = Zkey.of_powers_three angMom_a in - try Zmap.find map_1d.(m) key with + try Zmap.find map_1d key with | Not_found -> let result = let xyz = get_xyz angMom_a in let am = Powers.decr xyz angMom_a 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 if abs_float cab >= cutoff then - let v0 = vrr0_v m am in - Array.iteri (fun l result_l -> - let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab - and v0_l = v0.(l) - in - Array.iteri (fun k v0_lk -> - result_l.(k) <- v0_lk *. f0) v0_l ) result + Array.iteri (fun m result_m -> + let v0 = v_am.(m) in + Array.iteri (fun l result_ml -> + let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab + and v0_l = v0.(l) + in + Array.iteri (fun k v0_lk -> + result_ml.(k) <- v0_lk *. f0) v0_l + ) result_m + ) result end; let amxyz = Powers.get xyz am in if amxyz < 1 then - Array.iteri (fun l result_l -> - let expo_inv_p_l = expo_inv_p.(l) - and center_pq_xyz_l = (center_pq xyz).(l) - and result_l = result.(l) - and p0_l = p0.(l) - in - Array.iteri (fun k p0_lk -> - result_l.(k) <- result_l.(k) - +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk - ) p0_l ) result + Array.iteri (fun l expo_inv_p_l -> + let center_pq_xyz_l = (center_pq xyz).(l) in + Array.iteri (fun m result_m -> + let result_ml = result_m.(l) in + let p0 = v_am.(m+1) in + let p0_l = p0.(l) + in + Array.iteri (fun k p0_lk -> + result_ml.(k) <- result_ml.(k) + +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk + ) p0_l + ) result + ) expo_inv_p else begin 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 - 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 - and expo_inv_p_l = expo_inv_p.(l) and center_pq_xyz_l = (center_pq xyz).(l) - and v1_l = v1.(l) - and v2_l = v2.(l) - and result_l = result.(l) in - Array.iteri (fun k p0_lk -> - result_l.(k) <- result_l.(k) +. - expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +. - f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) - ) p0.(l) - ) result + Array.iteri (fun m result_m -> + 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 + Array.iteri (fun k p0_lk -> + result_ml.(k) <- result_ml.(k) +. + expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +. + f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) + ) p0.(l) + ) result + ) expo_inv_p end; result in - Zmap.add map_1d.(m) key result; + Zmap.add map_1d key result; result and vrr_v m angMom_a angMom_c = 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 @@ -783,7 +795,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q |> Array.concat 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)) in (* Compute the integral class from the primitive shell quartet *)