open Util open Constants module Am = AngularMomentum module Asp = AtomicShellPair module Aspc = AtomicShellPairCouple module Co = Coordinate module Cs = ContractedShell module Csp = ContractedShellPair module Cspc = ContractedShellPairCouple module Po = Powers module Psp = PrimitiveShellPair module Pspc = PrimitiveShellPairCouple module Ps = PrimitiveShell module Zp = Zero_m_parameters let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff exception NullQuartet type four_idx_intermediates = { zero : float array array; gp : float array; gg : float array; gq : float array; coef_g : float array ; center_ra : Co.t array ; center_rc : Co.t array ; center_ab : Co.t ; center_cd : Co.t ; } (** Horizontal and Vertical Recurrence Relations (HVRR) *) let rec hvrr angMom_a angMom_b angMom_c angMom_d abcd map_x map_y map_z = let gp = abcd.gp in let gg = abcd.gg in let gq = abcd.gq in let coef_g = abcd.coef_g in let run xyz map = let zero = match xyz with | Co.X -> abcd.zero.(0) | Co.Y -> abcd.zero.(1) | Co.Z -> abcd.zero.(2) in let angMom_a = Po.get xyz angMom_a in let angMom_b = Po.get xyz angMom_b in let angMom_c = Po.get xyz angMom_c in let angMom_d = Po.get xyz angMom_d in let center_ab = Co.get xyz abcd.center_ab in let center_cd = Co.get xyz abcd.center_cd in let center_ra = Array.map (fun x -> Co.get xyz x) abcd.center_ra in let center_rc = Array.map (fun x -> Co.get xyz x) abcd.center_rc in let rec vrr angMom_a angMom_c = match angMom_a, angMom_c with | 0, -1 | -1, 0 -> assert false | 0, 0 -> zero | 1, 0 -> let v1 = zero in Array.mapi (fun i v1i -> center_ra.(i) *. v1i ) v1 | 0, 1 -> let v1 = zero in Array.mapi (fun i v1i -> center_rc.(i) *. v1i ) v1 | 1, 1 -> let v1 = vrr 1 0 in let v2 = zero in Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. gg.(i) *. v2.(i) ) v1 | 2, 0 -> let v1 = vrr 1 0 in let v2 = zero in Array.mapi (fun i v1i -> center_ra.(i) *. v1i +. gp.(i) *. v2.(i)) v1 | _, 0 -> let v1 = vrr (angMom_a-1) 0 in let a = float_of_int (angMom_a-1) in let v2 = vrr (angMom_a-2) 0 in Array.mapi (fun i v1i -> center_ra.(i) *. v1i +. a *. gp.(i) *. v2.(i)) v1 | _, 1 -> let v1 = vrr angMom_a 0 in let a = float_of_int angMom_a in let v2 = vrr (angMom_a-1) 0 in Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. a *. gg.(i) *. v2.(i) ) v1 | 0, _ -> let v1 = vrr 0 (angMom_c-1) in let b = float_of_int (angMom_c-1) in let v3 = vrr 0 (angMom_c-2) in Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. b *. gq.(i) *. v3.(i)) v1 | _ -> let v1 = vrr angMom_a (angMom_c-1) in let a = float_of_int angMom_a in let b = float_of_int (angMom_c-1) in let v2 = vrr (angMom_a-1) (angMom_c-1) in let v3 = vrr angMom_a (angMom_c-2) in Array.mapi (fun i v1i -> center_rc.(i) *. v1i +. a *. gg.(i) *. v2.(i) +. b *. gq.(i) *. v3.(i)) v1 in let rec hrr angMom_a angMom_b angMom_c angMom_d = let key = Zkey.of_int_four angMom_a angMom_b angMom_c angMom_d in try Zmap.find map key with | Not_found -> let result = match angMom_b, angMom_d with | -1, 0 | 0, -1 -> assert false | 0, 0 -> vrr angMom_a angMom_c | _, 0 -> let h1 = hrr (angMom_a+1) (angMom_b-1) angMom_c angMom_d in if center_ab = 0. then h1 else let h2 = hrr angMom_a (angMom_b-1) angMom_c angMom_d in Array.mapi (fun i h1i -> h1i +. center_ab *. h2.(i)) h1 | _ -> let h1 = hrr angMom_a angMom_b (angMom_c+1) (angMom_d-1) in if center_cd = 0. then h1 else let h2 = hrr angMom_a angMom_b angMom_c (angMom_d-1) in Array.mapi (fun i h1i -> h1i +. center_cd *. h2.(i)) h1 in (Zmap.add map key result; result) in hrr angMom_a angMom_b angMom_c angMom_d in let x, y, z = (run Co.X map_x), (run Co.Y map_y), (run Co.Z map_z) in let rec aux accu = function | 0 -> accu +. coef_g.(0) *. x.(0) *. y.(0) *. z.(0) | i -> aux (accu +. coef_g.(i) *. x.(i) *. y.(i) *. z.(i)) (i-1) in aux 0. (Array.length x - 1) let contracted_class_shell_pair_couple expo_g_inv coef_g shell_pair_couple : float Zmap.t = (* Pre-computation of integral class indices *) let class_indices = Cspc.zkey_array shell_pair_couple 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 shell_p = Cspc.shell_pair_p shell_pair_couple and shell_q = Cspc.shell_pair_q shell_pair_couple in let center_ab = Csp.a_minus_b shell_p and center_cd = Csp.a_minus_b shell_q in let norm_scales = Cspc.norm_scales shell_pair_couple in List.iter (fun (coef_prod, spc) -> let sp_ab = Pspc.shell_pair_p spc and sp_cd = Pspc.shell_pair_q spc in let expo_p_inv = Psp.exponent_inv sp_ab in let expo_q_inv = Psp.exponent_inv sp_cd in let expo_pgq = Array.map (fun e -> 1. /. (expo_p_inv +. expo_q_inv +. e) ) expo_g_inv in let fp = Array.map (fun e -> expo_p_inv *. e) expo_pgq in let fq = Array.map (fun e -> expo_q_inv *. e) expo_pgq in let gp = let x = 0.5 *. expo_p_inv in Array.map (fun f -> (1. -. f) *. x) fp in let gq = let x = 0.5 *. expo_q_inv in Array.map (fun f -> (1. -. f) *. x) fq in let gg = let x = 0.5 *. expo_q_inv in Array.map (fun f -> f *. x) fp in let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in let center_qc = Psp.center_minus_a sp_cd in let center_pa = Psp.center_minus_a sp_ab in let center_ra = Array.map (fun f -> Co.(center_pa |- (f |. center_pq))) fp in let center_rc = Array.map (fun f -> Co.(center_qc |+ (f |. center_pq))) fq in (* zero_m is defined here *) let zero = Array.map (fun xyz -> let x = Co.get xyz center_pq in Array.mapi (fun i ei -> let fg = expo_g_inv.(i) *. ei in (sqrt fg) *. exp (-. x *. x *. ei )) expo_pgq ) Co.[| X ; Y ; Z |] in begin match Cspc.ang_mom shell_pair_couple with (* | Am.S -> let integral = zero.(0) *. zero.(1) *. zero.(2) in contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral *) | _ -> let map_x, map_y, map_z = Zmap.create (Array.length class_indices), Zmap.create (Array.length class_indices), Zmap.create (Array.length class_indices) 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 let norm = norm_scales.(i) in let coef_prod = coef_prod *. norm in let abcd = { zero ; gp ; gg ; gq ; center_ab ; center_cd ; center_ra ; center_rc ; coef_g ; } in let integral = hvrr angMom_a angMom_b angMom_c angMom_d abcd map_x map_y map_z in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral ) end ) (Cspc.coefs_and_shell_pair_couples shell_pair_couple); let result = Zmap.create (Array.length contracted_class) in Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; result