open Util open Constants open Coordinate type t = { expo : float; norm_coef : float; norm_coef_scale : float array lazy_t; center : Coordinate.t; totAngMom : AngularMomentum.t; } module Am = AngularMomentum let compute_norm_coef alpha totAngMom = let atot = Am.to_int totAngMom in let factor int_array = let dfa = Array.map (fun j -> ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) ) int_array in sqrt (dfa.(0) *.dfa.(1) *. dfa.(2)) in let expo = if atot mod 2 = 0 then let alpha_2 = alpha +. alpha in (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2)) else let alpha_2 = alpha +. alpha in (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot) in let f a = expo *. (factor a) in f let make totAngMom center expo = let norm_coef_func = compute_norm_coef expo totAngMom in let norm = 1. /. norm_coef_func [| Am.to_int totAngMom ; 0 ; 0 |] in let powers = Am.zkey_array (Am.Singlet totAngMom) in let norm_coef_scale = lazy ( Array.map (fun a -> (norm_coef_func (Zkey.to_int_array a)) *. norm ) powers ) in let norm_coef = 1. /. norm in { expo ; norm_coef ; norm_coef_scale ; center ; totAngMom } let to_string s = let coord = s.center in Printf.sprintf "%1s %8.3f %8.3f %8.3f %16.8e" (Am.to_string s.totAngMom) (get X coord) (get Y coord) (get Z coord) s.expo (** Normalization coefficient of contracted function i, which depends on the exponent and the angular momentum. Two conventions can be chosen : a single normalisation factor for all functions of the class, or a coefficient which depends on the powers of x,y and z. Returns, for each contracted function, an array of functions taking as argument the [|x;y;z|] powers. *) let expo x = x.expo let center x = x.center let totAngMom x = x.totAngMom let norm x = 1. /. x.norm_coef let norm_coef x = x.norm_coef let norm_coef_scale x = Lazy.force x.norm_coef_scale let size_of_shell x = Array.length (norm_coef_scale x)