#+begin_src elisp tangle: no :results none :exports none
(setq pwd (file-name-directory buffer-file-name))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq lib (concat pwd "lib/"))
(setq testdir (concat pwd "test/"))
(setq mli (concat lib name ".mli"))
(setq ml  (concat lib name ".ml"))
(setq test-ml  (concat testdir name ".ml"))
(org-babel-tangle)
#+end_src 

* Atomic shell
  :PROPERTIES:
  :header-args: :noweb yes :comments both
  :END:

  Set of contracted Gaussians differing only by the powers of $x$, $y$ and $z$, with a
  constant ~Angular_momentum.t~, all centered on the same point.

  In other words, it is the set of all contracted shells sharing the same center.

  \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) \\
          & = (x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, \exp \left( -\alpha_{ij} |r-R_A|^2 \right)
  \end{align*}

  where:

  - $g_{ij\,n_x,n_y,n_z}(r)$ is the $i$-th ~PrimitiveShell.t~ of the $j$-th ~Contracted_shell.t~ 
  - $n_x + n_y + n_z = l$, the total angular momentum
  - $\alpha_{ij}$ are the exponents (tabulated) of the $j$-th ~Contracted_shell.t~ 
  - $d_{ij}$ are the contraction coefficients of the $j$-th ~Contracted_shell.t~ 
  - $\mathcal{N}_{ij}$ is the normalization coefficient of the $i$-th primitive shell
    (~PrimitiveShell.norm_coef~) of the $j$-th ~Contracted_shell.t~ 
  - $f(n_x,n_y,n_z)$ is a scaling factor adjusting the normalization coefficient for the
    particular powers of $x,y,z$ (~PrimitiveShell.norm_coef_scale~)


** Type

   #+begin_src ocaml :tangle (eval mli)
type t

open Common
   #+end_src

   #+begin_src ocaml :tangle (eval ml) :exports none
open Common

type t = {
  expo      : float array array;
  coef      : float array array;
  norm_coef : float array array;
  norm_coef_scale : float array;
  contr     : Contracted_shell.t array;
  index     : int;
  center    : Coordinate.t;
  ang_mom   : Angular_momentum.t;
}

module Am = Angular_momentum
module Co = Coordinate
module Cs = Contracted_shell
   #+end_src

** Access

   #+begin_src ocaml :tangle (eval mli)
val index : t -> int
(** Index in the basis set, represented as an array of contracted shells. *)

val center : t -> Coordinate.t
(** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *)

val ang_mom : t -> Angular_momentum.t
(** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *)

val size : t -> int
(** Number of contracted functions, {% $n$ %} in the definition. *)

val contracted_shells: t -> Contracted_shell.t array
(** Array of contracted gaussians *)

val exponents : t -> float array array
(** Array of exponents {% $\alpha_{ij}$ %}. The first index is the index of
    the contracted function, and the second index is the index of the primitive.
*)

val coefficients : t -> float array array
(** Array of contraction coefficients {% $d_{ij}$ %}. The first index is the index of
    the contracted function, and the second index is the index of the primitive.
*)

val normalizations : t -> float array array
(** Normalization coefficients {% $\mathcal{N}_{ij}$ %}. The first index is the index of
    the contracted function, and the second index is the index of the primitive.
*)

val norm_scales : t -> float array
(** Scaling factors {% $f(n_x,n_y,n_z)$ %}, given in the same order as
    [Angular_momentum.zkey_array ang_mom]. *)

val size_of_shell : t -> int
(** Number of contracted functions in the shell: length of {!norm_coef_scale}. *)

   #+end_src

 | ~index~              | Index in the basis set, represented as an array of contracted shells.                                                                                       |
 | ~center~             | Coordinate of the center $\mathbf{A} = (X_A,Y_A,Z_A)$.                                                                                                      |
 | ~ang_mom~            | Total angular momentum : $l = n_x + n_y + n_z$.                                                                                                             |
 | ~size~               | Number of contracted functions, $n$ in the definition.                                                                                                      |
 | ~contracted_shells:~ | Array of contracted gaussians                                                                                                                               |
 | ~exponents~          | Array of exponents $\alpha_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive.              |
 | ~coefficients~       | Array of contraction coefficients $d_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive.    |
 | ~normalizations~     | Normalization coefficients $\mathcal{N}_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive. |
 | ~norm_scales~        | Scaling factors $f(n_x,n_y,n_z)$, given in the same order as ~Angular_momentum.zkey_array ang_mom~.                                                         |
 | ~size_of_shell~      | Number of contracted functions in the shell: length of ~norm_coef_scale~.                                                                                   |
 
   #+begin_example

   #+end_example

   #+begin_src ocaml :tangle (eval ml) :exports none
let exponents x = x.expo

let coefficients x = x.coef

let center x = x.center

let ang_mom x = x.ang_mom

let size x = Array.length x.contr

let normalizations x = x.norm_coef

let norm_scales x = x.norm_coef_scale

let index x = x.index

let size_of_shell x = Array.length x.norm_coef_scale

let contracted_shells x = x.contr
   #+end_src

** Creation

   #+begin_src ocaml :tangle (eval mli)
val make : ?index:int -> Contracted_shell.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. *)
   #+end_src

   | ~make~       | Creates a contracted shell from a list of coefficients and primitives. |
   | ~with_index~ | Returns a copy of the contracted shell with a modified index.          |

   #+begin_example

   #+end_example

   #+begin_src ocaml :tangle (eval ml) :exports none
let make ?(index=0) contr = 
  assert (Array.length contr > 0);

  let coef = Array.map Cs.coefficients contr
  and expo = Array.map Cs.exponents contr
  in

  let center = Cs.center contr.(0) in
  let rec unique_center = function
  | 0 -> true
  | i -> if Cs.center contr.(i) = center then (unique_center [@tailcall]) (i-1) else false
  in
  if not (unique_center (Array.length contr - 1)) then
    invalid_arg "ContractedAtomicShell.make Coordinate.t differ";
    
  let ang_mom = Cs.ang_mom contr.(0) in
  let rec unique_angmom = function
  | 0 -> true
  | i -> if Cs.ang_mom contr.(i) = ang_mom then (unique_angmom [@tailcall]) (i-1) else false
  in
  if not (unique_angmom (Array.length contr - 1)) then
    invalid_arg "Contracted_shell.make: Angular_momentum.t differ";
    
  let norm_coef =
    Array.map Cs.normalizations contr 
  in
  let norm_coef_scale = Cs.norm_scales contr.(0)
  in
  { index ; expo ; coef ; center ; ang_mom ; norm_coef ;
    norm_coef_scale ; contr  }


let with_index a i =
  { a with index = i }
   #+end_src


** Printers

   #+begin_src ocaml :tangle (eval mli)
val pp : Format.formatter -> t -> unit
   #+end_src

   #+begin_src ocaml :tangle (eval ml) :exports none
let pp ppf s =
  let open Format in
  fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+ (size_of_shell s)*(size s));
  fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.ang_mom Co.pp s.center;
  Array.iter2 (fun e_arr c_arr -> fprintf ppf "@[<v>";
    Array.iter2 (fun e c -> fprintf ppf "@[%16.8e  %16.8e@]@;" e c)
    e_arr c_arr;
    fprintf ppf "@;@]") s.expo s.coef;
  fprintf ppf "@]"
   #+end_src