From 6a62b28145782ed6404c8b0026bd37dd77f47e37 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 14 Mar 2018 16:22:08 +0100 Subject: [PATCH] Primitive Shell --- Basis/Basis.ml | 10 ++-- Basis/ContractedShell.ml | 95 +++++++++-------------------------- Basis/ContractedShell.mli | 8 +-- Basis/ContractedShellPair.mli | 8 +-- Basis/PrimitiveShell.ml | 82 ++++++++++++++++++++++++++++++ Basis/PrimitiveShell.mli | 54 ++++++++++++++++++++ 6 files changed, 173 insertions(+), 84 deletions(-) create mode 100644 Basis/PrimitiveShell.ml create mode 100644 Basis/PrimitiveShell.mli diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 38c7741..6c6af5d 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -6,6 +6,7 @@ type t = module Cs = ContractedShell module Gb = GeneralBasis +module Ps = PrimitiveShell (** Returns an array of the basis set per atom *) @@ -14,12 +15,11 @@ let of_nuclei_and_general_basis n b = Array.map (fun (e, center) -> List.assoc e b |> Array.map (fun (totAngMom, shell) -> - let expo = Array.map (fun Gb.{exponent ; coefficient} -> - exponent) shell - and coef = Array.map (fun Gb.{exponent ; coefficient} -> - coefficient) shell + let lc = + Array.map (fun Gb.{exponent ; coefficient} -> + coefficient, Ps.make totAngMom center exponent) shell in - Cs.make ~expo ~coef ~totAngMom ~center ~index:0) + Cs.make lc) ) n |> Array.to_list |> Array.concat diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml index ba16751..ad10c0f 100644 --- a/Basis/ContractedShell.ml +++ b/Basis/ContractedShell.ml @@ -7,59 +7,46 @@ type t = { coef : float array; (** Array of contraction coefficients {% $d_i$ %} *) center : Coordinate.t; (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %} *) totAngMom : AngularMomentum.t; (** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *) - size : int; (** Number of contracted functions, {% $m$ %} in the formula *) norm_coef : float array; (** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *) norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *) index : int; (** Index in the basis set, represented as an array of contracted shells. *) } module Am = AngularMomentum +module Ps = PrimitiveShell -let compute_norm_coef expo 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 - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2)) - ) expo - else - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot) - ) expo - in - Array.map (fun x -> let f a = x *. (factor a) in f) expo +let make ?(index=0) lc = + assert (Array.length lc > 0); + let coef = Array.map fst lc + and prim = Array.map snd lc + in -let make ~index ~expo ~coef ~center ~totAngMom = - assert (Array.length expo = Array.length coef); - assert (Array.length expo > 0); - let norm_coef_func = - compute_norm_coef expo totAngMom + let center = Ps.center prim.(0) in + let rec unique_center = function + | 0 -> true + | i -> if Ps.center prim.(i) = center then unique_center (i-1) else false in - let powers = - Am.zkey_array (Am.Singlet totAngMom) + if not (unique_center (Array.length prim - 1)) then + invalid_arg "ContractedShell.make Coordinate.t differ"; + + let totAngMom = Ps.totAngMom prim.(0) in + let rec unique_angmom = function + | 0 -> true + | i -> if Ps.totAngMom prim.(i) = totAngMom then unique_angmom (i-1) else false in + if not (unique_angmom (Array.length prim - 1)) then + invalid_arg "ContractedShell.make: AngularMomentum.t differ"; + + let expo = Array.map Ps.expo prim in + let norm_coef = - Array.map (fun f -> f [| Am.to_int totAngMom ; 0 ; 0 |]) norm_coef_func + Array.map Ps.norm_coef prim in - let norm_coef_scale = - Array.map (fun a -> - (norm_coef_func.(0) (Zkey.to_int_array a)) /. norm_coef.(0) - ) powers + let norm_coef_scale = Ps.norm_coef_scale prim.(0) in - { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ; + { index ; expo ; coef ; center ; totAngMom ; norm_coef ; norm_coef_scale } @@ -80,38 +67,6 @@ let to_string s = |> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) -(** 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 compute_norm_coef expo 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 - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2)) - ) expo - else - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot) - ) expo - in - Array.map (fun x -> let f a = x *. factor a in f) expo - let expo x = x.expo diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli index 5ebe8dd..7613d7a 100644 --- a/Basis/ContractedShell.mli +++ b/Basis/ContractedShell.mli @@ -33,12 +33,8 @@ type t val to_string : t -> string (** Pretty-printing of the contracted shell in a string *) -val make : - index:int -> - expo:float array -> - coef:float array -> - center:Coordinate.t -> totAngMom:AngularMomentum.t -> t -(** Creates a contracted shell *) +val make : ?index:int -> (float * PrimitiveShell.t) array -> t +(** Creates a contracted shell from a list of coefficients and primitives. *) val with_index : t -> int -> t (** Returns a copy of the contracted shell with a modified index *) diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli index 36f5021..fbf4cf9 100644 --- a/Basis/ContractedShellPair.mli +++ b/Basis/ContractedShellPair.mli @@ -61,8 +61,10 @@ val totAngMomInt : t -> int val monocentric : t -> bool (** If true, the two contracted shells have the same center. *) -val hash : 'a array -> int array -val cmp : 'a array -> 'a array -> int -val equivalent : 'a array -> 'a array -> bool +(* +val hash : t -> int array +val cmp : t -> t -> int +val equivalent : t -> t -> bool val unique : 'a array array array -> 'a array list val indices : 'a array array array -> (int array, int * int) Hashtbl.t +*) diff --git a/Basis/PrimitiveShell.ml b/Basis/PrimitiveShell.ml new file mode 100644 index 0000000..65877bd --- /dev/null +++ b/Basis/PrimitiveShell.ml @@ -0,0 +1,82 @@ +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_coef = + 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_coef + ) powers ) + 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_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) diff --git a/Basis/PrimitiveShell.mli b/Basis/PrimitiveShell.mli new file mode 100644 index 0000000..9efa8ca --- /dev/null +++ b/Basis/PrimitiveShell.mli @@ -0,0 +1,54 @@ +(** Set of Gaussians with a given {!AngularMomentum.t} + +{% \\[ +g(r) = (x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} \exp \left( -\alpha |r-R_A|^2 \right) + \\] %} + +where: + +- {% $n_x + n_y + n_z = l$ %}, the total angular momentum + +- {% $\alpha$ %} is the exponent + +*) + +type t + +val to_string : t -> string +(** Pretty-printing of the primitive shell in a string. *) + +val make : AngularMomentum.t -> Coordinate.t -> float -> t +(** Creates a primitive shell from the total angular momentum, the coordinates of the + center and the exponent. *) + +val expo : t -> float +(** Returns the exponent {% $\alpha$ %}. *) + +val center : t -> Coordinate.t +(** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *) + +val totAngMom : t -> AngularMomentum.t +(** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *) + +val norm_coef : t -> float +(** Normalization coefficient of the shell: + + {% \\[ + \mathcal{N} = \sqrt{\iiint \left[ (x-X_A)^{l} + \exp (-\alpha |r-R_A|^2) \right]^2 \, dx\, dy\, dz} + \\] %} +*) + +val norm_coef_scale : t -> float array +(** Scaling factors adjusting the normalization coefficient for the. + particular powers of {% $x,y,z$ %}. They are given in the same order as + [AngularMomentum.zkey_array totAngMom]: + + {% \\[ + f = \frac{1}{\mathcal{N}} \sqrt{\iiint [g(r)]^2 \, d^3r} + \\] %} +*) + +val size_of_shell : t -> int +(** Number of functions in the shell. *) +