diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index ea5f70e..46d5641 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -32,7 +32,6 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in let maxsze = maxm+1 in - let empty = Array.make (maxm+1) 0. in let get_xyz angMom = match angMom with @@ -57,36 +56,33 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let am = Powers.decr xyz angMom_a in let amxyz = Powers.get xyz am in - if amxyz >= 0 then - 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 - if amxyz < 1 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) - end - else - begin - let amm = Powers.decr xyz am in - 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 - result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) - +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) - done; - result.(maxm) <- f3 *. v3.(maxm) - end; - result + 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 + 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) + end else - empty + begin + let amm = Powers.decr xyz am in + 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 + result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) + +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) + done; + result.(maxm) <- f3 *. v3.(maxm) + end; + result in Zmap.add map_1d key result; result @@ -102,68 +98,84 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) try Zmap.find map_2d key with | Not_found -> let result = + (* Invariant : angMom_c.tot > 0 so cm.tot >= 0 *) let xyz = get_xyz angMom_c in let cm = Powers.decr xyz angMom_c in let cmxyz = Powers.get xyz cm in let axyz = Powers.get xyz angMom_a in - if cmxyz >= 0 then - let f1 = - -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) - and f2 = - expo_inv_q *. (Coordinate.get xyz center_pq) - in - let result = Array.make maxsze 0. in - if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then + let f1 = + -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) + and f2 = + expo_inv_q *. (Coordinate.get xyz center_pq) + in + let result = Array.make maxsze 0. in + if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then + begin + let v1 = + vrr angMom_a cm + in + for m=0 to maxm-1 do + result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ; + done; + result.(maxm) <- f1 *. v1.(maxm) ; + end; + if cmxyz > 0 then + begin + let f3 = + (float_of_int cmxyz) *. expo_inv_q *. 0.5 + in + if (abs_float f3 > cutoff) || + (abs_float (f3 *. expo_inv_q) > cutoff) then begin - let v1 = - vrr angMom_a cm + let v3 = + let cmm = Powers.decr xyz cm in + vrr angMom_a cmm in for m=0 to maxm-1 do - result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ; + result.(m) <- result.(m) +. + f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) done; - result.(maxm) <- f1 *. v1.(maxm) ; - end; - if cmxyz > 0 then - begin - let f3 = - (float_of_int cmxyz) *. expo_inv_q *. 0.5 - in - if (abs_float f3 > cutoff) || - (abs_float (f3 *. expo_inv_q) > cutoff) then - begin - let v3 = - let cmm = Powers.decr xyz cm in - vrr angMom_a cmm - in - for m=0 to maxm-1 do - result.(m) <- result.(m) +. - f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) - done; - result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm) - end - end; - if axyz > 0 && cmxyz >= 0 then - begin - let am = Powers.decr xyz angMom_a in - let f5 = - (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 - in - if (abs_float f5 > cutoff) then - let v5 = - vrr am cm - in - for m=0 to maxm-1 do - result.(m) <- result.(m) -. f5 *. v5.(m+1) - done - end; - result - else - empty + result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm) + end + end; + if axyz > 0 then + begin + let am = Powers.decr xyz angMom_a in + let f5 = + (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 + in + if (abs_float f5 > cutoff) then + let v5 = + vrr am cm + in + for m=0 to maxm-1 do + result.(m) <- result.(m) -. f5 *. v5.(m+1) + done + end; + result in Zmap.add map_2d key result; result + (* + and trr angMom_a angMom_c = + + match (angMom_a.tot, angMom_c.tot) with + | (i,0) -> if (i>0) then vrr0 angMom_a + else zero_m_array + | (_,_) -> + let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in + + try Zmap.find map_2d key with + | Not_found -> + let result = + let xyz = get_xyz angMom_c in + let cm = Powers.decr xyz angMom_c in + let cmxyz = Powers.get xyz cm in + let axyz = Powers.get xyz angMom_a in + *) + diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index e6e154c..117ac8c 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -35,7 +35,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (** Vertical recurrence relations *) let rec vrr0_v m angMom_a = match angMom_a.tot with - | 0 -> Some zero_m_array.(m) + | 0 -> zero_m_array.(m) | _ -> let key = Zkey.of_powers (Zkey.Three angMom_a) in @@ -45,83 +45,54 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = let xyz = get_xyz angMom_a in let am = Powers.decr xyz angMom_a in - let amxyz = Powers.get xyz am in - if amxyz >= 0 then - begin - let cab = Coordinate.get xyz center_ab in - let v1_top, p1_top = - if abs_float cab < cutoff then - None, - vrr0_v (m+1) am - else - vrr0_v m am, vrr0_v (m+1) am - in - let v1_top2, p1_top2 = - if amxyz < 1 then None, None else - let amm = Powers.decr xyz am in - vrr0_v m amm, vrr0_v (m+1) amm - in + let cab = Coordinate.get xyz center_ab in + let p0 = vrr0_v (m+1) am in - let result = Array.make_matrix np nq 0. in - let p0 = - match p1_top with - | Some p1_top -> p1_top - | _ -> assert false - in - begin - match v1_top with - | None -> () - | Some v0 -> - 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 - end; - let amxyz = Powers.get xyz am in - if amxyz < 1 then + 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 expo_inv_p_l = expo_inv_p.(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 + 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 + 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 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) - 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 - else - begin - let v1 = - match v1_top2 with - | Some v1_top2 -> v1_top2 - | None -> assert false - in - let v2 = - match p1_top2 with - | Some p1_top2 -> p1_top2 - | None -> assert false - in - Array.iteri (fun l result_l -> - let f = float_of_int 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 - end; - Some result - end - else - None + 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 + end; + result in Zmap.add map_1d.(m) key result; result @@ -129,7 +100,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) and vrr_v m angMom_a angMom_c = match (angMom_a.tot, angMom_c.tot) with - | (i,0) -> vrr0_v m angMom_a + | (i,0) -> Some (vrr0_v m angMom_a) | (_,_) -> let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c))