10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 04:13:33 +01:00
This commit is contained in:
Anthony Scemama 2018-03-20 15:16:24 +01:00
parent 5d73ec5144
commit 0da6453453
17 changed files with 153 additions and 154 deletions

View File

@ -90,15 +90,12 @@ let pp_debug ppf x =
fprintf ppf "}@,@]" fprintf ppf "}@,@]"
let pp ppf s = let pp ppf s =
(match s.totAngMom with fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+ (size_of_shell s)*(size s));
| Am.S -> fprintf ppf "@[%3d@] " (s.index+1)
| _ -> fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+(Array.length s.norm_coef_scale))
);
fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.totAngMom Co.pp s.center; fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.totAngMom Co.pp s.center;
Array.iter2 (fun e_arr c_arr -> fprintf ppf "[@[<v>"; Array.iter2 (fun e_arr c_arr -> fprintf ppf "@[<v>";
Array.iter2 (fun e c -> fprintf ppf "@[%16.8e %16.8e@]@;" e c) Array.iter2 (fun e c -> fprintf ppf "@[%16.8e %16.8e@]@;" e c)
e_arr c_arr; e_arr c_arr;
fprintf ppf "@]]@;") s.expo s.coef; fprintf ppf "@;@]") s.expo s.coef;
fprintf ppf "@]" fprintf ppf "@]"

View File

@ -1,6 +1,8 @@
(** Set of contracted Gaussians differing only by the powers of x, y and z, with a (** Set of contracted Gaussians differing only by the powers of x, y and z, with a
constant {!AngularMomentum.t}, all centered on the same center. constant {!AngularMomentum.t}, all centered on the same center.
In other words, it is the set of all contracted shells sharing the same center.
{% {%
\begin{align*} \begin{align*}
\chi_{n_x,n_y,n_z}(r) & = f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, g_{ij\,n_x,n_y,n_z}(r) \\ \chi_{n_x,n_y,n_z}(r) & = f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, g_{ij\,n_x,n_y,n_z}(r) \\

View File

@ -2,10 +2,10 @@ type t =
{ {
size : int; size : int;
contracted_shells : ContractedShell.t array ; contracted_shells : ContractedShell.t array ;
contracted_atomic_shells : ContractedAtomicShell.t array lazy_t; atomic_shells : AtomicShell.t array lazy_t;
} }
module Ca = ContractedAtomicShell module As = AtomicShell
module Cs = ContractedShell module Cs = ContractedShell
module Gb = GeneralBasis module Gb = GeneralBasis
module Ps = PrimitiveShell module Ps = PrimitiveShell
@ -30,7 +30,7 @@ let of_nuclei_and_general_basis nucl bas =
|> Array.to_list |> Array.to_list
|> Array.concat |> Array.concat
in in
let contracted_atomic_shells = lazy( let atomic_shells = lazy(
let uniq_center_angmom = let uniq_center_angmom =
Array.map (fun x -> Cs.center x, Cs.totAngMom x) contracted_shells Array.map (fun x -> Cs.center x, Cs.totAngMom x) contracted_shells
|> Array.to_list |> Array.to_list
@ -44,24 +44,24 @@ let of_nuclei_and_general_basis nucl bas =
List.filter (fun x -> Cs.center x = center && Cs.totAngMom x = totAngMom) csl List.filter (fun x -> Cs.center x = center && Cs.totAngMom x = totAngMom) csl
|> Array.of_list |> Array.of_list
in in
Ca.make ~index:(Cs.index a.(0)) a As.make ~index:(Cs.index a.(0)) a
) uniq_center_angmom ) uniq_center_angmom
|> List.sort (fun x y -> compare (Ca.index x) (Ca.index y)) |> List.sort (fun x y -> compare (As.index x) (As.index y))
|> Array.of_list |> Array.of_list
) in ) in
{ contracted_shells ; contracted_atomic_shells ; size = !index_ } { contracted_shells ; atomic_shells ; size = !index_ }
let size x = x.size let size x = x.size
let contracted_atomic_shells x = Lazy.force x.contracted_atomic_shells let atomic_shells x = Lazy.force x.atomic_shells
let contracted_shells x = x.contracted_shells let contracted_shells x = x.contracted_shells
let to_string b = let to_string b =
let b = contracted_atomic_shells b in let b = atomic_shells b in
let line =" let line ="
----------------------------------------------------------------------- -----------------------------------------------------------------------
" in " in
@ -75,7 +75,7 @@ let to_string b =
----------------------------------------------------------------------- -----------------------------------------------------------------------
" "
^ ^
( Array.map (fun p -> Format.(fprintf str_formatter "%a" Ca.pp p; ( Array.map (fun p -> Format.(fprintf str_formatter "%a" As.pp p;
flush_str_formatter ())) b flush_str_formatter ())) b
|> Array.to_list |> Array.to_list
|> String.concat line |> String.concat line

View File

@ -22,7 +22,7 @@ val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> string -> t
val size : t -> int val size : t -> int
(** Number of contracted basis functions. *) (** Number of contracted basis functions. *)
val contracted_atomic_shells : t -> ContractedAtomicShell.t array val atomic_shells : t -> AtomicShell.t array
(** Returns the contracted basis functions per atom. *) (** Returns the contracted basis functions per atom. *)
val contracted_shells : t -> ContractedShell.t array val contracted_shells : t -> ContractedShell.t array

View File

@ -40,12 +40,12 @@ let make ?(index=0) lc =
if not (unique_angmom (Array.length prim - 1)) then if not (unique_angmom (Array.length prim - 1)) then
invalid_arg "ContractedShell.make: AngularMomentum.t differ"; invalid_arg "ContractedShell.make: AngularMomentum.t differ";
let expo = Array.map Ps.expo prim in let expo = Array.map Ps.exponent prim in
let norm_coef = let norm_coef =
Array.map Ps.norm_coef prim Array.map Ps.normalization prim
in in
let norm_coef_scale = Ps.norm_coef_scale prim.(0) let norm_coef_scale = Ps.norm_scales prim.(0)
in in
{ index ; expo ; coef ; center ; totAngMom ; norm_coef ; { index ; expo ; coef ; center ; totAngMom ; norm_coef ;
norm_coef_scale ; prim } norm_coef_scale ; prim }

View File

@ -5,15 +5,15 @@ exception Null_contribution
type t = type t =
{ {
shell_a : ContractedShell.t; shell_a : ContractedShell.t;
shell_b : ContractedShell.t; shell_b : ContractedShell.t;
shell_pairs : PrimitiveShellPair.t array; shell_pairs : PrimitiveShellPair.t array;
coef : float array; coefficients : float array;
expo_inv : float array; exponents_inv : float array;
center_ab : Coordinate.t; (* A-B *) center_ab : Coordinate.t; (* A-B *)
norm_sq : float; (* |A-B|^2 *) norm_sq : float; (* |A-B|^2 *)
norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *) norm_scales : float array; (* norm_coef.(i) / norm_coef.(0) *)
totAngMomInt : int; (* Total angular Momentum *) totAngMom : AngularMomentum.t; (* Total angular Momentum *)
} }
@ -27,12 +27,7 @@ module Psp = PrimitiveShellPair
A contracted shell with N functions combined with a contracted A contracted shell with N functions combined with a contracted
shell with M functions generates a NxM array of shell pairs. shell with M functions generates a NxM array of shell pairs.
*) *)
let create ?(cutoff=1.e-32) s_a s_b = let make ?(cutoff=1.e-32) s_a s_b =
(*
Format.printf "@[<2>shell_a :@ %a@]@;" Cs.pp s_a;
Format.printf "@[<2>shell_b :@ %a@]@;" Cs.pp s_b;
*)
let make = Psp.create_make_of (Cs.primitives s_a).(0) (Cs.primitives s_b).(0) in let make = Psp.create_make_of (Cs.primitives s_a).(0) (Cs.primitives s_b).(0) in
@ -55,8 +50,8 @@ Format.printf "@[<2>shell_b :@ %a@]@;" Cs.pp s_b;
in in
let coef = Array.map (fun (c,y) -> c *. Psp.norm_coef y) shell_pairs let coefficients = Array.map (fun (c,y) -> c *. Psp.normalization y) shell_pairs
and expo_inv = Array.map (fun (_,y) -> Psp.expo_inv y) shell_pairs and exponents_inv = Array.map (fun (_,y) -> Psp.exponent_inv y) shell_pairs
in in
let shell_pairs = Array.map snd shell_pairs in let shell_pairs = Array.map snd shell_pairs in
if Array.length shell_pairs = 0 then if Array.length shell_pairs = 0 then
@ -64,23 +59,23 @@ Format.printf "@[<2>shell_b :@ %a@]@;" Cs.pp s_b;
else else
let root = shell_pairs.(0) in let root = shell_pairs.(0) in
Some { Some {
shell_a = s_a ; shell_b = s_b ; coef ; expo_inv ; shell_pairs ; shell_a = s_a ; shell_b = s_b ; coefficients ; exponents_inv ; shell_pairs ;
center_ab = Psp.a_minus_b root; center_ab = Psp.a_minus_b root;
norm_coef_scale = Psp.norm_coef_scale root; norm_scales = Psp.norm_scales root;
norm_sq=Psp.a_minus_b_sq root; norm_sq=Psp.a_minus_b_sq root;
totAngMomInt = Psp.totAngMom root |> Am.to_int; totAngMom = Psp.totAngMom root;
} }
let shell_a x = x.shell_a let shell_a x = x.shell_a
let shell_b x = x.shell_b let shell_b x = x.shell_b
let shell_pairs x = x.shell_pairs let shell_pairs x = x.shell_pairs
let coef x = x.coef let coefficients x = x.coefficients
let expo_inv x = x.expo_inv let exponents_inv x = x.exponents_inv
let center_ab x = x.center_ab let center_ab x = x.center_ab
let norm_sq x = x.norm_sq let norm_sq x = x.norm_sq
let totAngMomInt x = x.totAngMomInt let totAngMom x = x.totAngMom
let norm_coef_scale x = x.norm_coef_scale let norm_scales x = x.norm_scales
let monocentric x = Psp.monocentric x.shell_pairs.(0) let monocentric x = Psp.monocentric x.shell_pairs.(0)
@ -115,7 +110,7 @@ let cmp a b =
let of_contracted_shell_array basis = let of_contracted_shell_array basis =
Array.mapi (fun i shell_a -> Array.mapi (fun i shell_a ->
Array.mapi (fun j shell_b -> Array.mapi (fun j shell_b ->
create shell_a shell_b) make shell_a shell_b)
(Array.sub basis 0 (i+1)) (Array.sub basis 0 (i+1))
) basis ) basis

View File

@ -14,7 +14,7 @@ A contracted shell pair is a product of two {!ContractedShell.t}:
type t type t
val create : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option val make : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option
(** Creates an contracted shell pair {% $\varphi_{ab}$ %} from a contracted (** Creates an contracted shell pair {% $\varphi_{ab}$ %} from a contracted
shell {% $\chi_a$ %} (first argument) and a contracted shell {% $\chi_b$ %} shell {% $\chi_a$ %} (first argument) and a contracted shell {% $\chi_b$ %}
(second argument). (second argument).
@ -27,8 +27,12 @@ val create : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option
*) *)
val of_contracted_shell_array : ContractedShell.t array -> t option array array val of_contracted_shell_array : ContractedShell.t array -> t option array array
(** Creates all possible contracted shell pairs from the basis set. (** Creates all possible contracted shell pairs from a list of contracted shells.
If the shell pair is not significant, sets the value to [None]. If a shell pair is not significant, sets the value to [None]:
{[
(of_contracted_shell_array p).(i).(j) = create p.(i) p.(j)
]}
*) *)
val shell_a : t -> ContractedShell.t val shell_a : t -> ContractedShell.t
@ -46,9 +50,9 @@ val shell_pairs : t -> PrimitiveShellPair.t array
primitive functions used to build the contracted shell pair. primitive functions used to build the contracted shell pair.
*) *)
val coef : t -> float array val coefficients : t -> float array
val expo_inv : t -> float array val exponents_inv : t -> float array
val center_ab : t -> Coordinate.t val center_ab : t -> Coordinate.t
(* A-B *) (* A-B *)
@ -56,10 +60,10 @@ val center_ab : t -> Coordinate.t
val norm_sq : t -> float val norm_sq : t -> float
(* |A-B|^2 *) (* |A-B|^2 *)
val norm_coef_scale : t -> float array val norm_scales : t -> float array
(* norm_coef.(i) / norm_coef.(0) *) (* normalizations.(i) / normalizations.(0) *)
val totAngMomInt : t -> int val totAngMom : t -> AngularMomentum.t
(* Total angular Momentum *) (* Total angular Momentum *)
val monocentric : t -> bool val monocentric : t -> bool

View File

@ -24,7 +24,7 @@ let to_powers x =
(** Computes all the kinetic integrals of the contracted shell pair *) (** Computes all the kinetic integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t = let contracted_class shell_a shell_b : float Zmap.t =
match Csp.create shell_a shell_b with match Csp.make shell_a shell_b with
| Some shell_p -> | Some shell_p ->
begin begin
(* Pre-computation of integral class indices *) (* Pre-computation of integral class indices *)
@ -42,14 +42,14 @@ let contracted_class shell_a shell_b : float Zmap.t =
let center_ab = let center_ab =
Csp.center_ab shell_p Csp.center_ab shell_p
in in
let norm_coef_scale = let norm_coef_scales =
Csp.norm_coef_scale shell_p Csp.norm_scales shell_p
in in
for ab=0 to (Array.length sp - 1) for ab=0 to (Array.length sp - 1)
do do
let coef_prod = let coef_prod =
(Csp.coef shell_p).(ab) (Csp.coefficients shell_p).(ab)
in in
(** Screening on thr product of coefficients *) (** Screening on thr product of coefficients *)
if (abs_float coef_prod) > 1.e-4*.cutoff then if (abs_float coef_prod) > 1.e-4*.cutoff then
@ -58,12 +58,12 @@ let contracted_class shell_a shell_b : float Zmap.t =
Psp.center_minus_a sp.(ab) Psp.center_minus_a sp.(ab)
in in
let expo_inv = let expo_inv =
(Csp.expo_inv shell_p).(ab) (Csp.exponents_inv shell_p).(ab)
in in
let expo_a = let expo_a =
Ps.expo (Psp.shell_a sp.(ab)) Ps.exponent (Psp.shell_a sp.(ab))
and expo_b = and expo_b =
Ps.expo (Psp.shell_b sp.(ab)) Ps.exponent (Psp.shell_b sp.(ab))
in in
let xyz_of_int k = let xyz_of_int k =
@ -99,7 +99,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
let s = Array.init 3 f let s = Array.init 3 f
and k = Array.init 3 g and k = Array.init 3 g
in in
let norm = norm_coef_scale.(i) in let norm = norm_coef_scales.(i) in
let integral = chop norm (fun () -> let integral = chop norm (fun () ->
k.(0)*.s.(1)*.s.(2) +. k.(0)*.s.(1)*.s.(2) +.
s.(0)*.k.(1)*.s.(2) +. s.(0)*.k.(1)*.s.(2) +.

View File

@ -58,7 +58,7 @@ let of_basis_nuclei basis nuclei =
(* Pre-compute all shell pairs *) (* Pre-compute all shell pairs *)
let shell_pairs = let shell_pairs =
Array.mapi (fun i shell_a -> Array.map (fun shell_b -> Array.mapi (fun i shell_a -> Array.map (fun shell_b ->
ContractedShellPair.create shell_a shell_b) (Array.sub shell 0 (i+1)) ) shell ContractedShellPair.make shell_a shell_b) (Array.sub shell 0 (i+1)) ) shell
in in
(* Compute Integrals *) (* Compute Integrals *)

View File

@ -117,14 +117,14 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
(* Compute all integrals in the shell for each pair of significant shell pairs *) (* Compute all integrals in the shell for each pair of significant shell pairs *)
let norm_coef_scale_p = Csp.norm_coef_scale shell_p let norm_scales_p = Csp.norm_scales shell_p
in in
let sp = Csp.shell_pairs shell_p in let sp = Csp.shell_pairs shell_p in
for ab=0 to Array.length sp - 1 for ab=0 to Array.length sp - 1
do do
try try
begin begin
let coef_prod = (Csp.coef shell_p).(ab) in let coef_prod = (Csp.coefficients shell_p).(ab) in
(** Screening on the product of coefficients *) (** Screening on the product of coefficients *)
if abs_float coef_prod < 1.e-3 *. integrals_cutoff then if abs_float coef_prod < 1.e-3 *. integrals_cutoff then
@ -132,11 +132,11 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
let expo_pq_inv = let expo_pq_inv =
(Csp.expo_inv shell_p).(ab) (Csp.exponents_inv shell_p).(ab)
in in
let expo_b = let expo_b =
Ps.expo (Psp.shell_b sp.(ab)) Ps.exponent (Psp.shell_b sp.(ab))
in in
let center_ab = let center_ab =
@ -168,7 +168,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge
| _ -> | _ ->
let map = Zmap.create (2*maxm) in let map = Zmap.create (2*maxm) in
let norm_coef_scale = norm_coef_scale_p in let norm_scales = norm_scales_p in
(* Compute the integral class from the primitive shell quartet *) (* Compute the integral class from the primitive shell quartet *)
class_indices class_indices
|> Array.iteri (fun i key -> |> Array.iteri (fun i key ->
@ -177,7 +177,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
| Zkey.Six x -> x | Zkey.Six x -> x
| _ -> assert false | _ -> assert false
in in
let norm = norm_coef_scale.(i) in let norm = norm_scales.(i) in
let coef_prod = coef_prod *. norm in let coef_prod = coef_prod *. norm in
let integral = let integral =
hvrr_one_e hvrr_one_e

View File

@ -23,7 +23,7 @@ let to_powers x =
(** Computes all the overlap integrals of the contracted shell pair *) (** Computes all the overlap integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t = let contracted_class shell_a shell_b : float Zmap.t =
match Csp.create shell_a shell_b with match Csp.make shell_a shell_b with
| Some shell_p -> | Some shell_p ->
begin begin
@ -43,8 +43,8 @@ let contracted_class shell_a shell_b : float Zmap.t =
let center_ab = let center_ab =
Csp.center_ab shell_p Csp.center_ab shell_p
in in
let norm_coef_scale = let norm_coef_scales =
Csp.norm_coef_scale shell_p Csp.norm_scales shell_p
in in
(* Compute all integrals in the shell for each pair of significant shell pairs *) (* Compute all integrals in the shell for each pair of significant shell pairs *)
@ -58,13 +58,13 @@ let contracted_class shell_a shell_b : float Zmap.t =
for ab=0 to (Array.length sp - 1) for ab=0 to (Array.length sp - 1)
do do
let coef_prod = let coef_prod =
(Csp.coef shell_p).(ab) (Csp.coefficients shell_p).(ab)
in in
(** Screening on thr product of coefficients *) (** Screening on thr product of coefficients *)
if (abs_float coef_prod) > 1.e-3*.cutoff then if (abs_float coef_prod) > 1.e-3*.cutoff then
begin begin
let expo_inv = let expo_inv =
(Csp.expo_inv shell_p).(ab) (Csp.exponents_inv shell_p).(ab)
in in
let center_pa = let center_pa =
Psp.center_minus_a sp.(ab) Psp.center_minus_a sp.(ab)
@ -79,7 +79,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
(Co.get xyz center_ab, (Co.get xyz center_ab,
Co.get xyz center_pa) Co.get xyz center_pa)
in in
let norm = norm_coef_scale.(i) in let norm = norm_coef_scales.(i) in
let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
) class_indices ) class_indices

View File

@ -3,9 +3,9 @@ open Constants
open Coordinate open Coordinate
type t = { type t = {
expo : float; exponent : float;
norm_coef : float; normalization : float;
norm_coef_scale : float array lazy_t; norm_scales : float array lazy_t;
center : Coordinate.t; center : Coordinate.t;
totAngMom : AngularMomentum.t; totAngMom : AngularMomentum.t;
} }
@ -37,9 +37,9 @@ let compute_norm_coef alpha totAngMom =
in f in f
let make totAngMom center expo = let make totAngMom center exponent =
let norm_coef_func = let norm_coef_func =
compute_norm_coef expo totAngMom compute_norm_coef exponent totAngMom
in in
let norm = let norm =
1. /. norm_coef_func [| Am.to_int totAngMom ; 0 ; 0 |] 1. /. norm_coef_func [| Am.to_int totAngMom ; 0 ; 0 |]
@ -47,39 +47,40 @@ let make totAngMom center expo =
let powers = let powers =
Am.zkey_array (Am.Singlet totAngMom) Am.zkey_array (Am.Singlet totAngMom)
in in
let norm_coef_scale = lazy ( let norm_scales = lazy (
Array.map (fun a -> Array.map (fun a ->
(norm_coef_func (Zkey.to_int_array a)) *. norm (norm_coef_func (Zkey.to_int_array a)) *. norm
) powers ) ) powers )
in in
let norm_coef = 1. /. norm in let normalization = 1. /. norm in
{ expo ; norm_coef ; norm_coef_scale ; center ; totAngMom } { exponent ; normalization ; norm_scales ; center ; totAngMom }
let to_string s = let to_string s =
let coord = s.center in let coord = s.center in
Printf.sprintf "%1s %8.3f %8.3f %8.3f %16.8e" (Am.to_string s.totAngMom) 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 (get X coord) (get Y coord) (get Z coord) s.exponent
(** Normalization coefficient of contracted function i, which depends on the (** Normalization coefficient of contracted function i, which depends on the
exponent and the angular momentum. Two conventions can be chosen : a single exponent and the angular momentum. Two conventions can be chosen : a single
normalisation factor for all functions of the class, or a coefficient which normalization factor for all functions of the class, or a coefficient which
depends on the powers of x,y and z. depends on the powers of x,y and z.
Returns, for each contracted function, an array of functions taking as Returns, for each contracted function, an array of functions taking as
argument the [|x;y;z|] powers. argument the [|x;y;z|] powers.
*) *)
let expo x = x.expo let exponent x = x.exponent
let center x = x.center let center x = x.center
let totAngMom x = x.totAngMom let totAngMom x = x.totAngMom
let norm x = 1. /. x.norm_coef let norm x = 1. /. x.normalization
let norm_coef x = x.norm_coef let normalization x = x.normalization
let norm_coef_scale x = Lazy.force x.norm_coef_scale let norm_scales x = Lazy.force x.norm_scales
let size_of_shell x = Array.length (norm_scales x)
let size_of_shell x = Array.length (norm_coef_scale x)

View File

@ -27,7 +27,7 @@ val make : AngularMomentum.t -> Coordinate.t -> float -> t
(** Creates a primitive shell from the total angular momentum, the coordinates of the (** Creates a primitive shell from the total angular momentum, the coordinates of the
center and the exponent. *) center and the exponent. *)
val expo : t -> float val exponent : t -> float
(** Exponent {% $\alpha$ %}. *) (** Exponent {% $\alpha$ %}. *)
val center : t -> Coordinate.t val center : t -> Coordinate.t
@ -44,13 +44,13 @@ val norm : t -> float
\\] %} \\] %}
*) *)
val norm_coef : t -> float val normalization : t -> float
(** Normalization coefficient by which the shell has to be multiplied (** Normalization coefficient by which the shell has to be multiplied
to be normalized : to be normalized :
{% \\[ \mathcal{N} = \frac{1}{|| g_{l,0,0}(\mathbf{r}) ||} \\] %}. {% \\[ \mathcal{N} = \frac{1}{|| g_{l,0,0}(\mathbf{r}) ||} \\] %}.
*) *)
val norm_coef_scale : t -> float array val norm_scales : t -> float array
(** Scaling factors {% $f(n_x,n_y,n_z)$ %} adjusting the normalization coefficient (** Scaling factors {% $f(n_x,n_y,n_z)$ %} adjusting the normalization coefficient
for the powers of {% $x,y,z$ %}. The normalization coefficients of the for the powers of {% $x,y,z$ %}. The normalization coefficients of the
functions of the shell are given by {% $\mathcal{N}\times f$ %}. They are functions of the shell are given by {% $\mathcal{N}\times f$ %}. They are

View File

@ -3,14 +3,14 @@ open Constants
type t = { type t = {
expo : float; (* alpha + beta *) exponent : float; (* alpha + beta *)
expo_inv : float; (* 1/(alpha + beta) *) exponent_inv : float; (* 1/(alpha + beta) *)
center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *) center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *)
center_minus_a : Coordinate.t; (* P - A *) center_minus_a : Coordinate.t; (* P - A *)
a_minus_b : Coordinate.t; (* A - B *) a_minus_b : Coordinate.t; (* A - B *)
a_minus_b_sq : float; (* |A-B|^2 *) a_minus_b_sq : float; (* |A-B|^2 *)
norm_coef_scale : float array lazy_t; norm_scales : float array lazy_t;
norm_coef : float; (* norm_coef_a * norm_coef_b * g, with normalization : float; (* norm_coef_a * norm_coef_b * g, with
g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *) g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *)
totAngMom : AngularMomentum.t; totAngMom : AngularMomentum.t;
shell_a : PrimitiveShell.t; shell_a : PrimitiveShell.t;
@ -29,9 +29,9 @@ let hash a =
let equivalent a b = let equivalent a b =
a.expo = b.expo && a.exponent = b.exponent &&
a.totAngMom = b.totAngMom && a.totAngMom = b.totAngMom &&
a.norm_coef = b.norm_coef && a.normalization = b.normalization &&
a.center = b.center && a.center = b.center &&
a.center_minus_a = b.center_minus_a && a.center_minus_a = b.center_minus_a &&
a.a_minus_b = b.a_minus_b a.a_minus_b = b.a_minus_b
@ -51,10 +51,10 @@ let create_make_of p_a p_b =
Co.dot a_minus_b a_minus_b Co.dot a_minus_b a_minus_b
in in
let norm_coef_scale = lazy ( let norm_scales = lazy (
Array.map (fun v1 -> Array.map (fun v1 ->
Array.map (fun v2 -> v1 *. v2) (Ps.norm_coef_scale p_b) Array.map (fun v2 -> v1 *. v2) (Ps.norm_scales p_b)
) (Ps.norm_coef_scale p_a) ) (Ps.norm_scales p_a)
|> Array.to_list |> Array.to_list
|> Array.concat |> Array.concat
) in ) in
@ -66,42 +66,42 @@ let create_make_of p_a p_b =
function p_a -> function p_a ->
let norm_coef_a = let norm_coef_a =
Ps.norm_coef p_a Ps.normalization p_a
in in
let alpha_a = let alpha_a =
Co.( Ps.expo p_a |. Ps.center p_a ) Co.( Ps.exponent p_a |. Ps.center p_a )
in in
function p_b -> function p_b ->
let norm_coef = let normalization =
norm_coef_a *. Ps.norm_coef p_b norm_coef_a *. Ps.normalization p_b
in in
let expo = let exponent =
Ps.expo p_a +. Ps.expo p_b Ps.exponent p_a +. Ps.exponent p_b
in in
let expo_inv = 1. /. expo in let exponent_inv = 1. /. exponent in
let norm_coef = let normalization =
let argexpo = let argexpo =
Ps.expo p_a *. Ps.expo p_b *. a_minus_b_sq *. expo_inv Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv
in in
norm_coef *. (pi *. expo_inv)**1.5 *. exp (-. argexpo) normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo)
in in
function cutoff -> function cutoff ->
if abs_float norm_coef > cutoff then ( if abs_float normalization > cutoff then (
let beta_b = let beta_b =
Co.( Ps.expo p_b |. Ps.center p_b ) Co.( Ps.exponent p_b |. Ps.center p_b )
in in
let center = let center =
Co.(expo_inv |. (alpha_a |+ beta_b)) Co.(exponent_inv |. (alpha_a |+ beta_b))
in in
let center_minus_a = let center_minus_a =
@ -110,8 +110,8 @@ let create_make_of p_a p_b =
Some { Some {
totAngMom ; totAngMom ;
expo ; expo_inv ; center ; center_minus_a ; a_minus_b ; exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ;
a_minus_b_sq ; norm_coef ; norm_coef_scale ; shell_a = p_a; a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a;
shell_b = p_b } shell_b = p_b }
) )
@ -127,10 +127,10 @@ let make p_a p_b =
| None -> assert false | None -> assert false
let norm_coef_scale x = let norm_scales x =
Lazy.force x.norm_coef_scale Lazy.force x.norm_scales
let expo_inv x = x.expo_inv let exponent_inv x = x.exponent_inv
let monocentric x = let monocentric x =
Ps.center x.shell_a = Ps.center x.shell_b Ps.center x.shell_a = Ps.center x.shell_b
@ -144,9 +144,9 @@ let a_minus_b_sq x = x.a_minus_b_sq
let center_minus_a x = x.center_minus_a let center_minus_a x = x.center_minus_a
let norm_coef x = x.norm_coef let normalization x = x.normalization
let expo x = x.expo let exponent x = x.exponent
let center x = x.center let center x = x.center

View File

@ -69,20 +69,20 @@ val shell_a : t -> PrimitiveShell.t
val shell_b : t -> PrimitiveShell.t val shell_b : t -> PrimitiveShell.t
(** Returns the second primitive shell that was used to build the shell pair. *) (** Returns the second primitive shell that was used to build the shell pair. *)
val norm_coef : t -> float val normalization : t -> float
(** Normalization coefficient of the shell pair. *) (** Normalization coefficient of the shell pair. *)
val norm_coef_scale : t -> float array val norm_scales : t -> float array
(** Normalization factor, characteristic of the powers of x, y and z of (** Normalization factor, characteristic of the powers of x, y and z of
both shells of the pair. It is the outer product of the 2 both shells of the pair. It is the outer product of the 2
{!PrimitiveShell.norm_coef_scale} arrays of the shells consituting the {!PrimitiveShell.norm_coef_scale} arrays of the shells consituting the
pair. pair.
*) *)
val expo : t -> float val exponent : t -> float
(** Exponent of the Gaussian output of the Gaussian product : {% \\[ \alpha + \beta \\] %}*) (** Exponent of the Gaussian output of the Gaussian product : {% \\[ \alpha + \beta \\] %}*)
val expo_inv : t -> float val exponent_inv : t -> float
(** Inverse of the exponent : {% \\[ \sigma_P = \frac{1}{\alpha + \beta} \\] %}*) (** Inverse of the exponent : {% \\[ \sigma_P = \frac{1}{\alpha + \beta} \\] %}*)
val a_minus_b : t -> Coordinate.t val a_minus_b : t -> Coordinate.t

View File

@ -281,7 +281,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
and shell_d = Csp.shell_b shell_q and shell_d = Csp.shell_b shell_q
and sp = Csp.shell_pairs shell_p and sp = Csp.shell_pairs shell_p
in in
let maxm = Csp.totAngMomInt shell_p + Csp.totAngMomInt shell_q in let maxm = Am.(Csp.totAngMom shell_p + Csp.totAngMom shell_q |> to_int) in
(* Pre-computation of integral class indices *) (* Pre-computation of integral class indices *)
let class_indices = let class_indices =
@ -301,22 +301,22 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
(* Compute all integrals in the shell for each pair of significant shell pairs *) (* Compute all integrals in the shell for each pair of significant shell pairs *)
let norm_coef_scale_p_list = Array.to_list (Csp.norm_coef_scale shell_p) in let norm_coef_scale_p_list = Array.to_list (Csp.norm_scales shell_p) in
let norm_coef_scale_q = Csp.norm_coef_scale shell_q in let norm_coef_scale_q = Csp.norm_scales shell_q in
for ab=0 to (Array.length sp - 1) do for ab=0 to (Array.length sp - 1) do
let sp_ab = (Csp.shell_pairs shell_p).(ab) in let sp_ab = (Csp.shell_pairs shell_p).(ab) in
let c_ab = (Csp.coef shell_p).(ab) in let c_ab = (Csp.coefficients shell_p).(ab) in
let expo_b = Ps.expo (Psp.shell_b sp_ab) in let expo_b = Ps.exponent (Psp.shell_b sp_ab) in
let expo_inv_p = Psp.expo_inv sp_ab in let expo_inv_p = Psp.exponent_inv sp_ab in
let center_ab = Psp.a_minus_b sp_ab in let center_ab = Psp.a_minus_b sp_ab in
let center_pa = Psp.center_minus_a sp_ab in let center_pa = Psp.center_minus_a sp_ab in
for cd=0 to (Array.length (Csp.shell_pairs shell_q) - 1) do for cd=0 to (Array.length (Csp.shell_pairs shell_q) - 1) do
let sp_cd = (Csp.shell_pairs shell_q).(cd) in let sp_cd = (Csp.shell_pairs shell_q).(cd) in
let c_cd = (Csp.coef shell_q).(cd) in let c_cd = (Csp.coefficients shell_q).(cd) in
let coef_prod = c_ab *. c_cd in let coef_prod = c_ab *. c_cd in
(** Screening on the product of coefficients *) (** Screening on the product of coefficients *)
@ -326,7 +326,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in
let norm_pq_sq = Co.dot center_pq center_pq in let norm_pq_sq = Co.dot center_pq center_pq in
let expo_inv_q = Psp.expo_inv sp_cd in let expo_inv_q = Psp.exponent_inv sp_cd in
let expo_pq_inv = expo_inv_p +. expo_inv_q in let expo_pq_inv = expo_inv_p +. expo_inv_q in
let zero_m_array = let zero_m_array =
@ -342,7 +342,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
in in
contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
| _ -> | _ ->
let expo_d = Ps.expo (Psp.shell_b sp_cd) in let expo_d = Ps.exponent (Psp.shell_b sp_cd) in
let map_1d = Zmap.create (4*maxm) in let map_1d = Zmap.create (4*maxm) in
let map_2d = Zmap.create (Array.length class_indices) in let map_2d = Zmap.create (Array.length class_indices) in
let center_cd = Psp.a_minus_b sp_cd in let center_cd = Psp.a_minus_b sp_cd in
@ -429,8 +429,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
(** Computes all the two-electron integrals of the contracted shell quartet *) (** 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 contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
let shell_p = Csp.create ~cutoff shell_a shell_b let shell_p = Csp.make ~cutoff shell_a shell_b
and shell_q = Csp.create ~cutoff shell_c shell_d and shell_q = Csp.make ~cutoff shell_c shell_d
in in
match shell_p, shell_q with match shell_p, shell_q with
| Some shell_p, Some shell_q -> | Some shell_p, Some shell_q ->

View File

@ -563,7 +563,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
and sq = Csp.shell_pairs shell_q and sq = Csp.shell_pairs shell_q
in in
let maxm = let maxm =
Csp.totAngMomInt shell_p + Csp.totAngMomInt shell_q Am.(Csp.totAngMom shell_p + Csp.totAngMom shell_q |> to_int)
in in
(* Pre-computation of integral class indices *) (* Pre-computation of integral class indices *)
@ -587,11 +587,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let coef_max_p = let coef_max_p =
Array.fold_left (fun accu x -> Array.fold_left (fun accu x ->
if (abs_float x) > accu then (abs_float x) else accu) if (abs_float x) > accu then (abs_float x) else accu)
0. (Csp.coef shell_p) 0. (Csp.coefficients shell_p)
and coef_max_q = and coef_max_q =
Array.fold_left (fun accu x -> Array.fold_left (fun accu x ->
if (abs_float x) > accu then (abs_float x) else accu) if (abs_float x) > accu then (abs_float x) else accu)
0. (Csp.coef shell_q) 0. (Csp.coefficients shell_q)
in in
let rec build_list cutoff vec accu = function let rec build_list cutoff vec accu = function
@ -601,10 +601,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
else accu ) (k-1) else accu ) (k-1)
in in
let p_list = let p_list =
let vec = (Csp.coef shell_p) in let vec = (Csp.coefficients shell_p) in
build_list (cutoff /. coef_max_q) vec [] (Array.length vec - 1) build_list (cutoff /. coef_max_q) vec [] (Array.length vec - 1)
and q_list = and q_list =
let vec = (Csp.coef shell_q) in let vec = (Csp.coefficients shell_q) in
build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1) build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1)
in in
@ -630,16 +630,16 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
begin begin
try try
let expo_inv_p = let expo_inv_p =
Vec.init np (fun ab -> Psp.expo_inv sp.(ab-1)) Vec.init np (fun ab -> Psp.exponent_inv sp.(ab-1))
and expo_inv_q = and expo_inv_q =
Vec.init nq (fun cd -> Psp.expo_inv sq.(cd-1)) Vec.init nq (fun cd -> Psp.exponent_inv sq.(cd-1))
in in
let coef = let coef =
let result = Mat.make0 nq np in let result = Mat.make0 nq np in
Lacaml.D.ger Lacaml.D.ger
(Vec.of_array @@ filter_q (Csp.coef shell_q)) (Vec.of_array @@ filter_q (Csp.coefficients shell_q))
(Vec.of_array @@ filter_p (Csp.coef shell_p)) (Vec.of_array @@ filter_p (Csp.coefficients shell_p))
result; result;
result result
in in
@ -671,24 +671,24 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
| _ -> | _ ->
let coef = let coef =
let cp = filter_p (Csp.coef shell_p) let cp = filter_p (Csp.coefficients shell_p)
and cq = filter_q (Csp.coef shell_q) and cq = filter_q (Csp.coefficients shell_q)
in in
Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) )
in in
let expo_inv_p = let expo_inv_p =
Array.map (fun shell_ab -> Psp.expo_inv shell_ab) sp Array.map (fun shell_ab -> Psp.exponent_inv shell_ab) sp
and expo_inv_q = and expo_inv_q =
Array.map (fun shell_cd -> Psp.expo_inv shell_cd) sq Array.map (fun shell_cd -> Psp.exponent_inv shell_cd) sq
in in
let expo_b = let expo_b =
Array.map (fun shell_ab -> Ps.expo (Psp.shell_b shell_ab) ) sp Array.map (fun shell_ab -> Ps.exponent (Psp.shell_b shell_ab) ) sp
and expo_d = and expo_d =
Array.map (fun shell_cd -> Ps.expo (Psp.shell_b shell_cd) ) sq Array.map (fun shell_cd -> Ps.exponent (Psp.shell_b shell_cd) ) sq
in in
let norm_coef_scale_p = Csp.norm_coef_scale shell_p in let norm_coef_scale_p = Csp.norm_scales shell_p in
let center_pq = let center_pq =
let result = let result =
@ -792,7 +792,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let norm = let norm =
let norm_coef_scale_q = let norm_coef_scale_q =
Csp.norm_coef_scale shell_q Csp.norm_scales shell_q
in in
Array.to_list norm_coef_scale_p Array.to_list norm_coef_scale_p
|> List.map (fun v1 -> |> List.map (fun v1 ->
@ -870,8 +870,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
(** Computes all the two-electron integrals of the contracted shell quartet *) (** 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 contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
let shell_p = Csp.create ~cutoff shell_a shell_b let shell_p = Csp.make ~cutoff shell_a shell_b
and shell_q = Csp.create ~cutoff shell_c shell_d and shell_q = Csp.make ~cutoff shell_c shell_d
in in
match shell_p, shell_q with match shell_p, shell_q with
| Some shell_p, Some shell_q -> | Some shell_p, Some shell_q ->