diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 0a456d4..7f8c74d 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -3,9 +3,15 @@ open Util let cutoff2 = cutoff *. cutoff exception NullQuartet +exception Found + +let at_least_one_valid arr = + try + Array.fold_left (fun _ x -> if (abs_float x > cutoff) then raise Found else false ) false arr + with Found -> true (** Horizontal and Vertical Recurrence Relations (HVRR) *) -let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) +let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (totAngMom_a, totAngMom_b, totAngMom_c, totAngMom_d) (maxm, zero_m_array) (expo_b, expo_d) @@ -30,6 +36,7 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) | 1 -> let i = if angMom_a.(0) = 1 then 0 else if angMom_a.(1) = 1 then 1 else 2 in let f = expo_b *. (Coordinate.coord center_ab i) in + if (abs_float f < cutoff) then empty else Array.mapi (fun k c -> c *. expo_inv_p *. ( (Coordinate.coord center_pq.(k) i) *. zero_m_array.(k).(m+1) -. f *. zero_m_array.(k).(m) ) ) coef_prod @@ -70,7 +77,10 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) in if (abs_float f < cutoff) then empty else let v1 = vrr0_v m amm (totAngMom_a-2) - and v2 = vrr0_v (m+1) amm (totAngMom_a-2) + in + let v2 = + if (abs_float (f *. expo_inv_p)) < cutoff then empty else + vrr0_v (m+1) amm (totAngMom_a-2) in Array.mapi (fun k _ -> p1.(k) +. f *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) coef_prod @@ -111,43 +121,62 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) if cm.(xyz) < 0 then empty else + let f1 = + Array.mapi (fun k _ -> + expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) ) expo_inv_q + in + let f2 = + Array.mapi (fun k _ -> + expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz) ) expo_inv_q + in let v1 = - vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) + if (at_least_one_valid f1) then + vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) + else empty and v2 = - vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) + if (at_least_one_valid f2) then + vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) + else empty in let p1 = - Array.mapi (fun k _ -> - -. v1.(k) *. expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) - -. v2.(k) *. (expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz)) - ) coef_prod + Array.mapi (fun k _ -> -. v1.(k) *. f1.(k) -. v2.(k) *. f2.(k)) coef_prod in let p2 = if cmm.(xyz) < 0 then p1 else - let v1 = - vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) - and v2 = - vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) - and fcm = + let fcm = (float_of_int cm.(xyz)) *. 0.5 in - Array.mapi (fun k _ -> p1.(k) +. fcm *. expo_inv_q.(k) - *. (v1.(k) +. expo_inv_q.(k) *. v2.(k)) - ) coef_prod + let f1 = + Array.mapi (fun k _ -> fcm *. expo_inv_q.(k) ) coef_prod + and f2 = + Array.mapi (fun k _ -> f1.(k) *. expo_inv_q.(k) ) coef_prod + in + let v1 = + if (at_least_one_valid f1) then + vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) + else empty + in + let v2 = + if (at_least_one_valid f2) then + vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) + else empty + in + Array.mapi (fun k _ -> p1.(k) +. f1.(k) *. v1.(k) +. f2.(k) *. v2.(k)) coef_prod in if (am.(xyz) < 0) || (cm.(xyz) < 0) then p2 else let fa = (float_of_int angMom_a.(xyz)) *. expo_inv_p *. 0.5 in - (* - if (abs_float fa < cutoff) then empty else - *) + let f1 = + Array.mapi (fun k _ -> fa *. expo_inv_q.(k) ) coef_prod + in + if (at_least_one_valid f1) then let v = vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) in - Array.mapi (fun k _ -> - p2.(k) -. fa *. expo_inv_q.(k) *. v.(k) - ) coef_prod + Array.mapi (fun k _ -> p2.(k) -. f1.(k) *. v.(k)) coef_prod + else p2 + ) in if not found then @@ -158,27 +187,27 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) (** Horizontal recurrence relations *) - and hrr0_v m angMom_a angMom_b angMom_c + and hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c = match totAngMom_b with | 0 -> begin match (totAngMom_a, totAngMom_c) with - | (0,0) -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod - | (_,0) -> vrr0_v m angMom_a totAngMom_a - | (_,_) -> vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c + | (0,0) -> Array.mapi (fun k c -> c *. zero_m_array.(k).(0)) coef_prod + | (_,0) -> vrr0_v 0 angMom_a totAngMom_a + | (_,_) -> vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c end | 1 -> let xyz = if angMom_b.(0) = 1 then 0 else if angMom_b.(1) = 1 then 1 else 2 in let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in ap.(xyz) <- ap.(xyz) + 1; let f = Coordinate.coord center_ab xyz in let v1 = - vrr_v m ap angMom_c (totAngMom_a+1) totAngMom_c + vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c in if (abs_float f < cutoff) then v1 else let v2 = - vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c + vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c in Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2 | _ -> @@ -194,20 +223,20 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) bm.(xyz) <- bm.(xyz) - 1; if (bm.(xyz) < 0) then empty else let h1 = - hrr0_v m ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c + hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in let f = (Coordinate.coord center_ab xyz) in if (abs_float f < cutoff) then h1 else let h2 = - hrr0_v m angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c + hrr0_v angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c in Array.map2 (fun h1 h2 -> h1 +. h2 *. f) h1 h2 - and hrr_v m angMom_a angMom_b angMom_c angMom_d + and hrr_v angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d = match (totAngMom_b, totAngMom_d) with - | (0,0) -> vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c - | (_,0) -> hrr0_v m angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c + | (0,0) -> vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c + | (_,0) -> hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] and dm = [| angMom_d.(0) ; angMom_d.(1) ; angMom_d.(2) |] @@ -220,13 +249,13 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) cp.(xyz) <- cp.(xyz) + 1; dm.(xyz) <- dm.(xyz) - 1; let h1 = - hrr_v m angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) + hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) and h2 = - hrr_v m angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) + hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in Array.mapi (fun k center_cd -> h1.(k) +. h2.(k) *. (Coordinate.coord center_cd xyz)) center_cd in - hrr_v m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b + hrr_v angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d @@ -361,7 +390,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q ) shell_q in let integral = - hvrr_two_e_vector 0 (angMomA, angMomB, angMomC, angMomD) + hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD) (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) (maxm, zero_m_array)