open Util open Constants type t = { expo : float; expo_inv : float; center_ab: Coordinate.t; center_a : Coordinate.t; center : Coordinate.t; norm_sq : float; norm_coef: float; coef : float; norm_coef_scale : float array; i : int; j : int; shell_a : Contracted_shell.t; shell_b : Contracted_shell.t; } exception Null_contribution let create_array ?cutoff p_a p_b = let cutoff, log_cutoff = match cutoff with | None -> -1., max_float | Some cutoff -> cutoff, -. (log cutoff) in let center_ab = Coordinate.( Contracted_shell.center p_a |- Contracted_shell.center p_b ) in let norm_sq = Coordinate.dot center_ab center_ab in let norm_coef_scale_a = Contracted_shell.norm_coef_scale p_a and norm_coef_scale_b = Contracted_shell.norm_coef_scale p_b in let norm_coef_scale = Array.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_b ) norm_coef_scale_a |> Array.to_list |> Array.concat in Array.init (Contracted_shell.size p_a) (fun i -> let p_a_expo_center = Coordinate.( Contracted_shell.expo p_a i |. Contracted_shell.center p_a ) in let norm_coef_a = Contracted_shell.norm_coef p_a i in Array.init (Contracted_shell.size p_b) (fun j -> try let norm_coef_b = Contracted_shell.norm_coef p_b j in let norm_coef = norm_coef_a *. norm_coef_b in if (norm_coef < cutoff) then raise Null_contribution; let p_b_expo_center = Coordinate.( Contracted_shell.expo p_b j |. Contracted_shell.center p_b ) in let expo = Contracted_shell.(expo p_a i +. expo p_b j) in let expo_inv = 1. /. expo in let center = Coordinate.( expo_inv |. (p_a_expo_center |+ p_b_expo_center ) ) in let argexpo = Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv in if (argexpo > log_cutoff) then raise Null_contribution; let g = (pi *. expo_inv)**(1.5) *. exp(-. argexpo) in let coef = norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g in if (abs_float coef < cutoff) then raise Null_contribution; let center_a = Coordinate.(center |- Contracted_shell.center p_a) in Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale } with | Null_contribution -> None ) ) |> Array.to_list |> Array.concat |> Array.to_list |> List.filter (function Some _ -> true | None -> false) |> List.map (function Some x -> x | None -> assert false) |> Array.of_list open Util