10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-23 04:43:32 +01:00
QCaml/gaussian/lib/contracted_shell.ml

145 lines
3.5 KiB
OCaml
Raw Normal View History

2024-01-17 10:30:24 +01:00
[@@@landmark "auto-off"]
2020-10-09 09:47:57 +02:00
open Common
2018-02-23 18:44:31 +01:00
type t = {
2018-03-20 14:11:31 +01:00
expo : float array;
coef : float array;
norm_coef : float array;
norm_coef_scale : float array;
2020-09-26 12:02:53 +02:00
prim : Primitive_shell.t array;
2019-02-26 12:47:23 +01:00
center : Coordinate.t;
2020-09-26 12:02:53 +02:00
ang_mom : Angular_momentum.t;
2019-02-26 12:47:23 +01:00
index : int;
2018-02-23 18:44:31 +01:00
}
2020-09-26 12:02:53 +02:00
module Am = Angular_momentum
2018-03-15 15:25:49 +01:00
module Co = Coordinate
2020-09-26 12:02:53 +02:00
module Ps = Primitive_shell
2018-02-23 18:44:31 +01:00
2018-03-13 18:24:00 +01:00
2024-01-17 10:30:24 +01:00
let[@landmark] make ?(index=0) lc =
2018-03-14 16:22:08 +01:00
assert (Array.length lc > 0);
2018-02-23 18:44:31 +01:00
2018-03-14 16:22:08 +01:00
let coef = Array.map fst lc
and prim = Array.map snd lc
in
2018-02-23 18:44:31 +01:00
2018-03-14 16:22:08 +01:00
let center = Ps.center prim.(0) in
let rec unique_center = function
| 0 -> true
2019-09-10 18:39:14 +02:00
| i -> if Ps.center prim.(i) = center then (unique_center [@tailcall]) (i-1) else false
2018-02-23 18:44:31 +01:00
in
2018-03-14 16:22:08 +01:00
if not (unique_center (Array.length prim - 1)) then
invalid_arg "ContractedShell.make Coordinate.t differ";
2024-01-17 10:30:24 +01:00
2018-03-21 15:01:39 +01:00
let ang_mom = Ps.ang_mom prim.(0) in
2018-03-14 16:22:08 +01:00
let rec unique_angmom = function
| 0 -> true
2019-09-10 18:39:14 +02:00
| i -> if Ps.ang_mom prim.(i) = ang_mom then (unique_angmom [@tailcall]) (i-1) else false
2018-02-23 18:44:31 +01:00
in
2018-03-14 16:22:08 +01:00
if not (unique_angmom (Array.length prim - 1)) then
invalid_arg "ContractedShell.make: AngularMomentum.t differ";
2024-01-17 10:30:24 +01:00
2018-03-20 15:16:24 +01:00
let expo = Array.map Ps.exponent prim in
2018-03-14 16:22:08 +01:00
2018-02-23 18:44:31 +01:00
let norm_coef =
2018-03-20 15:16:24 +01:00
Array.map Ps.normalization prim
2018-02-23 18:44:31 +01:00
in
2018-03-20 15:16:24 +01:00
let norm_coef_scale = Ps.norm_scales prim.(0)
2018-02-23 18:44:31 +01:00
in
2018-03-21 15:01:39 +01:00
{ index ; expo ; coef ; center ; ang_mom ; norm_coef ;
2018-03-14 21:51:27 +01:00
norm_coef_scale ; prim }
2018-02-23 18:44:31 +01:00
let with_index a i =
{ a with index = i }
2018-03-20 14:11:31 +01:00
let exponents x = x.expo
2018-02-23 18:44:31 +01:00
2018-03-20 14:11:31 +01:00
let coefficients x = x.coef
2018-03-13 18:56:28 +01:00
let center x = x.center
2018-03-21 15:01:39 +01:00
let ang_mom x = x.ang_mom
2018-03-13 18:56:28 +01:00
2018-03-20 14:11:31 +01:00
let size x = Array.length x.prim
2018-03-13 18:56:28 +01:00
2018-03-20 14:11:31 +01:00
let normalizations x = x.norm_coef
2018-03-13 18:56:28 +01:00
2018-03-20 14:11:31 +01:00
let norm_scales x = x.norm_coef_scale
2018-03-13 18:56:28 +01:00
let index x = x.index
let size_of_shell x = Array.length x.norm_coef_scale
2018-03-14 21:51:27 +01:00
2018-03-20 14:11:31 +01:00
let primitives x = x.prim
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let zkey_array x = Ps.zkey_array x.prim.(0)
2024-01-17 10:30:24 +01:00
let[@landmark] values t point =
2020-04-17 00:36:50 +02:00
(* Radial part *)
let r = Co.( point |- t.center ) in
let r2 = Co.dot r r in
let radial =
let rec aux accu = function
| -1 -> accu
2024-01-17 10:30:24 +01:00
| i -> let new_accu =
2020-04-17 00:36:50 +02:00
t.norm_coef.(i) *. t.coef.(i) *. exp(-. t.expo.(i) *. r2) +. accu
in aux new_accu (i-1)
in
aux 0. (Array.length t.expo - 1)
in
(* Angular part *)
let n = Am.to_int t.ang_mom in
let x = Array.create_float (n+1) in
let y = Array.create_float (n+1) in
let z = Array.create_float (n+1) in
2024-01-17 10:30:24 +01:00
let fill arr v =
2020-04-17 00:36:50 +02:00
arr.(0) <- 1.;
for i=1 to n do
arr.(i) <- arr.(i-1) *. v
done;
in
fill x r.x; fill y r.y; fill z r.z;
let powers =
Am.zkey_array (Am.Singlet t.ang_mom)
in
Array.mapi (fun i a ->
let p = Zkey.to_int_array a in
2024-01-17 10:30:24 +01:00
t.norm_coef_scale.(i) *. x.(p.(0)) *. y.(p.(1)) *. z.(p.(2)) *. radial
2020-04-17 00:36:50 +02:00
) powers
2024-01-17 10:30:24 +01:00
2018-03-15 15:25:49 +01:00
(** {2 Printers} *)
2018-03-16 00:23:47 +01:00
open Format
2020-09-26 12:02:53 +02:00
(*
2018-03-16 00:23:47 +01:00
let pp_debug ppf x =
2018-03-15 15:25:49 +01:00
fprintf ppf "@[<2>{@ ";
fprintf ppf "@[<2>expo =@ %a ;@]@ " pp_float_array_size x.expo;
fprintf ppf "@[<2>coef =@ %a ;@]@ " pp_float_array_size x.coef;
fprintf ppf "@[<2>center =@ %a ;@]@ " Co.pp_angstrom x.center;
2018-03-21 15:01:39 +01:00
fprintf ppf "@[<2>ang_mom =@ %a ;@]@ " Am.pp_string x.ang_mom;
2018-03-15 15:25:49 +01:00
fprintf ppf "@[<2>norm_coef =@ %a ;@]@ " pp_float_array_size x.norm_coef;
fprintf ppf "@[<2>norm_coef_scale =@ %a ;@]@ " pp_float_array_size x.norm_coef_scale;
fprintf ppf "@[<2>index =@ %d ;@]@ " x.index;
fprintf ppf "}@,@]"
2020-09-26 12:02:53 +02:00
*)
2018-03-15 15:25:49 +01:00
2018-03-16 00:23:47 +01:00
let pp ppf s =
2018-03-21 15:01:39 +01:00
(match s.ang_mom with
2018-03-16 00:23:47 +01:00
| Am.S -> fprintf ppf "@[%3d@] " (s.index+1)
| _ -> fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+(Array.length s.norm_coef_scale))
);
2018-03-21 15:01:39 +01:00
fprintf ppf "@[%a@ %a@]@[" Am.pp_string s.ang_mom Co.pp s.center;
2018-03-16 00:23:47 +01:00
Array.iter2 (fun e c -> fprintf ppf "@[%16.8e %16.8e@]@;" e c) s.expo s.coef;
fprintf ppf "@]"