open Util open Lacaml.D open Bigarray let cutoff = Constants.cutoff let cutoff2 = cutoff *. cutoff exception NullQuartet exception Found let at_least_one_valid arr = try Array.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false with Found -> true (** Horizontal and Vertical Recurrence Relations (HVRR) *) let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (totAngMom_a, totAngMom_b, totAngMom_c, totAngMom_d) (maxm, zero_m_array) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) map_1d map_2d np nq = let totAngMom_a = Angular_momentum.to_int totAngMom_a and totAngMom_b = Angular_momentum.to_int totAngMom_b and totAngMom_c = Angular_momentum.to_int totAngMom_c and totAngMom_d = Angular_momentum.to_int totAngMom_d in (** Vertical recurrence relations *) let rec vrr0_v m angMom_a = function | 0 -> Some zero_m_array.(m) | totAngMom_a -> let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in try Zmap.find map_1d.(m) 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, 0 | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1 | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2 in if amxyz < 0 then None else begin let cab = Coordinate.coord center_ab xyz in let v1_top, p1_top = if abs_float cab < cutoff then None, vrr0_v (m+1) am (totAngMom_a-1) else vrr0_v m am (totAngMom_a-1), vrr0_v (m+1) am (totAngMom_a-1) in let v1_top2, p1_top2 = if amxyz < 1 then (None,None) else vrr0_v m amm (totAngMom_a-2), vrr0_v (m+1) amm (totAngMom_a-2) in let result = Array.make_matrix np nq 0. in if amxyz < 1 then begin let p0 = match p1_top with | Some p1_top -> p1_top | _ -> assert false in begin match v1_top with | None -> () | Some v0 -> for l=0 to np-1 do let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab in for k=0 to nq-1 do result.(l).(k) <- v0.(l).(k) *. f0 done done end; for l=0 to np-1 do for k=0 to nq-1 do result.(l).(k) <- result.(l).(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) done done end else begin let p0 = match p1_top with | Some p1_top -> p1_top | _ -> assert false in begin match v1_top with | None -> () | Some v0 -> for l=0 to np-1 do let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab in for k=0 to nq-1 do result.(l).(k) <- v0.(l).(k) *. f0 done done end; let v1 = match v1_top2 with | Some v1_top2 -> v1_top2 | None -> assert false in let v2 = match p1_top2 with | Some p1_top2 -> p1_top2 | None -> assert false in for l=0 to np-1 do let f = (float_of_int amxyz) *. expo_inv_p.(l) *. 0.5 in for k=0 to nq-1 do result.(l).(k) <- result.(l).(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) +. f *. (v1.(l).(k) +. v2.(l).(k) *. expo_inv_p.(l)) done done end; Some result end in Zmap.add map_1d.(m) key result; result and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c = match (totAngMom_a, totAngMom_c) with | (i,0) -> vrr0_v m angMom_a totAngMom_a | (_,_) -> let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c)) in try Zmap.find map_2d.(m) key with | Not_found -> let result = begin let am, cm, cmm, axyz, cxyz, xyz = let (aax, aay, aaz) = angMom_a and (acx, acy, acz) = angMom_c in if (acz > 0) then (aax, aay, aaz-1), (acx, acy, acz-1), (acx, acy, acz-2), aaz, acz, 2 else if (acy > 0) then (aax, aay-1,aaz), (acx, acy-1,acz), (acx, acy-2,acz), aay,acy, 1 else (aax-1,aay,aaz), (acx-1,acy,acz), (acx-2,acy,acz), aax,acx, 0 in (* let result = Array.make_matrix np nq 0. in *) let do_compute = ref false in let v1 = let f = -. (Coordinate.coord center_cd xyz) in let f1 = Array.init nq (fun k -> let x = expo_d.(k) *. expo_inv_q.(k) *. f in if ( (not !do_compute) && (abs_float x > cutoff) ) then do_compute := true; x) in if (!do_compute) then match vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) with | None -> None | Some v1 -> begin let result = Array.make_matrix np nq 0. in for l=0 to np-1 do for k=0 to nq-1 do result.(l).(k) <- v1.(l).(k) *. f1.(k) done done; Some result end else None in let v2 = let f2 = Array.init np (fun l -> Array.init nq (fun k -> let x = expo_inv_q.(k) *. center_pq.(xyz).(l).(k) in if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) ) ) in if (!do_compute) then match vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) with | None -> None | Some v2 -> begin for l=0 to np-1 do for k=0 to nq-1 do f2.(l).(k) <- -. v2.(l).(k) *. f2.(l).(k) done done; Some f2 end else None in let p1 = match v1, v2 with | None, None -> None | None, Some v2 -> Some v2 | Some v1, None -> Some v1 | Some v1, Some v2 -> begin for l=0 to np-1 do for k=0 to nq-1 do v2.(l).(k) <- v2.(l).(k) +. v1.(l).(k) done done; Some v2 end in let p2 = if cxyz < 2 then p1 else 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 if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) ) in let v1 = if (!do_compute) then match vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) with | None -> None | Some v1 -> begin let result = Array.make_matrix np nq 0. in for l=0 to np-1 do for k=0 to nq-1 do result.(l).(k) <- v1.(l).(k) *. f1.(k) done; done; Some result end else None in let v3 = let f2 = Array.init nq (fun k -> let x = expo_inv_q.(k) *. f1.(k) in if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) ) in if (!do_compute) then match vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) with | None -> None | Some v3 -> begin let result = Array.make_matrix np nq 0. in for l=0 to np-1 do for k=0 to nq-1 do result.(l).(k) <- v3.(l).(k) *. f2.(k) done done; Some result end else None in match p1, v1, v3 with | None, None, None -> None | Some p1, None, None -> Some p1 | None, Some v1, None -> Some v1 | None, None, Some v3 -> Some v3 | Some p1, Some v1, Some v3 -> begin for l=0 to np-1 do for k=0 to nq-1 do v3.(l).(k) <- p1.(l).(k) +. v1.(l).(k) +. v3.(l).(k) done done; Some v3 end | Some p1, Some v1, None -> begin for l=0 to np-1 do for k=0 to nq-1 do p1.(l).(k) <- v1.(l).(k) +. p1.(l).(k) done done; Some p1 end | Some p1, None, Some v3 -> begin for l=0 to np-1 do for k=0 to nq-1 do p1.(l).(k) <- p1.(l).(k) +. v3.(l).(k) done done; Some p1 end | None , Some v1, Some v3 -> begin for l=0 to np-1 do for k=0 to nq-1 do v3.(l).(k) <- v1.(l).(k) +. v3.(l).(k) done done; Some v3 end in if (axyz < 1) || (cxyz < 1) then p2 else let v = vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) in match (p2, v) with | None, None -> None | Some p2, None -> Some p2 | _, Some v -> begin let p2 = match p2 with | None -> Array.make_matrix np nq 0. | Some p2 -> p2 in for l=0 to np-1 do let fa = (float_of_int axyz) *. expo_inv_p.(l) *. 0.5 in for k=0 to nq-1 do p2.(l).(k) <- p2.(l).(k) -. fa *. expo_inv_q.(k) *. v.(l).(k) done done; Some p2 end end in Zmap.add map_2d.(m) key result; result (** Horizontal recurrence relations *) and hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c = match totAngMom_b with | 0 -> begin match (totAngMom_a, totAngMom_c) with | (0,0) -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. zero_m_array.(0) | (_,_) -> begin match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. end end | 1 -> let (aax, aay, aaz) = angMom_a in let ap, xyz = match angMom_b with | (_,_,1) -> (aax,aay,aaz+1), 2 | (_,1,_) -> (aax,aay+1,aaz), 1 | (_,_,_) -> (aax+1,aay,aaz), 0 in let f = Coordinate.coord center_ab xyz in let v1 = match vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. in if (abs_float f < cutoff) then v1 else let v2 = match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. in v1 +. v2 *. f | _ -> let (aax, aay, aaz) = angMom_a and (abx, aby, abz) = angMom_b in let bxyz, xyz = match angMom_b with | (0,0,_) -> abz, 2 | (0,_,_) -> aby, 1 | _ -> abx, 0 in if (bxyz < 1) then 0. else let ap, bm = match xyz with | 0 -> (aax+1,aay,aaz),(abx-1,aby,abz) | 1 -> (aax,aay+1,aaz),(abx,aby-1,abz) | _ -> (aax,aay,aaz+1),(abx,aby,abz-1) in let h1 = hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in let f = (Coordinate.coord center_ab xyz) in if (abs_float f < cutoff) then h1 else let h2 = hrr0_v angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c in h1 +. h2 *. f and hrr_v 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 begin match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. end else hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> let (acx, acy, acz) = angMom_c and (adx, ady, adz) = angMom_d in let cp, dm, xyz = match angMom_d with | (_,0,0) -> (acx+1, acy, acz), (adx-1, ady, adz), 0 | (_,_,0) -> (acx, acy+1, acz), (adx, ady-1, adz), 1 | _ -> (acx, acy, acz+1), (adx, ady, adz-1), 2 in let h1 = hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) in let f = (Coordinate.coord center_cd xyz) in if (abs_float f < cutoff) then h1 else let h2 = hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in h1 +. f *. h2 in hrr_v (angMom_a.(0),angMom_a.(1),angMom_a.(2)) (angMom_b.(0),angMom_b.(1),angMom_b.(2)) (angMom_c.(0),angMom_c.(1),angMom_c.(2)) (angMom_d.(0),angMom_d.(1),angMom_d.(2)) 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 (** Screening on the product of coefficients *) let coef_max_p = Array.fold_left (fun accu x -> if (abs_float x) > accu then (abs_float x) else accu) 0. shell_p.ContractedShellPair.coef and coef_max_q = Array.fold_left (fun accu x -> if (abs_float x) > accu then (abs_float x) else accu) 0. shell_q.ContractedShellPair.coef in let rec build_list cutoff vec accu = function | -1 -> Array.of_list accu | k -> build_list cutoff vec ( if (abs_float vec.(k) > cutoff) then (k::accu) else accu ) (k-1) in let p_list = let vec = shell_p.ContractedShellPair.coef in build_list (cutoff /. coef_max_q) vec [] (Array.length vec - 1) and q_list = let vec = shell_q.ContractedShellPair.coef in build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1) in (* if (Array.length p_list) < (Array.length shell_p.ContractedShellPair.coef) then begin Printf.printf "Reduced p from %d to %d:\n" (Array.length shell_p.ContractedShellPair.coef) (Array.length p_list); Array.iter (fun k -> Printf.printf "%d " k) p_list; print_newline (); Printf.printf "\n%!" end; if (Array.length q_list) < (Array.length shell_q.ContractedShellPair.coef) then begin Printf.printf "Reduced q from %d to %d:\n" (Array.length shell_q.ContractedShellPair.coef) (Array.length q_list); Array.iter (fun k -> Printf.printf "%d " k) q_list; print_newline (); Printf.printf "\n%!" end; *) let np, nq = Array.length p_list, Array.length q_list in let filter_p vec = Array.init np (fun k -> vec.(p_list.(k))) and filter_q vec = Array.init nq (fun k -> vec.(q_list.(k))) in let sp = filter_p sp and sq = filter_q sq in (* Compute all integrals in the shell for each pair of significant shell pairs *) begin match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d) with | Angular_momentum.(S,S,S,S) -> contracted_class.(0) <- let expo_inv_p = Vec.init np (fun ab -> sp.(ab-1).ShellPair.expo_inv) and expo_inv_q = Vec.init nq (fun cd -> sq.(cd-1).ShellPair.expo_inv) in let coef = let result = Mat.make0 nq np in Lacaml.D.ger (Vec.of_array @@ filter_q shell_q.ContractedShellPair.coef) (Vec.of_array @@ filter_p shell_p.ContractedShellPair.coef) result; result in let zm_array = Mat.init_cols np nq (fun i j -> try 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} in let center_pq = Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center) in let norm_pq_sq = Coordinate.dot center_pq center_pq in let zero_m_array = zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq in zero_m_array.(0) with NullQuartet -> 0. ) in Mat.gemm_trace zm_array coef | _ -> let coef = let cp = filter_p shell_p.ContractedShellPair.coef and cq = filter_q shell_q.ContractedShellPair.coef in Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) in let expo_inv_p = Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp and expo_inv_q = Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq in let expo_b = Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp and expo_d = Array.map (fun shell_cd -> Contracted_shell.expo shell_d shell_cd.ShellPair.j) sq in let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in let center_pq = Array.init 3 (fun xyz -> Array.init np (fun ab -> let shell_ab = sp.(ab) in Array.init nq (fun cd -> let shell_cd = sq.(cd) in let cpq = Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) in match xyz with | 0 -> Coordinate.x cpq; | 1 -> Coordinate.y cpq; | 2 -> Coordinate.z cpq; | _ -> assert false ) ) ) in let zero_m_array = let result = Array.init (maxm+1) (fun _ -> Array.init np (fun _ -> Array.make nq 0. ) ) in let empty = Array.make (maxm+1) 0. in Array.iteri (fun ab shell_ab -> let zero_m_array_tmp = Array.mapi (fun cd shell_cd -> if (abs_float coef.(ab).(cd) < cutoff) then empty else let expo_pq_inv = expo_inv_p.(ab) +. expo_inv_q.(cd) in let norm_pq_sq = center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +. center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +. center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd) in zero_m ~maxm ~expo_pq_inv ~norm_pq_sq ) sq in (* Transpose result *) for m=0 to maxm do for cd=0 to nq-1 do result.(m).(ab).(cd) <- zero_m_array_tmp.(cd).(m) *. coef.(ab).(cd) done done ) sp; result in let norm = let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in Array.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q ) norm_coef_scale_p |> Array.to_list |> Array.concat in let map_1d = Array.init (maxm+1) (fun _ -> Zmap.create (4*maxm)) and map_2d = Array.init (maxm+1) (fun _ -> Zmap.create (Array.length class_indices)) in (* Compute the integral class from the primitive shell quartet *) Array.iteri (fun i key -> let a = Zkey.to_int_array Zkey.Kind_12 key in let (angMomA,angMomB,angMomC,angMomD) = ( [| a.(0) ; a.(1) ; a.(2) |], [| a.(3) ; a.(4) ; a.(5) |], [| a.(6) ; a.(7) ; a.(8) |], [| a.(9) ; a.(10) ; a.(11) |] ) in let integral = hvrr_two_e_vector (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) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (shell_p.ContractedShellPair.center_ab, shell_q.ContractedShellPair.center_ab, center_pq) map_1d map_2d np nq in contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i) ) class_indices end; 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