open Qcaml_common open Qcaml_particles open Util open Constants open Qcaml_gaussian_basis exception NullPair module Am = Angular_momentum module Co = Coordinate module Cs = Contracted_shell module Csp = Contracted_shell_pair module Po = Powers module Ps = Primitive_shell module Psp = Primitive_shell_pair (** Horizontal and Vertical Recurrence Relations (HVRR) *) let hvrr_one_e angMom_a angMom_b zero_m_array expo_inv_p center_ab center_pa center_pc map = let maxm = angMom_a.Po.tot + angMom_b.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 vrr angMom_a = let { Po.x=ax ; y=ay ; z=az ; _ } = angMom_a in if ax < 0 || ay < 0 || az < 0 then raise Exit else match angMom_a.Po.tot with | 0 -> zero_m_array | _ -> let key = Zkey.of_powers_three angMom_a in try Zmap.find map 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 = Co.get xyz center_pa in let f2 = expo_inv_p *. Co.get xyz center_pc in if amxyz < 1 then let v1 = vrr am in Array.init (maxsze - angMom_a.Po.tot) (fun m -> f1 *. v1.(m) -. f2 *. v1.(m+1)) else let v3 = let amm = Po.decr xyz am in vrr amm in let v1 = vrr am in let f3 = float_of_int_fast amxyz *. expo_inv_p *. 0.5 in Array.init (maxsze - angMom_a.Po.tot) (fun m -> f1 *. v1.(m) -. f2 *. v1.(m+1) +. f3 *. (v3.(m) -. expo_inv_p *. v3.(m+1)) ) in Zmap.add map key result; result (* Horizontal recurrence relations *) and hrr angMom_a angMom_b = match angMom_b.Po.tot with | 0 -> (vrr angMom_a).(0) | _ -> let xyz = get_xyz angMom_b in let bxyz = Po.get xyz angMom_b in if (bxyz < 1) then 0. else let ap = Po.incr xyz angMom_a in let bm = Po.decr xyz angMom_b in let h1 = hrr ap bm in let f2 = Co.get xyz center_ab in if abs_float f2 < integrals_cutoff then h1 else let h2 = hrr angMom_a bm in h1 +. f2 *. h2 in hrr angMom_a angMom_b (** Computes all the one-electron integrals of the contracted shell pair *) let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let shell_a = Csp.shell_a shell_p and shell_b = Csp.shell_b shell_p in let maxm = Am.to_int (Csp.ang_mom shell_p) in (* Pre-computation of integral class indices *) let class_indices = Csp.zkey_array shell_p in let contracted_class = Array.make (Array.length class_indices) 0.; in (* Compute all integrals in the shell for each pair of significant shell pairs *) let norm_scales_p = Csp.norm_scales shell_p in let center_ab = Csp.a_minus_b shell_p in List.iter (fun (coef_prod, psp) -> try begin (* Screening on the product of coefficients *) if abs_float coef_prod < 1.e-3 *. integrals_cutoff then raise NullPair; let expo_pq_inv = Psp.exponent_inv psp and center_p = Psp.center psp and center_pa = Psp.center_minus_a psp in Array.iter (fun (element, nucl_coord) -> let charge = Element.to_charge element |> Charge.to_float in let center_pc = Co.(center_p |- nucl_coord ) in let norm_pq_sq = Co.dot center_pc center_pc in let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in match Cs.(ang_mom shell_a, ang_mom shell_b) with | Am.(S,S) -> let integral = zero_m_array.(0) in contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge | _ -> let map = Zmap.create (2*maxm) in let norm_scales = norm_scales_p in (* Compute the integral class from the primitive shell quartet *) class_indices |> Array.iteri (fun i key -> let (angMomA,angMomB) = match Zkey.to_powers key with | Zkey.Six x -> x | _ -> assert false in let norm = norm_scales.(i) in let coef_prod = coef_prod *. norm in let integral = hvrr_one_e angMomA angMomB zero_m_array expo_pq_inv center_ab center_pa center_pc map in contracted_class.(i) <- contracted_class.(i) -. coef_prod *. integral *. charge ) ) geometry end with NullPair -> () ) (Csp.coefs_and_shell_pairs shell_p); let result = Zmap.create (Array.length contracted_class) in Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; result