QCaml/gaussian_basis/lib/contracted_shell.ml

open Qcaml_common

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

module Am = Angular_momentum
module Co = Coordinate
module Ps = Primitive_shell


let make ?(index=0) lc = 
  assert (Array.length lc > 0);

  let coef = Array.map fst lc
  and prim = Array.map snd lc
  in

  let center = Ps.center prim.(0) in
  let rec unique_center = function
  | 0 -> true
  | i -> if Ps.center prim.(i) = center then (unique_center [@tailcall]) (i-1) else false
  in
  if not (unique_center (Array.length prim - 1)) then
    invalid_arg "ContractedShell.make Coordinate.t differ";
    
  let ang_mom = Ps.ang_mom prim.(0) in
  let rec unique_angmom = function
  | 0 -> true
  | i -> if Ps.ang_mom prim.(i) = ang_mom then (unique_angmom [@tailcall]) (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.exponent prim in

  let norm_coef =
    Array.map Ps.normalization prim
  in
  let norm_coef_scale = Ps.norm_scales prim.(0)
  in
  { index ; expo ; coef ; center ; ang_mom ; norm_coef ;
    norm_coef_scale ; prim }


let with_index a i =
  { a with index = i }


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.prim

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 primitives x = x.prim

let zkey_array x = Ps.zkey_array x.prim.(0)

let values t point =
  (* 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
    | i -> let new_accu = 
             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
  let fill arr v = 
    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
    t.norm_coef_scale.(i) *. x.(p.(0)) *. y.(p.(1)) *. z.(p.(2)) *. radial 
    ) powers

  

(** {2 Printers} *)

open Format

(*
let pp_debug ppf x =
  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;
  fprintf ppf "@[<2>ang_mom =@ %a ;@]@ " Am.pp_string x.ang_mom;
  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 "}@,@]"
*)

let pp ppf s =
  (match s.ang_mom with
   | 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.ang_mom Co.pp s.center;
  Array.iter2 (fun e c -> fprintf ppf "@[%16.8e  %16.8e@]@;" e c) s.expo s.coef;
  fprintf ppf "@]"