diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 8e7638b..eff125f 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -7,7 +7,8 @@ module Zm = struct let name = "Electron repulsion integrals" - let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = + let zero_m ~maxm ~expo_p_inv ~expo_q_inv ~norm_pq_sq = + let expo_pq_inv = expo_p_inv +. expo_q_inv in assert (expo_pq_inv <> 0.); let norm_pq_sq = if norm_pq_sq > integrals_cutoff then norm_pq_sq else 0. diff --git a/Basis/TwoElectronIntegrals.ml b/Basis/TwoElectronIntegrals.ml index 34dcf7b..ade4e5f 100644 --- a/Basis/TwoElectronIntegrals.ml +++ b/Basis/TwoElectronIntegrals.ml @@ -17,7 +17,8 @@ module Fis = FourIdxStorage module type Zero_mType = sig val name : string - val zero_m : maxm:int -> expo_pq_inv:float -> norm_pq_sq:float -> float array + val zero_m : maxm:int -> expo_p_inv:float -> expo_q_inv:float -> + norm_pq_sq:float -> float array end diff --git a/Basis/TwoElectronIntegrals.mli b/Basis/TwoElectronIntegrals.mli index db90cff..fc827b8 100644 --- a/Basis/TwoElectronIntegrals.mli +++ b/Basis/TwoElectronIntegrals.mli @@ -11,7 +11,8 @@ module type Zero_mType = val name : string (** Name of the kind of integrals, for printing purposes. *) -val zero_m : maxm:int -> expo_pq_inv:float -> norm_pq_sq:float -> float array +val zero_m : maxm:int -> expo_p_inv:float -> expo_q_inv:float -> + norm_pq_sq:float -> float array (** The returned float array contains all the {% $(00|00)^m$ %} values, where [m] is the index of the array. diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index ca71c9e..fb626d5 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -21,8 +21,8 @@ type four_idx_intermediates = { expo_b : float ; expo_d : float ; - expo_inv_p : float ; - expo_inv_q : float ; + expo_p_inv : float ; + expo_q_inv : float ; center_ab : Co.t ; center_cd : Co.t ; center_pq : Co.t ; @@ -41,8 +41,8 @@ let rec hvrr_two_e let abcd = { expo_b = abcd.expo_d ; expo_d = abcd.expo_b ; - expo_inv_p = abcd.expo_inv_q ; - expo_inv_q = abcd.expo_inv_p ; + expo_p_inv = abcd.expo_q_inv ; + expo_q_inv = abcd.expo_p_inv ; center_ab = abcd.center_cd ; center_cd = abcd.center_ab ; center_pq = Co.neg abcd.center_pq ; @@ -67,8 +67,8 @@ let rec hvrr_two_e | _ -> Co.Z in - let expo_inv_p = abcd.expo_inv_p - and expo_inv_q = abcd.expo_inv_q + let expo_p_inv = abcd.expo_p_inv + and expo_q_inv = abcd.expo_q_inv and center_ab = abcd.center_ab and center_cd = abcd.center_cd and center_pq = abcd.center_pq @@ -89,8 +89,8 @@ let rec hvrr_two_e let am = Po.decr xyz angMom_a in let amxyz = Po.get xyz am in - let f1 = expo_inv_p *. Co.get xyz center_pq - and f2 = abcd.expo_b *. expo_inv_p *. Co.get xyz center_ab + let f1 = expo_p_inv *. Co.get xyz center_pq + and f2 = abcd.expo_b *. expo_p_inv *. Co.get xyz center_ab in let result = Array.create_float (maxsze - angMom_a.Po.tot) in if amxyz = 0 then @@ -104,10 +104,10 @@ let rec hvrr_two_e let amm = Po.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 + let f3 = (float_of_int amxyz) *. expo_p_inv *. 0.5 in Array.iteri (fun m _ -> result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) - +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) ) result + +. f3 *. (v3.(m) +. expo_p_inv *. v3.(m+1)) ) result end; result in Zmap.add map_1d key result; @@ -132,16 +132,16 @@ let rec hvrr_two_e let axyz = Po.get xyz angMom_a in let f1 = - -. abcd.expo_d *. expo_inv_q *. Co.get xyz center_cd + -. abcd.expo_d *. expo_q_inv *. Co.get xyz center_cd and f2 = - expo_inv_q *. Co.get xyz center_pq + expo_q_inv *. Co.get xyz center_pq in let result = Array.make (maxsze - angMom_a.Po.tot - angMom_c.Po.tot) 0. in if axyz > 0 then begin let am = Po.decr xyz angMom_a in let f5 = - (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 + (float_of_int axyz) *. expo_p_inv *. expo_q_inv *. 0.5 in if (abs_float f5 > cutoff) then let v5 = @@ -153,10 +153,10 @@ let rec hvrr_two_e if cmxyz > 0 then begin let f3 = - (float_of_int cmxyz) *. expo_inv_q *. 0.5 + (float_of_int cmxyz) *. expo_q_inv *. 0.5 in if (abs_float f3 > cutoff) || - (abs_float (f3 *. expo_inv_q) > cutoff) then + (abs_float (f3 *. expo_q_inv) > cutoff) then begin let v3 = let cmm = Po.decr xyz cm in @@ -164,7 +164,7 @@ let rec hvrr_two_e in Array.iteri (fun m _ -> result.(m) <- result.(m) +. - f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) ) result + f3 *. (v3.(m) +. expo_q_inv *. v3.(m+1)) ) result end end; if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then @@ -196,7 +196,7 @@ let rec hvrr_two_e let cm = Po.decr xyz angMom_c in let cmxyz = Po.get xyz cm in - let expo_inv_q_over_p = expo_inv_q /. expo_inv_p in + let expo_inv_q_over_p = expo_q_inv /. expo_p_inv in let f = Co.get xyz center_qc +. expo_inv_q_over_p *. Co.get xyz center_pa @@ -205,7 +205,7 @@ let rec hvrr_two_e let result = if cmxyz < 1 then result else - let f = 0.5 *. (float_of_int cmxyz) *. expo_inv_q in + let f = 0.5 *. (float_of_int cmxyz) *. expo_q_inv in if abs_float f < cutoff then 0. else let cmm = Po.decr xyz cm in let v3 = trr angMom_a cmm in @@ -225,7 +225,7 @@ let rec hvrr_two_e in let result = if axyz < 1 then result else - let f = 0.5 *. (float_of_int axyz) *. expo_inv_q in + let f = 0.5 *. (float_of_int axyz) *. expo_q_inv in if abs_float f < cutoff then result else let am = Po.decr xyz angMom_a in let v2 = trr am cm in @@ -332,16 +332,15 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t and sp_cd = Pspc.shell_pair_q spc in - let expo_inv_p = Psp.exponent_inv sp_ab + let expo_p_inv = Psp.exponent_inv sp_ab in let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in let norm_pq_sq = Co.dot center_pq center_pq in - let expo_inv_q = Psp.exponent_inv sp_cd in - let expo_pq_inv = expo_inv_p +. expo_inv_q in + let expo_q_inv = Psp.exponent_inv sp_cd in let zero_m_array = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + zero_m ~maxm ~expo_p_inv ~expo_q_inv ~norm_pq_sq in begin @@ -382,7 +381,7 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t let coef_prod = coef_prod *. norm in let abcd = { - expo_b ; expo_d ; expo_inv_p ; expo_inv_q ; + expo_b ; expo_d ; expo_p_inv ; expo_q_inv ; center_ab ; center_cd ; center_pq ; center_pa ; center_qc ; zero_m_array ; } in @@ -444,16 +443,15 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : and sp_cd = Pspc.shell_pair_q spc in - let expo_inv_p = Psp.exponent_inv sp_ab + let expo_p_inv = Psp.exponent_inv sp_ab in let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in let norm_pq_sq = Co.dot center_pq center_pq in - let expo_inv_q = Psp.exponent_inv sp_cd in - let expo_pq_inv = expo_inv_p +. expo_inv_q in + let expo_q_inv = Psp.exponent_inv sp_cd in let zero_m_array = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + zero_m ~maxm ~expo_p_inv ~expo_q_inv ~norm_pq_sq in begin @@ -494,7 +492,7 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : let coef_prod = coef_prod *. norm in let abcd = { - expo_b ; expo_d ; expo_inv_p ; expo_inv_q ; + expo_b ; expo_d ; expo_p_inv ; expo_q_inv ; center_ab ; center_cd ; center_pq ; center_pa ; center_qc ; zero_m_array ; } in diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index af043a0..f232cbc 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -23,17 +23,36 @@ let at_least_one_valid arr = Array.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false with Found -> true +type fou_idx_intermediate = +{ + expo_b : float array; + expo_d : float array; + expo_p_inv : float array; + expo_q_inv : float array; + center_ab : Co.t ; + center_cd : Co.t ; + center_pq : Co.axis -> float array array; + center_pa : Co.axis -> float array; + center_qc : Co.axis -> float array; + zero_m_array : float array array array; +} + + (** Horizontal and Vertical Recurrence Relations (HVRR) *) let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) - zero_m_array - (expo_b, expo_d) - (expo_inv_p, expo_inv_q) - (center_ab, center_cd, center_pq) - (center_pa, center_qc) - map_1d map_2d np nq + abcd map_1d map_2d np nq = + let expo_p_inv = abcd.expo_p_inv + and expo_q_inv = abcd.expo_q_inv + and center_ab = abcd.center_ab + and center_cd = abcd.center_cd + and center_pq = abcd.center_pq + in + + let zero_m_array = abcd.zero_m_array in + let maxm = Array.length zero_m_array - 1 in let get_xyz angMom = @@ -62,10 +81,11 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) begin if abs_float cab >= cutoff then + let expo_b = abcd.expo_b in 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 + let f0 = -. expo_b.(l) *. expo_p_inv.(l) *. cab and v0_l = v0.(l) in Array.iteri (fun k v0_lk -> @@ -87,14 +107,14 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk ) p0_l ) result - ) expo_inv_p + ) expo_p_inv else begin let amm = Po.decr xyz am in let amxyz = float_of_int amxyz in 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_p_inv.(l) *. 0.5 and center_pq_xyz_l = (center_pq xyz).(l) in Array.iteri (fun m result_m -> @@ -111,7 +131,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) ) p0.(l) ) result - ) expo_inv_p + ) expo_p_inv end; result in @@ -139,8 +159,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let f = -. (Co.get xyz center_cd) in let f1 = + let expo_d = abcd.expo_d in Array.init nq (fun k -> - let x = expo_d.(k) *. expo_inv_q.(k) *. f in + let x = expo_d.(k) *. expo_q_inv.(k) *. f in if ( (not !do_compute) && (abs_float x > cutoff) ) then do_compute := true; x) @@ -163,7 +184,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) Array.init np (fun l -> let cpq_l = (center_pq xyz).(l) in Array.init nq (fun k -> - let x = expo_inv_q.(k) *. cpq_l.(k) in + let x = expo_q_inv.(k) *. cpq_l.(k) in if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) ) ) @@ -212,7 +233,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let fcm = (float_of_int (cxyz-1)) *. 0.5 in let f1 = Array.init nq (fun k -> - let x = fcm *. expo_inv_q.(k) in + let x = fcm *. expo_q_inv.(k) in if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) ) @@ -240,7 +261,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let v3 = let f2 = Array.init nq (fun k -> - let x = expo_inv_q.(k) *. f1.(k) in + let x = expo_q_inv.(k) *. f1.(k) in if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) ) @@ -334,12 +355,12 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | Some p2 -> p2 in for l=0 to np-1 do - let fa = (float_of_int axyz) *. expo_inv_p.(l) *. 0.5 in + let fa = (float_of_int axyz) *. expo_p_inv.(l) *. 0.5 in let p2_l = p2.(l) and v_l = v.(l) in for k=0 to nq-1 do - p2_l.(k) <- p2_l.(k) -. fa *. expo_inv_q.(k) *. v_l.(k) + p2_l.(k) <- p2_l.(k) -. fa *. expo_q_inv.(k) *. v_l.(k) done done; Some p2 @@ -367,7 +388,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) Array.mapi (fun l expo_inv_p_l -> let expo_p_l = 1./.expo_inv_p_l in Array.mapi (fun k expo_inv_q_k -> - expo_inv_q_k *. expo_p_l) expo_inv_q ) expo_inv_p + expo_inv_q_k *. expo_p_l) expo_q_inv ) expo_p_inv in let result = None in @@ -382,7 +403,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) Some (Array.init np (fun l -> let v3_l = v3.(l) in Array.mapi (fun k v3_lk -> - expo_inv_q.(k) *. f *. v3_lk) v3_l + expo_q_inv.(k) *. f *. v3_lk) v3_l ) ) | Some result, None -> Some result | Some result, Some v3 -> @@ -390,7 +411,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result_l = result.(l) in Array.iteri (fun k v3_lk -> result_l.(k) <- result_l.(k) +. - expo_inv_q.(k) *. f *. v3_lk) v3_l + expo_q_inv.(k) *. f *. v3_lk) v3_l ) v3 ; Some result) end in @@ -460,14 +481,14 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result_l = result.(l) in Array.iteri (fun k v2_lk -> result_l.(k) <- result_l.(k) +. - expo_inv_q.(k) *. f *. v2_lk) v2_l + expo_q_inv.(k) *. f *. v2_lk) v2_l ) v2; Some result) | None, None -> None | None, Some v2 -> Some (Array.init np (fun l -> let v2_l = v2.(l) in Array.mapi (fun k v2_lk -> - expo_inv_q.(k) *. f *. v2_lk) v2_l + expo_q_inv.(k) *. f *. v2_lk) v2_l ) ) end in @@ -599,9 +620,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q contracted_class.(0) <- begin try - let expo_inv_p = + let expo_p_inv = Vec.init np (fun ab -> Psp.exponent_inv sp.(ab-1)) - and expo_inv_q = + and expo_q_inv = Vec.init nq (fun cd -> Psp.exponent_inv sq.(cd-1)) in @@ -615,8 +636,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then raise NullQuartet; - let expo_pq_inv = - expo_inv_p.{i} +. expo_inv_q.{j} + let expo_p_inv, expo_q_inv = + expo_p_inv.{i}, expo_q_inv.{j} in let center_pq = @@ -627,7 +648,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let zero_m_array = - zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq + zero_m ~maxm:0 ~expo_p_inv ~expo_q_inv ~norm_pq_sq in zero_m_array.(0) with NullQuartet -> 0. @@ -644,9 +665,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm = Cspc.norm_scales shell_pair_couple in - let expo_inv_p = + let expo_p_inv = Array.map (fun shell_ab -> Psp.exponent_inv shell_ab) sp - and expo_inv_q = + and expo_q_inv = Array.map (fun shell_cd -> Psp.exponent_inv shell_cd) sq in @@ -729,8 +750,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q if (abs_float coef.(ab).(cd) < cutoff) then empty else - let expo_pq_inv = - expo_inv_p.(ab) +. expo_inv_q.(cd) + let expo_p_inv, expo_q_inv = + expo_p_inv.(ab), expo_q_inv.(cd) in let norm_pq_sq = let x = (center_pq X).(ab).(cd) @@ -740,7 +761,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q x *. x +. y *. y +. z *. z in - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + zero_m ~maxm ~expo_p_inv ~expo_q_inv ~norm_pq_sq ) sq in (* Transpose result *) @@ -799,15 +820,17 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; ); + let abcd = + { expo_b ; expo_d ; expo_p_inv ; expo_q_inv ; + center_ab = Csp.a_minus_b shell_p; + center_cd = Csp.a_minus_b shell_q ; + center_pq ; center_pa ; + center_qc ; zero_m_array } + in + let integral = hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) - zero_m_array - (expo_b, expo_d) - (expo_inv_p, expo_inv_q) - (Csp.a_minus_b shell_p, - Csp.a_minus_b shell_q, center_pq) - (center_pa, center_qc) - map_1d map_2d np nq + abcd map_1d map_2d np nq in contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i) with NullQuartet -> ()