QCaml/Basis/Contracted_shell.ml

open Util

type t = {
  expo      : float array;
  coef      : float array;
  center    : Coordinate.t;
  totAngMom : Angular_momentum.t;
  size      : int;
  norm_coef : (int array -> float) array;
  indice    : int;
  powers   : Zkey.t array;
}

let size a = a.size
let expo a i = a.expo.(i)
let coef a i = a.coef.(i)
let center a = a.center
let totAngMom a = a.totAngMom
let norm_coef a i = a.norm_coef.(i)
let indice a  = a.indice
let powers a = a.powers


let to_string s =
  let coord =
    Coordinate.to_Bohr s.center
  in
  let open Printf in
  (match s.totAngMom with
   | Angular_momentum.S -> sprintf "%3d    " (s.indice+1)
   | _                  -> sprintf "%3d-%-3d" (s.indice+1) (s.indice+(Array.length s.powers))
  ) ^
  ( sprintf "%1s %8.3f %8.3f %8.3f " (Angular_momentum.to_string s.totAngMom)
      (Coordinate.x coord) (Coordinate.y coord) (Coordinate.z coord) ) ^
  (Array.map2 (fun e c -> sprintf "%16.8e  %16.8e" e c) s.expo s.coef
   |> 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.
*)
let compute_norm_coef s =
  let atot =
    Angular_momentum.to_int s.totAngMom
  in
  Array.mapi (fun i alpha ->
      let alpha_2 =
        alpha +. alpha
      in
      let c =
        (alpha_2 *. pi_inv)**(1.5) *. (pow (alpha_2 +. alpha_2) atot)
      in
      let result a = 
        let dfa = Array.map (fun j ->
            ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
          ) a
        in sqrt (c *.  dfa.(0) *.dfa.(1) *. dfa.(2))
      in
      result
    ) s.expo

let create ~indice ~expo ~coef ~center ~totAngMom =
  assert (Array.length expo = Array.length coef);
  assert (Array.length expo > 0);
  let tmp =
    { indice ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef = [||];
      powers = Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom) }
  in
  { tmp with norm_coef = compute_norm_coef tmp }