open Util open Constants open Coordinate let debug=false 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) (totAngMom_a_in, totAngMom_b_in, totAngMom_c_in, totAngMom_d_in) (maxm, zero_m_array) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) map_1d map_2d map_1d' map_2d' = let maxsze = maxm+1 in let totAngMom_a = Angular_momentum.to_int totAngMom_a_in and totAngMom_b = Angular_momentum.to_int totAngMom_b_in and totAngMom_c = Angular_momentum.to_int totAngMom_c_in and totAngMom_d = Angular_momentum.to_int totAngMom_d_in in (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *) if (totAngMom_a+totAngMom_b < totAngMom_c+totAngMom_d) then hvrr_two_e (angMom_c, angMom_d, angMom_a, angMom_b) (totAngMom_c_in, totAngMom_d_in, totAngMom_a_in, totAngMom_b_in) (maxm, zero_m_array) (expo_d, expo_b) (expo_inv_q, expo_inv_p) (center_cd, center_ab, (Coordinate.neg center_pq) ) map_1d' map_2d' map_1d map_2d else let maxm = totAngMom_a + totAngMom_b + totAngMom_c + totAngMom_d in let empty = Array.make (maxm+1) 0. in if debug then begin Printf.printf "\n---- %d %d %d %d ----\n" totAngMom_a totAngMom_b totAngMom_c totAngMom_d; let (x,y,z) = angMom_a in Printf.printf "%d %d %d\n" x y z; let (x,y,z) = angMom_b in Printf.printf "%d %d %d\n" x y z; let (x,y,z) = angMom_c in Printf.printf "%d %d %d\n" x y z; let (x,y,z) = angMom_d in Printf.printf "%d %d %d\n" x y z; Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d expo_inv_p expo_inv_q (get X center_ab) (get Y center_ab) (get Z center_ab) (get X center_cd) (get Y center_cd) (get Z center_cd) (get X center_pq) (get Y center_pq) (get Z center_pq) end; (** Vertical recurrence relations *) let rec vrr0 angMom_a totAngMom_a = if debug then begin let (x,y,z) = angMom_a in Printf.printf "vrr0: %d : %d %d %d\n" totAngMom_a x y z end; match totAngMom_a with | 0 -> zero_m_array | _ -> let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in try Zmap.find map_1d key with | Not_found -> let result = let am, amm, amxyz, xyz = match angMom_a with | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z in if amxyz < 0 then empty else 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 (totAngMom_a-1) 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 v3 = vrr0 amm (totAngMom_a-2) in let v1 = vrr0 am (totAngMom_a-1) 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 and vrr angMom_a angMom_c totAngMom_a totAngMom_c = if debug then begin let angMom_ax, angMom_ay, angMom_az = angMom_a in let angMom_cx, angMom_cy, angMom_cz = angMom_c in Printf.printf "vrr : %d %d : %d %d %d %d %d %d\n" totAngMom_a totAngMom_c angMom_ax angMom_ay angMom_az angMom_cx angMom_cy angMom_cz end; match (totAngMom_a, totAngMom_c) with | (i,0) -> if (i>0) then vrr0 angMom_a totAngMom_a (* OneElectronRR.hvrr_one_e (angMom_a, angMom_b) (totAngMom_a_in, totAngMom_b_in) (maxm, zero_m_array) (expo_b) (expo_inv_p) (center_ab, center_pq, center_ab) map_1d *) else zero_m_array | (_,_) -> let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c) ) in try Zmap.find map_2d key with | Not_found -> let result = let am, cm, cmm, axyz, cmxyz, xyz = let angMom_ax, angMom_ay, angMom_az = angMom_a and angMom_cx, angMom_cy, angMom_cz = angMom_c in match angMom_c with | (_,0,0) -> (angMom_ax-1, angMom_ay, angMom_az), (angMom_cx-1, angMom_cy, angMom_cz), (angMom_cx-2, angMom_cy, angMom_cz), angMom_ax,angMom_cx-1, X | (_,_,0) -> (angMom_ax, angMom_ay-1, angMom_az), (angMom_cx, angMom_cy-1, angMom_cz), (angMom_cx, angMom_cy-2, angMom_cz), angMom_ay,angMom_cy-1, Y | _ -> (angMom_ax, angMom_ay, angMom_az-1), (angMom_cx, angMom_cy, angMom_cz-1), (angMom_cx, angMom_cy, angMom_cz-2), angMom_az,angMom_cz-1, Z in if cmxyz < 0 then empty else 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 totAngMom_a (totAngMom_c-1) 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 v3 = vrr angMom_a cmm totAngMom_a (totAngMom_c-2) 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 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 (totAngMom_a-1) (totAngMom_c-1) 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 (** Horizontal recurrence relations *) and hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c = (* if debug then begin let angMom_ax, angMom_ay, angMom_az = angMom_a and angMom_bx, angMom_by, angMom_bz = angMom_b and angMom_cx, angMom_cy, angMom_cz = angMom_c in Printf.printf "hrr0: %d %d %d : %d %d %d %d %d %d %d %d %d\n" totAngMom_a totAngMom_b totAngMom_c angMom_ax angMom_ay angMom_az angMom_bx angMom_by angMom_bz angMom_cx angMom_cy angMom_cz end; *) match totAngMom_b with | 0 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) | 1 -> let angMom_ax, angMom_ay, angMom_az = angMom_a in let ap, xyz = match angMom_b with | (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), X | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), Y | _ -> (angMom_ax,angMom_ay,angMom_az+1), Z in let v1 = vrr ap angMom_c (totAngMom_a+1) totAngMom_c in let f2 = (Coordinate.get xyz center_ab) in if (abs_float f2 < cutoff) then v1.(0) else let v2 = vrr angMom_a angMom_c totAngMom_a totAngMom_c in v1.(0) +. f2 *. v2.(0) | _ -> let angMom_ax, angMom_ay, angMom_az = angMom_a and angMom_bx, angMom_by, angMom_bz = angMom_b in let bxyz, xyz = match angMom_b with | (_,0,0) -> angMom_bx, X | (_,_,0) -> angMom_by, Y | (_,_,_) -> angMom_bz, Z in if (bxyz < 1) then 0. else let ap, bm = match xyz with | X -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) | Y -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) | Z -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) in let h1 = hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in let f2 = (Coordinate.get xyz center_ab) in if (abs_float f2 < cutoff) then h1 else let h2 = hrr0 angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c in h1 +. f2 *. h2 and hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d = match (totAngMom_b, totAngMom_d) with | (_,0) -> if (totAngMom_b = 0) then (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) else hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> let (angMom_cx, angMom_cy, angMom_cz) = angMom_c and (angMom_dx, angMom_dy, angMom_dz) = angMom_d in let cp, dm, xyz = match angMom_d with | (_,0,0) -> (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), X | (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), Y | _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), Z in let h1 = hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) in let f2 = Coordinate.get xyz center_cd in if (abs_float f2 < cutoff) then h1 else let h2 = hrr angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in h1 +. f2 *. h2 in hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = let shell_a = shell_p.ContractedShellPair.shell_a and shell_b = shell_p.ContractedShellPair.shell_b and shell_c = shell_q.ContractedShellPair.shell_a and shell_d = shell_q.ContractedShellPair.shell_b and sp = shell_p.ContractedShellPair.shell_pairs and sq = shell_q.ContractedShellPair.shell_pairs in let maxm = shell_p.ContractedShellPair.totAngMomInt + shell_q.ContractedShellPair.totAngMomInt in (* Pre-computation of integral class indices *) let class_indices = Angular_momentum.zkey_array (Angular_momentum.Quartet Contracted_shell.(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 = shell_p.ContractedShellPair.monocentric && shell_q.ContractedShellPair.monocentric && Contracted_shell.center shell_p.ContractedShellPair.shell_a = Contracted_shell.center shell_q.ContractedShellPair.shell_a in (* Compute all integrals in the shell for each pair of significant shell pairs *) let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in for ab=0 to (Array.length sp - 1) do let cab = shell_p.ContractedShellPair.coef.(ab) in let b = sp.(ab).ShellPair.j in for cd=0 to (Array.length shell_q.ContractedShellPair.shell_pairs - 1) do let coef_prod = cab *. shell_q.ContractedShellPair.coef.(cd) in (** Screening on the product of coefficients *) try if (abs_float coef_prod) < 1.e-3*.cutoff then raise NullQuartet; let center_pq = sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center in let norm_pq_sq = Coordinate.dot center_pq center_pq in let expo_pq_inv = shell_p.ContractedShellPair.expo_inv.(ab) +. shell_q.ContractedShellPair.expo_inv.(cd) in let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in begin match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d) with | Angular_momentum.(S,S,S,S) -> let integral = zero_m_array.(0) in contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral | _ -> let d = shell_q.ContractedShellPair.shell_pairs.(cd).ShellPair.j in let map_1d = Zmap.create (4*maxm) in let map_2d = Zmap.create (Array.length class_indices) 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 (angMomA,angMomB,angMomC,angMomD) = match Zkey.to_int_tuple ~kind:Zkey.Kind_12 key with | Zkey.Twelve x -> x | _ -> assert false in try if monocentric then begin let ax,ay,az = angMomA and bx,by,bz = angMomB and cx,cy,cz = angMomC and dx,dy,dz = angMomD in if ( ((1 land ax+bx+cx+dx)=1) || ((1 land ay+by+cy+dy)=1) || ((1 land az+bz+cz+dz)=1) ) then raise NullQuartet end; (* (* Schwartz screening *) if (maxm > 2) then ( let schwartz_p = let key = Zkey.of_int_tuple (Zkey.Twelve (angMomA, angMomB, angMomA, angMomB) ) 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_int_tuple (Zkey.Twelve (angMomC, angMomD, angMomC, angMomD) ) 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 (angMomA, angMomB, angMomC, angMomD) (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) (maxm, zero_m_array) (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d) (shell_p.ContractedShellPair.expo_inv.(ab), shell_q.ContractedShellPair.expo_inv.(cd) ) (sp.(ab).ShellPair.center_ab, sq.(cd).ShellPair.center_ab, center_pq) map_1d map_2d 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 = ContractedShellPair.create ~cutoff shell_a shell_b and shell_q = ContractedShellPair.create ~cutoff shell_c shell_d in contracted_class_shell_pairs ~zero_m shell_p shell_q