open Util module Am = AngularMomentum module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair module Sp = ShellPair module Po = Powers let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff exception NullQuartet (** Horizontal and Vertical Recurrence Relations (HVRR) *) let rec hvrr_two_e 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 = (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *) if angMom_a.Po.tot + angMom_b.Po.tot < angMom_c.Po.tot + angMom_d.Po.tot then hvrr_two_e angMom_c angMom_d angMom_a angMom_b zero_m_array expo_d expo_b expo_inv_q expo_inv_p center_cd center_ab (Co.neg center_pq) center_qc center_pa map_1d map_2d else let maxm = angMom_a.Po.tot + angMom_b.Po.tot + angMom_c.Po.tot + angMom_d.Po.tot in let maxsze = maxm+1 in let get_xyz angMom = match angMom with | { Po.y=0 ; z=0 ; _ } -> Co.X | { z=0 ; _ } -> Co.Y | _ -> Co.Z in (** Vertical recurrence relations *) let rec vrr0 angMom_a = match angMom_a.Po.tot with | 0 -> zero_m_array | _ -> let key = Zkey.of_powers_three angMom_a in try Zmap.find map_1d key with | Not_found -> let result = let xyz = get_xyz angMom_a in 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 = expo_b *. expo_inv_p *. Co.get xyz center_ab in let result = Array.create_float (maxsze - angMom_a.Po.tot) in if amxyz = 0 then begin let v1 = vrr0 am in Array.iteri (fun m _ -> result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)) result end else begin 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 Array.iteri (fun m _ -> result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) ) result end; result in Zmap.add map_1d key result; result and vrr angMom_a angMom_c = match angMom_a.Po.tot, angMom_c.Po.tot with | (i,0) -> if (i>0) then vrr0 angMom_a else zero_m_array | (_,_) -> let key = Zkey.of_powers_six angMom_a angMom_c in try Zmap.find map_2d key with | Not_found -> let result = (* angMom_c.Po.tot > 0 so cm.Po.tot >= 0 *) let xyz = get_xyz angMom_c in let cm = Po.decr xyz angMom_c in let cmxyz = Po.get xyz cm in let axyz = Po.get xyz angMom_a in let f1 = -. expo_d *. expo_inv_q *. Co.get xyz center_cd and f2 = expo_inv_q *. 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 in if (abs_float f5 > cutoff) then let v5 = vrr am cm in Array.iteri (fun m _ -> result.(m) <- result.(m) -. f5 *. v5.(m+1)) result 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 = Po.decr xyz cm in vrr angMom_a cmm in Array.iteri (fun m _ -> result.(m) <- result.(m) +. f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) ) result end end; if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then begin let v1 = vrr angMom_a cm in 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; result (* and trr angMom_a angMom_c = match (angMom_a.Po.tot, angMom_c.Po.tot) with | (i,0) -> if (i>0) then (vrr0 angMom_a).(0) else zero_m_array.(0) | (_,_) -> let key = Zkey.of_powers_six angMom_a angMom_c in try (Zmap.find map_2d key).(0) with | Not_found -> let result = let xyz = get_xyz angMom_c in let axyz = Po.get xyz angMom_a in 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 f = Co.get xyz center_qc +. expo_inv_q_over_p *. Co.get xyz center_pa in let result = 0. in let result = if cmxyz < 1 then result else let f = 0.5 *. (float_of_int cmxyz) *. expo_inv_q in if abs_float f < cutoff then 0. else let cmm = Po.decr xyz cm in let v3 = trr angMom_a cmm in result +. f *. v3 in let result = if abs_float f < cutoff then result else let v1 = trr angMom_a cm in result +. f *. v1 in let result = if cmxyz < 0 then result else let f = -. expo_inv_q_over_p in let ap = Po.incr xyz angMom_a in let v4 = trr ap cm in result +. v4 *. f in let result = if axyz < 1 then result else let f = 0.5 *. (float_of_int axyz) *. expo_inv_q in if abs_float f < cutoff then result else let am = Po.decr xyz angMom_a in let v2 = trr am cm in result +. f *. v2 in result in Zmap.add map_2d key [|result|]; result *) in let vrr a c = (vrr a c).(0) (* if maxm < 10 then (vrr a c).(0) else trr a c *) in (** Horizontal recurrence relations *) let rec hrr0 angMom_a angMom_b angMom_c = match angMom_b.Po.tot with | 1 -> let xyz = get_xyz angMom_b in let ap = Po.incr xyz angMom_a in let v1 = vrr ap angMom_c in let f2 = Co.get xyz center_ab in if (abs_float f2 < cutoff) then v1 else let v2 = vrr angMom_a angMom_c in v1 +. f2 *. v2 | 0 -> vrr angMom_a angMom_c | _ -> let xyz = get_xyz angMom_b in let bxyz = Po.get xyz angMom_b in if bxyz > 0 then let ap = Po.incr xyz angMom_a in let bm = Po.decr xyz angMom_b in let h1 = hrr0 ap bm angMom_c in let f2 = Co.get xyz center_ab in if abs_float f2 < cutoff then h1 else let h2 = hrr0 angMom_a bm angMom_c in h1 +. f2 *. h2 else 0. and hrr angMom_a angMom_b angMom_c angMom_d = match (angMom_b.Po.tot, angMom_d.Po.tot) with | (_,0) -> if (angMom_b.Po.tot = 0) then vrr angMom_a angMom_c else hrr0 angMom_a angMom_b angMom_c | (_,_) -> let xyz = get_xyz angMom_d in let cp = Po.incr xyz angMom_c in let dm = Po.decr xyz angMom_d in let h1 = hrr angMom_a angMom_b cp dm in let f2 = Co.get xyz center_cd in if abs_float f2 < cutoff then h1 else let h2 = hrr angMom_a angMom_b angMom_c dm in h1 +. f2 *. h2 in hrr angMom_a angMom_b angMom_c angMom_d let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = let shell_a = Csp.shell_a shell_p and shell_b = Csp.shell_b shell_p and shell_c = Csp.shell_a shell_q and shell_d = Csp.shell_b shell_q and sp = Csp.shell_pairs shell_p and sq = Csp.shell_pairs shell_q in let maxm = Csp.totAngMomInt shell_p + Csp.totAngMomInt shell_q in (* Pre-computation of integral class indices *) let class_indices = Am.zkey_array (Am.Quartet Cs.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d )) in let contracted_class = Array.make (Array.length class_indices) 0.; in let monocentric = Csp.monocentric shell_p && Csp.monocentric shell_q && Cs.center (Csp.shell_a shell_p) = Cs.center (Csp.shell_a shell_q) in (* Compute all integrals in the shell for each pair of significant shell pairs *) let norm_coef_scale_p = Csp.norm_coef_scale shell_p in let norm_coef_scale_q = Csp.norm_coef_scale shell_q in for ab=0 to (Array.length sp - 1) do let cab = (Csp.coef shell_p).(ab) in let b = sp.(ab).Sp.j in for cd=0 to (Array.length (Csp.shell_pairs shell_q) - 1) do let coef_prod = cab *. (Csp.coef shell_q).(cd) in (** Screening on the product of coefficients *) try if (abs_float coef_prod) < 1.e-3 *. cutoff then raise NullQuartet; let center_pq = Co.(sp.(ab).Sp.center |- sq.(cd).Sp.center) in let norm_pq_sq = Co.dot center_pq center_pq in let expo_pq_inv = (Csp.expo_inv shell_p).(ab) +. (Csp.expo_inv shell_q).(cd) in let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in begin match Cs.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d) with | Am.(S,S,S,S) -> let integral = zero_m_array.(0) in contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral | _ -> let d = (Csp.shell_pairs shell_q).(cd).Sp.j in let map_1d = Zmap.create (4*maxm) in let map_2d = Zmap.create (Array.length class_indices) in let norm_coef_scale = Array.to_list norm_coef_scale_p |> List.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) |> Array.concat in (* Compute the integral class from the primitive shell quartet *) class_indices |> Array.iteri (fun i key -> let (angMom_a,angMom_b,angMom_c,angMom_d) = match Zkey.to_powers key with | Zkey.Twelve x -> x | _ -> assert false in try if monocentric then begin if ( ((1 land angMom_a.Po.x + angMom_b.Po.x + angMom_c.Po.x + angMom_d.Po.x)=1) || ((1 land angMom_a.Po.y + angMom_b.Po.y + angMom_c.Po.y + angMom_d.Po.y)=1) || ((1 land angMom_a.Po.z + angMom_b.Po.z + angMom_c.Po.z + angMom_d.Po.z)=1) ) then raise NullQuartet end; (* Schwartz screening *) (* if (maxm > 8) then ( let schwartz_p = let key = Zkey.of_powers_twelve angMom_a angMom_b angMom_a angMom_b in match schwartz_p with | None -> 1. | Some schwartz_p -> Zmap.find schwartz_p key in if schwartz_p < cutoff then raise NullQuartet; let schwartz_q = let key = Zkey.of_powers_twelve angMom_c angMom_d angMom_c angMom_d in match schwartz_q with | None -> 1. | Some schwartz_q -> Zmap.find schwartz_q key in if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; ); *) let norm = norm_coef_scale.(i) in let coef_prod = coef_prod *. norm in let integral = hvrr_two_e angMom_a angMom_b angMom_c angMom_d zero_m_array (Cs.expo shell_b).(b) (Cs.expo shell_d).(d) (Csp.expo_inv shell_p).(ab) (Csp.expo_inv shell_q).(cd) sp.(ab).Sp.center_ab sq.(cd).Sp.center_ab center_pq sp.(ab).Sp.center_a sq.(cd).Sp.center_a map_1d map_2d in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral with NullQuartet -> () ) end with NullQuartet -> () done done; let result = Zmap.create (Array.length contracted_class) in Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; result (** Computes all the two-electron integrals of the contracted shell quartet *) let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = let shell_p = Csp.create ~cutoff shell_a shell_b and shell_q = Csp.create ~cutoff shell_c shell_d in contracted_class_shell_pairs ~zero_m shell_p shell_q