QCaml/common/angular_momentum.org

#+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


* Angular Momentum
  :PROPERTIES:
  :header-args: :noweb yes :comments both
  :END:

  Azimuthal quantum number, repsesented as \( s,p,d,\dots \) .

** Type

   #+NAME: types
   #+begin_src ocaml :tangle (eval mli)
type t =
  | S | P | D | F | G | H | I | J | K | L | M | N | O
  | Int of int

exception Angular_momentum_error of string

type kind =
    Singlet of t
  | Doublet of (t * t)
  | Triplet of (t * t * t)
  | Quartet of (t * t * t * t)

   #+end_src

   An exception is raised when the ~Angular_momentum.t~ element can't
   be created.

   The ~kind~ is used to build shells, shell doublets, triplets or
   quartets, use in the two-electron operators.

   #+begin_src ocaml :tangle (eval ml) :exports none
<<types>>
open Powers
   #+end_src

** Conversions


*** ~of_char~

    #+begin_src ocaml :tangle (eval mli)
val of_char : char -> t
    #+end_src

    Returns an ~Angular_momentum.t~ when a shell is given as a character
    (case insensitive):

    #+begin_example
Angular_momentum.of_char 'p' -> Angular_momentum.P
    #+end_example

    #+begin_src ocaml :tangle (eval ml) :exports none
let of_char = function
  | 's' | 'S' -> S      | 'p' | 'P' -> P
  | 'd' | 'D' -> D      | 'f' | 'F' -> F
  | 'g' | 'G' -> G      | 'h' | 'H' -> H
  | 'i' | 'I' -> I      | 'j' | 'J' -> J
  | 'k' | 'K' -> K      | 'l' | 'L' -> L
  | 'm' | 'M' -> M      | 'n' | 'N' -> N
  | 'o' | 'O' -> O
  | c -> raise (Angular_momentum_error (String.make 1 c))
    #+end_src

*** ~to_string~

    #+begin_src ocaml :tangle (eval mli)
val to_string : t -> string
    #+end_src

    Converts the angular momentum into a string:

    #+begin_example
Angular_momentum.(to_string D) -> "D"
    #+end_example

    #+begin_src ocaml :tangle (eval ml) :exports none
let to_string = function
  | S -> "S"       | P -> "P"
  | D -> "D"       | F -> "F"
  | G -> "G"       | H -> "H"
  | I -> "I"       | J -> "J"
  | K -> "K"       | L -> "L"
  | M -> "M"       | N -> "N"
  | O -> "O"       | Int i -> string_of_int i
    #+end_src

*** ~to_char~

    #+begin_src ocaml :tangle (eval mli)
val to_char : t -> char
    #+end_src

    Converts the angular momentum into a char:

    #+begin_example
Angular_momentum.(to_char D) -> 'D'
    #+end_example

    #+begin_src ocaml :tangle (eval ml) :exports none
let to_char = function
  | S -> 'S'       | P -> 'P'
  | D -> 'D'       | F -> 'F'
  | G -> 'G'       | H -> 'H'
  | I -> 'I'       | J -> 'J'
  | K -> 'K'       | L -> 'L'
  | M -> 'M'       | N -> 'N'
  | O -> 'O'       | Int _ -> '_'
    #+end_src

*** ~to_int~

    #+begin_src ocaml :tangle (eval mli)
val to_int : t -> int
    #+end_src

    Returns the $l_{max}$ value of the shell:

    #+begin_example
Angular_momentum.(to_char D) -> 2
    #+end_example

    #+begin_src ocaml :tangle (eval ml) :exports none
let to_int = function
  | S -> 0         | P -> 1
  | D -> 2         | F -> 3
  | G -> 4         | H -> 5
  | I -> 6         | J -> 7
  | K -> 8         | L -> 9
  | M -> 10        | N -> 11
  | O -> 12        | Int i -> i
    #+end_src

*** ~of_int~

    #+begin_src ocaml :tangle (eval mli)
val of_int : int -> t
    #+end_src

    Returns a shell given an $l$ value.

    #+begin_example
Angular_momentum.of_int 3 -> Angular_momentum.F
    #+end_example

    #+begin_src ocaml :tangle (eval ml) :exports none
let of_int = function
  | 0 -> S         | 1 -> P
  | 2 -> D         | 3 -> F
  | 4 -> G         | 5 -> H
  | 6 -> I         | 7 -> J
  | 8 -> K         | 9 -> L
  | 10 -> M        | 11 -> N
  | 12 -> O        | i -> Int i
    #+end_src

** Shell functions


*** ~n_functions~

    #+begin_src ocaml :tangle (eval mli)
val n_functions : t -> int
    #+end_src

    Returns the number of cartesian functions in a shell.

    #+begin_example
Angular_momentum.n_functions D -> 6
    #+end_example

    #+begin_src ocaml :tangle (eval ml) :exports none
let n_functions a =
  let a =
    to_int a
  in
  (a*a + 3*a + 2)/2
    #+end_src


*** ~zkey_array~

    #+begin_src ocaml :tangle (eval mli)
val zkey_array : kind -> Zkey.t array
    #+end_src

    Array of ~Zkey.t~, where each element is a a key associated with the
    the powers of $x,y,z$.

    #+begin_example
  Angular_momentum.( zkey_array Doublet (P,S) ) ->
    [| {Zkey.left = 0; right = 1125899906842624} ;
       {Zkey.left = 0; right = 1099511627776}    ;
       {Zkey.left = 0; right = 1073741824}       |]
  =

  let s,x,y,z =
      Powers.( of_int_tuple (0,0,0),
               of_int_tuple (1,0,0),
               of_int_tuple (0,1,0),
               of_int_tuple (0,0,1) )
  in
  Array.map (fun (a,b) -> {!Zkey.of_powers_six} a b)
      [| (x,s) ; (y,s) ; (z,s) |]
    #+end_example

    #+begin_src ocaml :tangle (eval ml) :exports none
let zkey_array_memo : (kind, Zkey.t array) Hashtbl.t =
  Hashtbl.create 13

let zkey_array a =

  let keys_1d l =
    let create_z { x ; y ; _ } =
      Powers.of_int_tuple (x,y,l-(x+y))
    in
    let rec create_y accu xyz =
      let { x ; y ; z ;_ } = xyz in
      match y with
      | 0 -> (create_z xyz)::accu
      | _ -> let ynew = y-1 in
          (create_y [@tailcall]) ( (create_z xyz)::accu) (Powers.of_int_tuple (x,ynew,z))
    in
    let rec create_x accu xyz =
      let { x ; z ;_ } = xyz in
      match x with
      | 0 -> (create_y [] xyz)@accu
      | _ -> let xnew = x-1 in
          let ynew = l-xnew in
          (create_x [@tailcall]) ((create_y [] xyz)@accu) (Powers.of_int_tuple (xnew, ynew, z))
    in
    create_x [] (Powers.of_int_tuple (l,0,0))
    |> List.rev
  in

  try
    Hashtbl.find zkey_array_memo a

  with Not_found ->

    let result =
      begin
        match a with
        | Singlet l1 ->
            List.rev_map (fun x -> Zkey.of_powers_three x) (keys_1d @@ to_int l1)

        | Doublet (l1, l2) ->
            List.rev_map (fun a ->
              List.rev_map (fun b -> Zkey.of_powers_six a b) (keys_1d @@ to_int l2)
            ) (keys_1d @@ to_int l1)
            |> List.concat

        | Triplet (l1, l2, l3) ->

            List.rev_map (fun a ->
              List.rev_map (fun b ->
                List.rev_map (fun c ->
                  Zkey.of_powers_nine a b c) (keys_1d @@ to_int l3)
              ) (keys_1d @@ to_int l2)
              |> List.concat
            ) (keys_1d @@ to_int l1)
            |> List.concat

        | Quartet (l1, l2, l3, l4) ->

            List.rev_map (fun a ->
              List.rev_map (fun b ->
                List.rev_map (fun c ->
                  List.rev_map (fun d ->
                    Zkey.of_powers_twelve a b c d) (keys_1d @@ to_int l4)
                ) (keys_1d @@ to_int l3)
                |> List.concat
              ) (keys_1d @@ to_int l2)
              |> List.concat
            ) (keys_1d @@ to_int l1)
            |> List.concat
      end
      |> List.rev
      |> Array.of_list
    in
    Hashtbl.add zkey_array_memo a result;
    result
    #+end_src

** Arithmetic

   #+begin_src ocaml :tangle (eval mli)
val ( + ) : t -> t -> t
val ( - ) : t -> t -> t
   #+end_src

   #+begin_src ocaml :tangle (eval ml) :exports none
let ( + ) a b =
  of_int ( (to_int a) + (to_int b) )

let ( - ) a b =
  of_int ( (to_int a) - (to_int b) )
   #+end_src

** Printers

   Printers can print as a string (~pp_string~) or as an integer (~pp_int~).

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

   #+begin_src ocaml :tangle (eval ml) :exports none
let pp_string ppf x =
  Format.fprintf ppf "@[%s@]" (to_string x)

let pp_int ppf x =
  Format.fprintf ppf "@[%d@]" (to_int x)
   #+end_src


** TODO Tests