From c4b022aeeecf4063f7daabff02f17fca7442e88e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Jan 2024 10:30:24 +0100 Subject: [PATCH] Remove org-mode --- ao/lib/dune | 1 + common/charge.org | 96 --- common/command_line.org | 354 --------- common/constants.org | 79 -- common/coordinate.org | 218 ------ common/lib/charge.ml | 23 +- common/lib/charge.mli | 27 +- common/lib/command_line.ml | 71 +- common/lib/command_line.mli | 98 ++- common/lib/constants.ml | 29 +- common/lib/constants.mli | 36 +- common/lib/coordinate.ml | 80 +-- common/lib/coordinate.mli | 91 ++- common/lib/util.ml | 152 +--- common/lib/util.mli | 151 ++-- common/util.org | 672 ------------------ gaussian/lib/contracted_shell.ml | 18 +- gaussian/lib/contracted_shell_pair.ml | 39 +- gaussian/lib/contracted_shell_pair_couple.ml | 14 +- gaussian/lib/primitive_shell_pair.ml | 244 +++---- gaussian/lib/primitive_shell_pair_couple.ml | 17 +- gaussian_integrals/lib/eri.ml | 30 +- .../lib/two_electron_rr_vectorized.ml | 331 ++++----- gaussian_integrals/lib/zero_m_parameters.ml | 6 +- 24 files changed, 710 insertions(+), 2167 deletions(-) delete mode 100644 common/charge.org delete mode 100644 common/command_line.org delete mode 100644 common/constants.org delete mode 100644 common/coordinate.org delete mode 100644 common/util.org diff --git a/ao/lib/dune b/ao/lib/dune index 9397380..98a1304 100644 --- a/ao/lib/dune +++ b/ao/lib/dune @@ -12,6 +12,7 @@ qcaml.gaussian_integrals qcaml.operators ) + (instrumentation (backend landmarks)) (modules_without_implementation ao_dim) ) diff --git a/common/charge.org b/common/charge.org deleted file mode 100644 index cef7097..0000000 --- a/common/charge.org +++ /dev/null @@ -1,96 +0,0 @@ -#+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 - - -* Charge - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - -** Type - - <<<~Charge.t~>>> - #+begin_src ocaml :tangle (eval mli) -type t - #+end_src - - This type should be used for all charges in the program (electrons, nuclei,...). - - #+begin_src ocaml :tangle (eval ml) :exports none -type t = float - #+end_src - -** Conversions - - #+begin_src ocaml :tangle (eval mli) -val of_float : float -> t -val to_float : t -> float - -val of_int : int -> t -val to_int : t -> int - -val of_string: string -> t -val to_string: t -> string - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -external of_float : float -> t = "%identity" -external to_float : t -> float = "%identity" - -let of_int = float_of_int -let to_int = int_of_float - -let of_string = float_of_string - -let to_string x = - if x > 0. then - Printf.sprintf "+%f" x - else if x < 0. then - Printf.sprintf "%f" x - else - "0.0" - #+end_src - -** Simple operations - - #+begin_src ocaml :tangle (eval mli) -val ( + ) : t -> t -> t -val ( - ) : t -> t -> t -val ( * ) : t -> float -> t -val ( / ) : t -> float -> t -val is_null : t -> bool - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -let gen_op op = - fun a b -> - op (to_float a) (to_float b) - |> of_float - -let ( + ) = gen_op ( +. ) -let ( - ) = gen_op ( -. ) -let ( * ) = gen_op ( *. ) -let ( / ) = gen_op ( /. ) - -let is_null t = t == 0. - #+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 x = - Format.fprintf ppf "@[%s@]" (to_string x) - #+end_src - diff --git a/common/command_line.org b/common/command_line.org deleted file mode 100644 index 17ac72e..0000000 --- a/common/command_line.org +++ /dev/null @@ -1,354 +0,0 @@ -#+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 - -* Command line - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - This module is a wrapper around the ~Getopt~ library and helps to - define command-line arguments. - - Here is an example of how to use this module. - First, define the specification: - #+begin_src ocaml :tangle no -let open Command_line in -begin - set_header_doc (Sys.argv.(0) ^ " - One-line description"); - set_description_doc "Long description of the command."; - set_specs - [ { short='c'; long="check"; opt=Optional; - doc="Checks the input data"; - arg=Without_arg; }; - - { short='b' ; long="basis" ; opt=Mandatory; - arg=With_arg ""; - doc="Name of the file containing the basis set"; } ; - - { short='m' ; long="multiplicity" ; opt=Optional; - arg=With_arg ""; - doc="Spin multiplicity (2S+1). Default is singlet"; } ; - ] -end; - #+end_src - - Then, define what to do with the arguments: - #+begin_src ocaml :tangle no -let c = - Command_line.get_bool "check" -in - -let basis = - match Command_line.get "basis" with - | Some x -> x - | None -> assert false -in - -let multiplicity = - match Command_line.get "multiplicity" with - | None -> 1 - | Some n -> int_of_string n -in - #+end_src - -** Types - - #+NAME:type - #+begin_src ocaml :tangle (eval mli) -type short_opt = char -type long_opt = string -type optional = Mandatory | Optional -type documentation = string -type argument = | With_arg of string - | Without_arg - | With_opt_arg of string - -type description = { - short: short_opt ; - long : long_opt ; - opt : optional ; - doc : documentation ; - arg : argument ; -} - - #+end_src - - - <<>>: in the command line, a dash with a single character - (ex: =ls -l=) - - <<>>: in the command line, two dashes with a word - (ex: =ls --directory=) - - Command-line options can be ~Mandatory~ or ~Optional~ - - Documentation of the option is used in the help function - - Some options require an argument (~ls --ignore="*.ml"~ ), some - don't (~ls -l~) and for some arguments the argument is optional - (~git --log[=]~) - - #+begin_src ocaml :tangle (eval ml) :exports none -<> - #+end_src - -** Mutable attributes - - All the options are stored in the hash table ~dict~ where the key - is the long option and the value is a value of type ~description~. - - #+begin_src ocaml :tangle (eval ml) :exports none -let header_doc = ref "" -let description_doc = ref "" -let footer_doc = ref "" -let anon_args_ref = ref [] -let specs = ref [] -let dict = Hashtbl.create 67 - #+end_src - - #+begin_src ocaml :tangle (eval mli) -val set_header_doc : string -> unit -val set_description_doc : string -> unit -val set_footer_doc : string -> unit - #+end_src - - Functions to set the header, footer and main description of the - documentation provided by the ~help~ function: - - #+begin_src ocaml :tangle (eval ml) :exports none -let set_header_doc s = header_doc := s -let set_description_doc s = description_doc := s -let set_footer_doc s = footer_doc := s - #+end_src - - - #+begin_src ocaml :tangle (eval mli) -val anonymous : long_opt -> optional -> documentation -> description - #+end_src - - Function to create an anonymous argument. - - #+begin_src ocaml :tangle (eval ml) :exports none -let anonymous name opt doc = - { short=' ' ; long=name; opt; doc; arg=Without_arg; } - #+end_src - -** Text formatting functions :noexport: - - Function to print some text such that it fits on the screen - - #+begin_src ocaml :tangle (eval ml) :exports none -let output_text t = - Format.printf "@["; - begin - match Str.split (Str.regexp "\n") t with - | x :: [] -> - Format.printf "@["; - Str.split (Str.regexp " ") x - |> List.iter (fun y -> Format.printf "@[%s@]@ " y) ; - Format.printf "@]" - | t -> - List.iter (fun x -> - Format.printf "@["; - Str.split (Str.regexp " ") x - |> List.iter (fun y -> Format.printf "@[%s@]@ " y) ; - Format.printf "@]@;" - ) t - end; - Format.printf "@]" - #+end_src - - - Function to build the short description of the command-line - arguments, such as - - Example: - #+begin_example -my_program -b [-h] [-u ] -x [--] - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -let output_short x = - match x.short, x.opt, x.arg with - | ' ', Mandatory, _ -> Format.printf "@[%s@]" x.long - | ' ', Optional , _ -> Format.printf "@[[%s]@]" x.long - | _ , Mandatory, Without_arg -> Format.printf "@[-%c@]" x.short - | _ , Optional , Without_arg -> Format.printf "@[[-%c]@]" x.short - | _ , Mandatory, With_arg arg -> Format.printf "@[-%c %s@]" x.short arg - | _ , Optional , With_arg arg -> Format.printf "@[[-%c %s]@]" x.short arg - | _ , Mandatory, With_opt_arg arg -> Format.printf "@[-%c [%s]@]" x.short arg - | _ , Optional , With_opt_arg arg -> Format.printf "@[[-%c [%s]]@]" x.short arg - #+end_src - - - Function to build the long description of the command-line - arguments, such as - - Example: - #+begin_example - -x --xyz= Name of the file containing the nuclear - coordinates in xyz format - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -let output_long max_width x = - let arg = - match x.short, x.arg with - | ' ' , _ -> x.long - | _ , Without_arg -> x.long - | _ , With_arg arg -> Printf.sprintf "%s=%s" x.long arg - | _ , With_opt_arg arg -> Printf.sprintf "%s[=%s]" x.long arg - in - let long = - let l = String.length arg in - arg^(String.make (max_width-l) ' ') - in - Format.printf "@["; - begin - match x.short with - | ' ' -> Format.printf "@[%s @]" long - | short -> Format.printf "@[-%c --%s @]" short long - end; - Format.printf "@]"; - output_text x.doc - #+end_src - -** Query functions - - #+begin_src ocaml :tangle (eval mli) -val get : long_opt -> string option -val get_bool : long_opt -> bool -val anon_args : unit -> string list - #+end_src - - | ~get~ | Returns the argument associated with a long option | - | ~get_bool~ | True if the ~Optional~ argument is present in the command-line | - | ~anon_args~ | Returns the list of anonymous arguments | - - - #+begin_src ocaml :tangle (eval ml) :exports none -let anon_args () = !anon_args_ref - -let help () = - - (* Print the header *) - output_text !header_doc; - Format.printf "@.@."; - - (* Find the anonymous arguments *) - let anon = - List.filter (fun x -> x.short = ' ') !specs - in - - (* Find the options *) - let options = - List.filter (fun x -> x.short <> ' ') !specs - |> List.sort (fun x y -> Char.compare x.short y.short) - in - - (* Find column lengths *) - let max_width = - List.map (fun x -> - ( match x.arg with - | Without_arg -> String.length x.long - | With_arg arg -> String.length x.long + String.length arg - | With_opt_arg arg -> String.length x.long + String.length arg + 2 - ) - + ( if x.opt = Optional then 2 else 0) - ) !specs - |> List.fold_left max 0 - in - - - (* Print usage *) - Format.printf "@[@[Usage:@,@,@[@[%s@]" Sys.argv.(0); - List.iter (fun x -> Format.printf "@ "; output_short x) options; - Format.printf "@ @[[--]@]"; - List.iter (fun x -> Format.printf "@ "; output_short x;) anon; - Format.printf "@]@,@]@,"; - - - (* Print arguments and doc *) - Format.printf "@[Arguments:@,"; - Format.printf "@[" ; - List.iter (fun x -> Format.printf "@ "; output_long max_width x) anon; - Format.printf "@]@,@]@,"; - - - (* Print options and doc *) - Format.printf "@[Options:@,"; - - Format.printf "@[" ; - List.iter (fun x -> Format.printf "@ "; output_long max_width x) options; - Format.printf "@]@,@]@,"; - - - (* Print footer *) - if !description_doc <> "" then - begin - Format.printf "@[Description:@,@,"; - output_text !description_doc; - Format.printf "@," - end; - - (* Print footer *) - output_text !footer_doc; - Format.printf "@." - - -let get x = - try Some (Hashtbl.find dict x) - with Not_found -> None - - -let get_bool x = Hashtbl.mem dict x - #+end_src - -** Specification - - #+begin_src ocaml :tangle (eval mli) -val set_specs : description list -> unit - #+end_src - - Sets the specifications of the current program from a list of - ~descrption~ variables. - - #+begin_src ocaml :tangle (eval ml) :exports none -let set_specs specs_in = - specs := { short = 'h' ; - long = "help" ; - doc = "Prints the help message." ; - arg = Without_arg ; - opt = Optional ; - } :: specs_in; - - let cmd_specs = - List.filter (fun x -> x.short != ' ') !specs - |> List.map (fun { short ; long ; arg ; _ } -> - match arg with - | With_arg _ -> - (short, long, None, Some (fun x -> Hashtbl.replace dict long x) ) - | Without_arg -> - (short, long, Some (fun () -> Hashtbl.replace dict long ""), None) - | With_opt_arg _ -> - (short, long, Some (fun () -> Hashtbl.replace dict long ""), - Some (fun x -> Hashtbl.replace dict long x) ) - ) - in - - Getopt.parse_cmdline cmd_specs (fun x -> anon_args_ref := !anon_args_ref @ [x]); - - if (get_bool "help") then - (help () ; exit 0); - - (* Check that all mandatory arguments are set *) - List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs - |> List.iter (fun x -> - match get x.long with - | Some _ -> () - | None -> failwith ("Error: --"^x.long^" option is missing.") - ) - #+end_src diff --git a/common/constants.org b/common/constants.org deleted file mode 100644 index c5aa782..0000000 --- a/common/constants.org +++ /dev/null @@ -1,79 +0,0 @@ -#+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 - -* Constants - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - All constants used in the program. - -** Thresholds - - #+begin_src ocaml :tangle (eval mli) -val epsilon : float -val integrals_cutoff : float - #+end_src - - | ~epsilon~ | Value below which a float is considered null. Default is \epsilon = 2.10^{-15} | - | ~integrals_cutoff~ | Cutoff value for integrals. Default is \epsilon | - - #+begin_src ocaml :tangle (eval ml) :exports none -let epsilon = 2.e-15 -let integrals_cutoff = epsilon - #+end_src - -** Mathematical constants - - #+begin_src ocaml :tangle (eval mli) -val pi : float -val two_pi : float -val sq_pi : float -val sq_pi_over_two : float -val pi_inv : float -val two_over_sq_pi : float - #+end_src - - | ~pi~ | $\pi = 3.141~592~653~589~793~12$ | - | ~two_pi~ | $2 \pi$ | - | ~sq_pi~ | $\sqrt{\pi}$ | - | ~sq_pi_over_two~ | $\sqrt{\pi} / 2$ | - | ~pi_inv~ | $1 / \pi$ | - | ~two_over_sq_pi~ | $2 / \sqrt{\pi}$ | - - #+begin_src ocaml :tangle (eval ml) :exports none -let pi = acos (-1.) -let two_pi = 2. *. pi -let sq_pi = sqrt pi -let sq_pi_over_two = sq_pi *. 0.5 -let pi_inv = 1. /. pi -let two_over_sq_pi = 2. /. sq_pi - #+end_src - -** Physical constants - - #+begin_src ocaml :tangle (eval mli) -val a0 : float -val a0_inv : float -val ha_to_ev : float -val ev_to_ha : float - #+end_src - - | ~a0~ | Bohr's radius : $a_0 = 0.529~177~210~67(23)$ angstrom | - | ~a0_inv~ | $1 / a_0$ | - | ~ha_to_ev~ | Hartree to eV conversion factor : $27.211~386~02(17)$ | - | ~ev_to_ha~ | eV to Hartree conversion factor : 1 / ~ha_to_ev~ | - - #+begin_src ocaml :tangle (eval ml) :exports none -let a0 = 0.529_177_210_67 -let a0_inv = 1. /. a0 -let ha_to_ev = 27.211_386_02 -let ev_to_ha = 1. /. ha_to_ev - #+end_src diff --git a/common/coordinate.org b/common/coordinate.org deleted file mode 100644 index 0ded6d7..0000000 --- a/common/coordinate.org +++ /dev/null @@ -1,218 +0,0 @@ -#+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 - -* Coordinate - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - Coordinates in 3D space. - - All operations on points are done in atomic units. Therefore, - all coordinates are given in bohr and manipulated with this - module. - -** Type - - <<<~Coordinate.t~>>> - #+NAME: types - #+begin_src ocaml :tangle (eval mli) -type bohr -type angstrom - -type xyz = { - x : float ; - y : float ; - z : float ; -} - -type 'a point = xyz - -type t = bohr point - -type axis = X | Y | Z - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -type bohr -type angstrom - -type xyz = { - x : float ; - y : float ; - z : float ; -} - -type 'a point = xyz - -type t = bohr point - -type axis = X | Y | Z - #+end_src - -** Creation - - #+begin_src ocaml :tangle (eval mli) -val make : 'a point -> t -val make_angstrom : 'a point -> angstrom point -val zero : bohr point - #+end_src - - | ~make~ | Creates a point in atomic units | - | ~make_angstrom~ | Creates a point in angstrom | - | ~zero~ | $(0., 0., 0.)$ | - - #+begin_src ocaml :tangle (eval ml) :exports none -external make : 'a point -> t = "%identity" -external make_angstrom : 'a point -> angstrom point = "%identity" - -let zero = - make { x = 0. ; y = 0. ; z = 0. } - #+end_src - -** Conversion - - #+begin_src ocaml :tangle (eval mli) -val bohr_to_angstrom : bohr point -> angstrom point -val angstrom_to_bohr : angstrom point -> bohr point - #+end_src - - | ~bohr_to_angstrom~ | Converts a point in bohr to angstrom | - | ~angstrom_to_bohr~ | Converts a point in angstrom to bohr | - - #+begin_src ocaml :tangle (eval ml) :exports none -let b_to_a b = Constants.a0 *. b -let bohr_to_angstrom { x ; y ; z } = - make { x = b_to_a x ; - y = b_to_a y ; - z = b_to_a z ; } - - -let a_to_b a = Constants.a0_inv *. a -let angstrom_to_bohr { x ; y ; z } = - make { x = a_to_b x ; - y = a_to_b y ; - z = a_to_b z ; } - #+end_src - -** Vector operations - - #+begin_src ocaml :tangle (eval mli) -val neg : t -> t -val get : axis -> bohr point -> float -val dot : t -> t -> float -val norm : t -> float -val ( |. ) : float -> t -> t -val ( |+ ) : t -> t -> t -val ( |- ) : t -> t -> t - #+end_src - - | ~neg~ | Negative of a point | - | ~get~ | Extracts the projection of the coordinate on an axis | - | ~dot~ | Dot product | - | ~norm~ | $\ell{^2}$ norm of the vector | - | $\vert .$ | Scales the vector by a constant | - | $\vert +$ | Adds two vectors | - | $\vert -$ | Subtracts two vectors | - - Example: - #+begin_example -Coordinate.neg { x=1. ; y=2. ; z=3. } ;; -- : Coordinate.t = -1.0000 -2.0000 -3.0000 - -Coordinate.(get Y { x=1. ; y=2. ; z=3. }) ;; -- : float = 2. - -Coordinate.( - 2. |. { x=1. ; y=2. ; z=3. } -) ;; -- : Coordinate.t = 2.0000 4.0000 6.0000 - -Coordinate.( - { x=1. ; y=2. ; z=3. } |+ { x=2. ; y=3. ; z=1. } -);; -- : Coordinate.t = 3.0000 5.0000 4.0000 - -Coordinate.( - { x=1. ; y=2. ; z=3. } |- { x=2. ; y=3. ; z=1. } -);; -- : Coordinate.t = -1.0000 -1.0000 2.0000 - #+end_example - - - #+begin_src ocaml :tangle (eval ml) :exports none -let get axis { x ; y ; z } = - match axis with - | X -> x - | Y -> y - | Z -> z - - -let ( |. ) s { x ; y ; z } = - make { x = s *. x ; - y = s *. y ; - z = s *. z ; } - - -let ( |+ ) - { x = x1 ; y = y1 ; z = z1 } - { x = x2 ; y = y2 ; z = z2 } = - make { x = x1 +. x2 ; - y = y1 +. y2 ; - z = z1 +. z2 ; } - - -let ( |- ) - { x = x1 ; y = y1 ; z = z1 } - { x = x2 ; y = y2 ; z = z2 } = - make { x = x1 -. x2 ; - y = y1 -. y2 ; - z = z1 -. z2 ; } - - -let neg a = -1. |. a - - -let dot - { x = x1 ; y = y1 ; z = z1 } - { x = x2 ; y = y2 ; z = z2 } = - x1 *. x2 +. - y1 *. y2 +. - z1 *. z2 - - -let norm u = - sqrt ( dot u u ) - #+end_src - -** Printers - - #+begin_src ocaml :tangle (eval mli) -val pp : Format.formatter -> t -> unit -val pp_bohr: Format.formatter -> t -> unit -val pp_angstrom : Format.formatter -> t -> unit - #+end_src - - Coordinates can be printed in bohr or angstrom. - - #+begin_src ocaml :tangle (eval ml) :exports none -open Format -let pp ppf c = - fprintf ppf "@[@[%8.4f@] @[%8.4f@] @[%8.4f@]@]" c.x c.y c.z - -let pp_bohr ppf c = - fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Bohr)@]" c.x c.y c.z - -let pp_angstrom ppf c = - let c = bohr_to_angstrom c in - fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Angs)@]" c.x c.y c.z - #+end_src - diff --git a/common/lib/charge.ml b/common/lib/charge.ml index fb5cc13..1df38d3 100644 --- a/common/lib/charge.ml +++ b/common/lib/charge.ml @@ -1,31 +1,24 @@ - - (* This type should be used for all charges in the program (electrons, nuclei,...). *) - - -(* [[file:~/QCaml/common/charge.org::*Type][Type:2]] *) type t = float -(* Type:2 ends here *) -(* [[file:~/QCaml/common/charge.org::*Conversions][Conversions:2]] *) +(** Conversions *) external of_float : float -> t = "%identity" external to_float : t -> float = "%identity" let of_int = float_of_int let to_int = int_of_float -let of_string = float_of_string +let of_string = float_of_string -let to_string x = +let to_string x = if x > 0. then Printf.sprintf "+%f" x else if x < 0. then Printf.sprintf "%f" x - else + else "0.0" -(* Conversions:2 ends here *) -(* [[file:~/QCaml/common/charge.org::*Simple operations][Simple operations:2]] *) +(** Simple operations *) let gen_op op = fun a b -> op (to_float a) (to_float b) @@ -37,9 +30,7 @@ let ( * ) = gen_op ( *. ) let ( / ) = gen_op ( /. ) let is_null t = t == 0. -(* Simple operations:2 ends here *) -(* [[file:~/QCaml/common/charge.org::*Printers][Printers:2]] *) -let pp ppf x = +(** Printers *) +let pp ppf x = Format.fprintf ppf "@[%s@]" (to_string x) -(* Printers:2 ends here *) diff --git a/common/lib/charge.mli b/common/lib/charge.mli index 7b316a8..becad1f 100644 --- a/common/lib/charge.mli +++ b/common/lib/charge.mli @@ -1,15 +1,10 @@ -(* Type - * - * <<<~Charge.t~>>> *) +(** Type *) -(* [[file:~/QCaml/common/charge.org::*Type][Type:1]] *) type t -(* Type:1 ends here *) - -(* Conversions *) -(* [[file:~/QCaml/common/charge.org::*Conversions][Conversions:1]] *) +(** Conversions *) + val of_float : float -> t val to_float : t -> float @@ -18,22 +13,18 @@ val to_int : t -> int val of_string: string -> t val to_string: t -> string -(* Conversions:1 ends here *) - -(* Simple operations *) -(* [[file:~/QCaml/common/charge.org::*Simple operations][Simple operations:1]] *) +(** Simple operations *) + val ( + ) : t -> t -> t val ( - ) : t -> t -> t val ( * ) : t -> float -> t -val ( / ) : t -> float -> t +val ( / ) : t -> float -> t val is_null : t -> bool -(* Simple operations:1 ends here *) - -(* Printers *) -(* [[file:~/QCaml/common/charge.org::*Printers][Printers:1]] *) +(** Printers *) + val pp : Format.formatter -> t -> unit -(* Printers:1 ends here *) + diff --git a/common/lib/command_line.ml b/common/lib/command_line.ml index f1081a1..2461b86 100644 --- a/common/lib/command_line.ml +++ b/common/lib/command_line.ml @@ -1,17 +1,4 @@ - - -(* - <<>>: in the command line, a dash with a single character - * (ex: =ls -l=) - * - <<>>: in the command line, two dashes with a word - * (ex: =ls --directory=) - * - Command-line options can be ~Mandatory~ or ~Optional~ - * - Documentation of the option is used in the help function - * - Some options require an argument (~ls --ignore="*.ml"~ ), some - * don't (~ls -l~) and for some arguments the argument is optional - * (~git --log[=]~) *) - - -(* [[file:~/QCaml/common/command_line.org::*Types][Types:2]] *) +(** Types *) type short_opt = char type long_opt = string type optional = Mandatory | Optional @@ -27,51 +14,34 @@ type description = { doc : documentation ; arg : argument ; } -(* Types:2 ends here *) -(* Mutable attributes - * +(** Mutable attributes + * * All the options are stored in the hash table ~dict~ where the key * is the long option and the value is a value of type ~description~. *) -(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:1]] *) let header_doc = ref "" let description_doc = ref "" let footer_doc = ref "" let anon_args_ref = ref [] let specs = ref [] let dict = Hashtbl.create 67 -(* Mutable attributes:1 ends here *) - -(* Functions to set the header, footer and main description of the - * documentation provided by the ~help~ function: *) - - -(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:3]] *) let set_header_doc s = header_doc := s let set_description_doc s = description_doc := s let set_footer_doc s = footer_doc := s -(* Mutable attributes:3 ends here *) - -(* Function to create an anonymous argument. *) - - -(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:5]] *) let anonymous name opt doc = { short=' ' ; long=name; opt; doc; arg=Without_arg; } -(* Mutable attributes:5 ends here *) -(* Text formatting functions :noexport: + +(* Text formatting functions * * Function to print some text such that it fits on the screen *) - -(* [[file:~/QCaml/common/command_line.org::*Text formatting functions][Text formatting functions:1]] *) let output_text t = Format.printf "@["; begin @@ -90,8 +60,6 @@ let output_text t = ) t end; Format.printf "@]" -(* Text formatting functions:1 ends here *) - @@ -99,12 +67,11 @@ let output_text t = * arguments, such as * * Example: - * #+begin_example + * my_program -b [-h] [-u ] -x [--] - * #+end_example *) + *) -(* [[file:~/QCaml/common/command_line.org::*Text formatting functions][Text formatting functions:2]] *) let output_short x = match x.short, x.opt, x.arg with | ' ', Mandatory, _ -> Format.printf "@[%s@]" x.long @@ -115,8 +82,6 @@ let output_short x = | _ , Optional , With_arg arg -> Format.printf "@[[-%c %s]@]" x.short arg | _ , Mandatory, With_opt_arg arg -> Format.printf "@[-%c [%s]@]" x.short arg | _ , Optional , With_opt_arg arg -> Format.printf "@[[-%c [%s]]@]" x.short arg -(* Text formatting functions:2 ends here *) - @@ -124,13 +89,12 @@ let output_short x = * arguments, such as * * Example: - * #+begin_example + * * -x --xyz= Name of the file containing the nuclear * coordinates in xyz format - * #+end_example *) + *) -(* [[file:~/QCaml/common/command_line.org::*Text formatting functions][Text formatting functions:3]] *) let output_long max_width x = let arg = match x.short, x.arg with @@ -151,17 +115,9 @@ let output_long max_width x = end; Format.printf "@]"; output_text x.doc -(* Text formatting functions:3 ends here *) -(* | ~get~ | Returns the argument associated with a long option | - * | ~get_bool~ | True if the ~Optional~ argument is present in the command-line | - * | ~anon_args~ | Returns the list of anonymous arguments | *) - - - -(* [[file:~/QCaml/common/command_line.org::*Query functions][Query functions:2]] *) let anon_args () = !anon_args_ref let help () = @@ -237,15 +193,8 @@ let get x = let get_bool x = Hashtbl.mem dict x -(* Query functions:2 ends here *) - -(* Sets the specifications of the current program from a list of - * ~descrption~ variables. *) - - -(* [[file:~/QCaml/common/command_line.org::*Specification][Specification:2]] *) let set_specs specs_in = specs := { short = 'h' ; long = "help" ; @@ -280,4 +229,4 @@ let set_specs specs_in = | Some _ -> () | None -> failwith ("Error: --"^x.long^" option is missing.") ) -(* Specification:2 ends here *) + diff --git a/common/lib/command_line.mli b/common/lib/command_line.mli index 12e9a0e..8a18fce 100644 --- a/common/lib/command_line.mli +++ b/common/lib/command_line.mli @@ -1,8 +1,64 @@ -(* Types - * - * #+NAME:type *) +(** Command line *) + +(* This module is a wrapper around the ~Getopt~ library and helps to + * define command-line arguments. + * + * Here is an example of how to use this module. + * First, define the specification: + * + * let open Command_line in + * begin + * set_header_doc (Sys.argv.(0) ^ " - One-line description"); + * set_description_doc "Long description of the command."; + * set_specs + * [ { short='c'; long="check"; opt=Optional; + * doc="Checks the input data"; + * arg=Without_arg; }; + * + * { short='b' ; long="basis" ; opt=Mandatory; + * arg=With_arg ""; + * doc="Name of the file containing the basis set"; } ; + * + * { short='m' ; long="multiplicity" ; opt=Optional; + * arg=With_arg ""; + * doc="Spin multiplicity (2S+1). Default is singlet"; } ; + * ] + * end; + * + * Then, define what to do with the arguments: + * + * let c = + * Command_line.get_bool "check" + * in + * + * let basis = + * match Command_line.get "basis" with + * | Some x -> x + * | None -> assert false + * in + * + * let multiplicity = + * match Command_line.get "multiplicity" with + * | None -> 1 + * | Some n -> int_of_string n + * in + * ... + *) + + +(* - <<>>: in the command line, a dash with a single character + * (ex: =ls -l=) + * - <<>>: in the command line, two dashes with a word + * (ex: =ls --directory=) + * - Command-line options can be `Mandatory` or `Optional` + * - Documentation of the option is used in the help function + * - Some options require an argument (`ls --ignore="*.ml"`), some + * don't (`ls -l`) and for some arguments the argument is optional + * (`git --log[=]`) *) + + +(** Types *) -(* [[file:~/QCaml/common/command_line.org::type][type]] *) type short_opt = char type long_opt = string type optional = Mandatory | Optional @@ -18,30 +74,38 @@ type description = { doc : documentation ; arg : argument ; } -(* type ends here *) -(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:2]] *) + +(** Mutable attributes *) + val set_header_doc : string -> unit +(** Sets the header of the documentation provided by the ~help~ function: *) + val set_description_doc : string -> unit +(** Sets the description of the documentation provided by the ~help~ function: *) + val set_footer_doc : string -> unit -(* Mutable attributes:2 ends here *) +(** Sets the footer of the documentation provided by the ~help~ function: *) -(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:4]] *) val anonymous : long_opt -> optional -> documentation -> description -(* Mutable attributes:4 ends here *) - -(* Query functions *) +(** Creates an anonymous argument. *) -(* [[file:~/QCaml/common/command_line.org::*Query functions][Query functions:1]] *) +(** Query functions *) + val get : long_opt -> string option +(** Returns the argument associated with a long option *) + val get_bool : long_opt -> bool +(** True if the ~Optional~ argument is present in the command-line *) + val anon_args : unit -> string list -(* Query functions:1 ends here *) - -(* Specification *) +(** Returns the list of anonymous arguments *) -(* [[file:~/QCaml/common/command_line.org::*Specification][Specification:1]] *) +(** Specification *) + val set_specs : description list -> unit -(* Specification:1 ends here *) +(** Sets the specifications of the current program from a list of + * ~descrption~ variables. *) + diff --git a/common/lib/constants.ml b/common/lib/constants.ml index db2233c..beeafd2 100644 --- a/common/lib/constants.ml +++ b/common/lib/constants.ml @@ -1,44 +1,21 @@ - - -(* | ~epsilon~ | Value below which a float is considered null. Default is \epsilon = 2.10^{-15} | - * | ~integrals_cutoff~ | Cutoff value for integrals. Default is \epsilon | *) - - -(* [[file:~/QCaml/common/constants.org::*Thresholds][Thresholds:2]] *) +(** Thresholds *) let epsilon = 2.e-15 let integrals_cutoff = epsilon -(* Thresholds:2 ends here *) +(** Mathematical constants *) -(* | ~pi~ | $\pi = 3.141~592~653~589~793~12$ | - * | ~two_pi~ | $2 \pi$ | - * | ~sq_pi~ | $\sqrt{\pi}$ | - * | ~sq_pi_over_two~ | $\sqrt{\pi} / 2$ | - * | ~pi_inv~ | $1 / \pi$ | - * | ~two_over_sq_pi~ | $2 / \sqrt{\pi}$ | *) - - -(* [[file:~/QCaml/common/constants.org::*Mathematical constants][Mathematical constants:2]] *) let pi = acos (-1.) let two_pi = 2. *. pi let sq_pi = sqrt pi let sq_pi_over_two = sq_pi *. 0.5 let pi_inv = 1. /. pi let two_over_sq_pi = 2. /. sq_pi -(* Mathematical constants:2 ends here *) +(** Physical constants *) -(* | ~a0~ | Bohr's radius : $a_0 = 0.529~177~210~67(23)$ angstrom | - * | ~a0_inv~ | $1 / a_0$ | - * | ~ha_to_ev~ | Hartree to eV conversion factor : $27.211~386~02(17)$ | - * | ~ev_to_ha~ | eV to Hartree conversion factor : 1 / ~ha_to_ev~ | *) - - -(* [[file:~/QCaml/common/constants.org::*Physical constants][Physical constants:2]] *) let a0 = 0.529_177_210_67 let a0_inv = 1. /. a0 let ha_to_ev = 27.211_386_02 let ev_to_ha = 1. /. ha_to_ev -(* Physical constants:2 ends here *) diff --git a/common/lib/constants.mli b/common/lib/constants.mli index 0c85ba6..ed8a6c4 100644 --- a/common/lib/constants.mli +++ b/common/lib/constants.mli @@ -1,29 +1,43 @@ (* Thresholds *) - -(* [[file:~/QCaml/common/constants.org::*Thresholds][Thresholds:1]] *) val epsilon : float +(** Value below which a float is considered null. Default is \epsilon = 2.10^{-15}. *) + val integrals_cutoff : float -(* Thresholds:1 ends here *) - -(* Mathematical constants *) +(** Cutoff value for integrals. Default is \epsilon. *) -(* [[file:~/QCaml/common/constants.org::*Mathematical constants][Mathematical constants:1]] *) +(** Mathematical constants *) val pi : float +(** $\pi = 3.141~592~653~589~793~12$ *) + val two_pi : float +(** $2 \pi$ *) + val sq_pi : float +(** $\sqrt{\pi}$ *) + val sq_pi_over_two : float +(** $\sqrt{\pi} / 2$ *) + val pi_inv : float +(** $1 / \pi$ *) + val two_over_sq_pi : float -(* Mathematical constants:1 ends here *) +(** $2 / \sqrt{\pi}$ *) -(* Physical constants *) - -(* [[file:~/QCaml/common/constants.org::*Physical constants][Physical constants:1]] *) +(** Physical constants *) + val a0 : float +(** Bohr's radius : $a_0 = 0.529~177~210~67(23)$ angstrom. *) + val a0_inv : float +(** $1 / a_0$ *) + val ha_to_ev : float +(** Hartree to eV conversion factor : $27.211~386~02(17)$ *) + val ev_to_ha : float -(* Physical constants:1 ends here *) +(** eV to Hartree conversion factor : 1 / ~ha_to_ev~ *) + diff --git a/common/lib/coordinate.ml b/common/lib/coordinate.ml index cbf11d7..b27ea8a 100644 --- a/common/lib/coordinate.ml +++ b/common/lib/coordinate.ml @@ -1,6 +1,7 @@ -(* [[file:~/QCaml/common/coordinate.org::*Type][Type:2]] *) -type bohr -type angstrom +(** Types *) + +type bohr +type angstrom type xyz = { x : float ; @@ -13,86 +14,42 @@ type 'a point = xyz type t = bohr point type axis = X | Y | Z -(* Type:2 ends here *) -(* | ~make~ | Creates a point in atomic units | - * | ~make_angstrom~ | Creates a point in angstrom | - * | ~zero~ | $(0., 0., 0.)$ | *) -(* [[file:~/QCaml/common/coordinate.org::*Creation][Creation:2]] *) +(** Creation *) external make : 'a point -> t = "%identity" external make_angstrom : 'a point -> angstrom point = "%identity" let zero = make { x = 0. ; y = 0. ; z = 0. } -(* Creation:2 ends here *) -(* | ~bohr_to_angstrom~ | Converts a point in bohr to angstrom | - * | ~angstrom_to_bohr~ | Converts a point in angstrom to bohr | *) +(** Conversion *) - -(* [[file:~/QCaml/common/coordinate.org::*Conversion][Conversion:2]] *) -let b_to_a b = Constants.a0 *. b +let b_to_a b = Constants.a0 *. b let bohr_to_angstrom { x ; y ; z } = make { x = b_to_a x ; y = b_to_a y ; z = b_to_a z ; } -let a_to_b a = Constants.a0_inv *. a -let angstrom_to_bohr { x ; y ; z } = +let a_to_b a = Constants.a0_inv *. a +let angstrom_to_bohr { x ; y ; z } = make { x = a_to_b x ; y = a_to_b y ; z = a_to_b z ; } -(* Conversion:2 ends here *) -(* | ~neg~ | Negative of a point | - * | ~get~ | Extracts the projection of the coordinate on an axis | - * | ~dot~ | Dot product | - * | ~norm~ | $\ell{^2}$ norm of the vector | - * | $\vert .$ | Scales the vector by a constant | - * | $\vert +$ | Adds two vectors | - * | $\vert -$ | Subtracts two vectors | - * - * Example: - * #+begin_example - * Coordinate.neg { x=1. ; y=2. ; z=3. } ;; - * - : Coordinate.t = -1.0000 -2.0000 -3.0000 - * - * Coordinate.(get Y { x=1. ; y=2. ; z=3. }) ;; - * - : float = 2. - * - * Coordinate.( - * 2. |. { x=1. ; y=2. ; z=3. } - * ) ;; - * - : Coordinate.t = 2.0000 4.0000 6.0000 - * - * Coordinate.( - * { x=1. ; y=2. ; z=3. } |+ { x=2. ; y=3. ; z=1. } - * );; - * - : Coordinate.t = 3.0000 5.0000 4.0000 - * - * Coordinate.( - * { x=1. ; y=2. ; z=3. } |- { x=2. ; y=3. ; z=1. } - * );; - * - : Coordinate.t = -1.0000 -1.0000 2.0000 - * #+end_example *) - - - -(* [[file:~/QCaml/common/coordinate.org::*Vector operations][Vector operations:2]] *) -let get axis { x ; y ; z } = +let get axis { x ; y ; z } = match axis with - | X -> x - | Y -> y - | Z -> z + | X -> x + | Y -> y + | Z -> z let ( |. ) s { x ; y ; z } = @@ -130,22 +87,19 @@ let dot let norm u = sqrt ( dot u u ) -(* Vector operations:2 ends here *) -(* Coordinates can be printed in bohr or angstrom. *) +(** Printers *) -(* [[file:~/QCaml/common/coordinate.org::*Printers][Printers:2]] *) open Format -let pp ppf c = +let pp ppf c = fprintf ppf "@[@[%8.4f@] @[%8.4f@] @[%8.4f@]@]" c.x c.y c.z -let pp_bohr ppf c = +let pp_bohr ppf c = fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Bohr)@]" c.x c.y c.z -let pp_angstrom ppf c = +let pp_angstrom ppf c = let c = bohr_to_angstrom c in fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Angs)@]" c.x c.y c.z -(* Printers:2 ends here *) diff --git a/common/lib/coordinate.mli b/common/lib/coordinate.mli index 9e6de64..e752ff0 100644 --- a/common/lib/coordinate.mli +++ b/common/lib/coordinate.mli @@ -1,11 +1,16 @@ -(* Type - * - * <<<~Coordinate.t~>>> - * #+NAME: types *) +(** Coordinates + * + * Coordinates in 3D space. + * + * All operations on points are done in atomic units. + * Therefore, all coordinates are given in bohr and + * manipulated with this module. + *) -(* [[file:~/QCaml/common/coordinate.org::types][types]] *) -type bohr -type angstrom +(** Types *) + +type bohr +type angstrom type xyz = { x : float ; @@ -18,43 +23,85 @@ type 'a point = xyz type t = bohr point type axis = X | Y | Z -(* types ends here *) - -(* Creation *) -(* [[file:~/QCaml/common/coordinate.org::*Creation][Creation:1]] *) +(** Creation *) + val make : 'a point -> t +(** Creates a point in atomic units *) + val make_angstrom : 'a point -> angstrom point +(** Creates a point in angstrom *) + val zero : bohr point -(* Creation:1 ends here *) - -(* Conversion *) +(** (0., 0., 0.) *) -(* [[file:~/QCaml/common/coordinate.org::*Conversion][Conversion:1]] *) +(** Conversion *) + val bohr_to_angstrom : bohr point -> angstrom point +(** Converts a point in bohr to angstrom *) + val angstrom_to_bohr : angstrom point -> bohr point -(* Conversion:1 ends here *) +(** Converts a point in angstrom to bohr *) (* Vector operations *) -(* [[file:~/QCaml/common/coordinate.org::*Vector operations][Vector operations:1]] *) +(** Vector operations *) + val neg : t -> t +(** Negative of a point *) + val get : axis -> bohr point -> float +(** Extracts the projection of the coordinate on an axis *) + val dot : t -> t -> float +(** Dot product *) + val norm : t -> float +(** $\ell{^2}$ norm of the vector *) + val ( |. ) : float -> t -> t +(** Scales the vector by a constant *) + val ( |+ ) : t -> t -> t +(** Adds two vectors *) + val ( |- ) : t -> t -> t -(* Vector operations:1 ends here *) - -(* Printers *) +(** Subtracts two vectors *) -(* [[file:~/QCaml/common/coordinate.org::*Printers][Printers:1]] *) +(** Example + + * Coordinate.neg { x=1. ; y=2. ; z=3. } ;; + * - : Coordinate.t = -1.0000 -2.0000 -3.0000 + * + * Coordinate.(get Y { x=1. ; y=2. ; z=3. }) ;; + * - : float = 2. + * + * Coordinate.( + * 2. |. { x=1. ; y=2. ; z=3. } + * ) ;; + * - : Coordinate.t = 2.0000 4.0000 6.0000 + * + * Coordinate.( + * { x=1. ; y=2. ; z=3. } |+ { x=2. ; y=3. ; z=1. } + * );; + * - : Coordinate.t = 3.0000 5.0000 4.0000 + * + * Coordinate.( + * { x=1. ; y=2. ; z=3. } |- { x=2. ; y=3. ; z=1. } + * );; + * - : Coordinate.t = -1.0000 -1.0000 2.0000 + *) + + + +(** Printers *) + +(* Coordinates can be printed in bohr or angstrom. *) + val pp : Format.formatter -> t -> unit val pp_bohr: Format.formatter -> t -> unit val pp_angstrom : Format.formatter -> t -> unit -(* Printers:1 ends here *) diff --git a/common/lib/util.ml b/common/lib/util.ml index 2a9bab8..0534011 100644 --- a/common/lib/util.ml +++ b/common/lib/util.ml @@ -1,51 +1,31 @@ -(* [[file:~/QCaml/common/util.org::*Erf][Erf:3]] *) +(** Utility functions *) + external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] -(* Erf:3 ends here *) -(* [[file:~/QCaml/common/util.org::*Erfc][Erfc:3]] *) external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] -(* Erfc:3 ends here *) -(* [[file:~/QCaml/common/util.org::*Gamma][Gamma:3]] *) external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] -(* Gamma:3 ends here *) -(* [[file:~/QCaml/common/util.org::*Popcnt][Popcnt:3]] *) external popcnt : int64 -> int32 = "popcnt_bytecode" "popcnt" [@@unboxed] [@@noalloc] let popcnt i = (popcnt [@inlined] ) i |> Int32.to_int -(* Popcnt:3 ends here *) -(* [[file:~/QCaml/common/util.org::*Trailz][Trailz:3]] *) external trailz : int64 -> int32 = "trailz_bytecode" "trailz" "int" [@@unboxed] [@@noalloc] let trailz i = trailz i |> Int32.to_int -(* Trailz:3 ends here *) -(* [[file:~/QCaml/common/util.org::*Leadz][Leadz:3]] *) external leadz : int64 -> int32 = "leadz_bytecode" "leadz" "int" [@@unboxed] [@@noalloc] let leadz i = leadz i |> Int32.to_int -(* Leadz:3 ends here *) -(* | ~fact~ | Factorial function. | - * | ~binom~ | Binomial coefficient. ~binom n k~ = $C_n^k$ | - * | ~binom_float~ | float variant of ~binom~ | - * | ~pow~ | Fast implementation of the power function for small integer powers | - * | ~chop~ | In ~chop a f~, evaluate ~f~ only if the absolute value of ~a~ is larger than ~Constants.epsilon~, and return ~a *. f ()~. | - * | ~float_of_int_fast~ | Faster implementation of float_of_int for small positive ints | - * | ~not_implemented~ | Fails if some functionality is not implemented | - * | ~of_some~ | Extracts the value of an option | *) - -(* [[file:~/QCaml/common/util.org::*General functions][General functions:2]] *) let memo_float_of_int = Array.init 64 float_of_int @@ -79,7 +59,7 @@ let fact = function let binom = - let memo = + let memo = let m = Array.make_matrix 64 64 0 in for n=0 to Array.length m - 1 do m.(n).(0) <- 1; @@ -95,7 +75,7 @@ let binom = assert (n >= k); if k = 0 || k = n then 1 - else if n < 64 then + else if n < 64 then memo.(n).(k) else f (n-1) (k-1) + f (n-1) k @@ -126,34 +106,19 @@ let chop f g = exception Not_implemented of string -let not_implemented string = +let not_implemented string = raise (Not_implemented string) let of_some = function | Some a -> a | None -> assert false -(* General functions:2 ends here *) -(* The lower [[https://en.wikipedia.org/wiki/Incomplete_gamma_function][Incomplete Gamma function]] is implemented : - * \[ - * \gamma(\alpha,x) = \int_0^x e^{-t} t^{\alpha-1} dt - * \] - * - * p: $\frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ - * - * q: $\frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ - * - * Reference : Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten - * (New Algorithm handbook in C language) (Gijyutsu hyouron sha, - * Tokyo, 1991) p.227 [in Japanese] *) - -(* [[file:~/QCaml/common/util.org::*Functions related to the Boys function][Functions related to the Boys function:2]] *) -let incomplete_gamma ~alpha x = +let incomplete_gamma ~alpha x = assert (alpha >= 0.); assert (x >= 0.); let a = alpha in @@ -161,7 +126,7 @@ let incomplete_gamma ~alpha x = let gf = gamma_float alpha in let loggamma_a = log gf in let rec p_gamma x = - if x >= 1. +. a then 1. -. q_gamma x + if x >= 1. +. a then 1. -. q_gamma x else if x = 0. then 0. else let rec pg_loop prev res term k = @@ -175,7 +140,7 @@ let incomplete_gamma ~alpha x = pg_loop min_float r0 r0 1. and q_gamma x = - if x < 1. +. a then 1. -. p_gamma x + if x < 1. +. a then 1. -. p_gamma x else let rec qg_loop prev res la lb w k = if k > 1000. then failwith "q_gamma did not converge." @@ -195,26 +160,8 @@ let incomplete_gamma ~alpha x = qg_loop min_float (w /. lb) 1. lb w 2.0 in gf *. p_gamma x -(* Functions related to the Boys function:2 ends here *) - -(* The [[https://link.springer.com/article/10.1007/s10910-005-9023-3][Generalized Boys function]] is implemented, - * ~maxm~ is the maximum total angular momentum. - * - * \[ - * F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}} - * \] - * where $\gamma$ is the incomplete gamma function. - * - * - $F_0(0.) = 1$ - * - $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ - * - $F_m(0.) = \frac{1}{2m+1}$ - * - $F_m(t) = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}}$ - * - $F_m(t) = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ *) - - -(* [[file:~/QCaml/common/util.org::*Functions related to the Boys function][Functions related to the Boys function:4]] *) let boys_function ~maxm t = assert (t >= 0.); match maxm with @@ -231,8 +178,8 @@ let boys_function ~maxm t = Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) in let power_t_inv = (maxm+maxm+1) in - try - let fmax = + try + let fmax = let t_inv = sqrt (1. /. t) in let n = float_of_int maxm in let dm = 0.5 +. n in @@ -251,23 +198,15 @@ let boys_function ~maxm t = result with Exit -> result end -(* Functions related to the Boys function:4 ends here *) - -(* | ~list_some~ | Filters out all ~None~ elements of the list, and returns the elements without the ~Some~ | - * | ~list_range~ | Creates a list of consecutive integers | - * | ~list_pack~ | ~list_pack n l~ Creates a list of ~n~-elements lists | *) - - -(* [[file:~/QCaml/common/util.org::*List functions][List functions:2]] *) let list_some l = List.filter (function None -> false | _ -> true) l |> List.rev_map (function Some x -> x | _ -> assert false) |> List.rev -let list_range first last = +let list_range first last = if last < first then [] else let rec aux accu = function | 0 -> first :: accu @@ -281,7 +220,7 @@ let list_pack n l = let rec aux i accu1 accu2 = function | [] -> if accu1 = [] then List.rev accu2 - else + else List.rev ((List.rev accu1) :: accu2) | a :: rest -> match i with @@ -289,42 +228,27 @@ let list_pack n l = | _ -> (aux [@tailcall]) (i-1) (a::accu1) accu2 rest in aux (n-1) [] [] l -(* List functions:2 ends here *) -(* | ~array_range~ | Creates an array of consecutive integers | - * | ~array_sum~ | Returns the sum of all the elements of the array | - * | ~array_product~ | Returns the product of all the elements of the array | *) - - -(* [[file:~/QCaml/common/util.org::*Array functions][Array functions:2]] *) -let array_range first last = +let array_range first last = if last < first then [| |] else Array.init (last-first+1) (fun i -> i+first) -let array_sum a = +let array_sum a = Array.fold_left ( +. ) 0. a -let array_product a = +let array_product a = Array.fold_left ( *. ) 1. a -(* Array functions:2 ends here *) - -(* | ~seq_range~ | Creates a sequence returning consecutive integers | - * | ~seq_to_list~ | Read a sequence and put items in a list | - * | ~seq_fold~ | Apply a fold to the elements of the sequence | *) - - -(* [[file:~/QCaml/common/util.org::*Seq functions][Seq functions:2]] *) -let seq_range first last = +let seq_range first last = Seq.init (last-first) (fun i -> i+first) -let seq_to_list seq = - let rec aux accu xs = +let seq_to_list seq = + let rec aux accu xs = match Seq.uncons xs with | Some (x, xs) -> aux (x::accu) xs | None -> List.rev accu @@ -334,55 +258,24 @@ let seq_to_list seq = let seq_fold f init seq = Seq.fold_left f init seq -(* Seq functions:2 ends here *) - -(* | ~pp_float_array~ | Printer for float arrays | - * | ~pp_float_array_size~ | Printer for float arrays with size | - * | ~pp_float_2darray~ | Printer for matrices | - * | ~pp_float_2darray_size~ | Printer for matrices with size | - * | ~pp_bitstring~ | Printer for bit strings (used by ~Bitstring~ module) | - * - * Example: - * #+begin_example - * pp_float_array_size: - * [ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - * - * pp_float_array: - * [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - * - * pp_float_2darray_size - * [ - * 2:[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - * [ 4: 1.000000 2.000000 3.000000 4.000000 ] ] - * - * pp_float_2darray: - * [ [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - * [ 1.000000 2.000000 3.000000 4.000000 ] ] - * - * pp_bitstring 14: - * +++++------+-- - * #+end_example *) - - -(* [[file:~/QCaml/common/util.org::*Printers][Printers:2]] *) -let pp_float_array ppf a = +let pp_float_array ppf a = Format.fprintf ppf "@[<2>[@ "; Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; Format.fprintf ppf "]@]" -let pp_float_array_size ppf a = +let pp_float_array_size ppf a = Format.fprintf ppf "@[<2>@[ %d:@[<2>" (Array.length a); Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; Format.fprintf ppf "]@]@]" -let pp_float_2darray ppf a = +let pp_float_2darray ppf a = Format.fprintf ppf "@[<2>[@ "; Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array f) a; Format.fprintf ppf "]@]" -let pp_float_2darray_size ppf a = +let pp_float_2darray_size ppf a = Format.fprintf ppf "@[<2>@[ %d:@[" (Array.length a); Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array_size f) a; Format.fprintf ppf "]@]@]" @@ -390,4 +283,3 @@ let pp_float_2darray_size ppf a = let pp_bitstring n ppf bs = String.init n (fun i -> if (Z.testbit bs i) then '+' else '-') |> Format.fprintf ppf "@[%s@]" -(* Printers:2 ends here *) diff --git a/common/lib/util.mli b/common/lib/util.mli index 2474edb..c30f0d0 100644 --- a/common/lib/util.mli +++ b/common/lib/util.mli @@ -1,96 +1,157 @@ -(* [[file:~/QCaml/common/util.org::*Erf][Erf:2]] *) +(** Error function *) external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] -(* Erf:2 ends here *) -(* [[file:~/QCaml/common/util.org::*Erfc][Erfc:2]] *) external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] -(* Erfc:2 ends here *) -(* [[file:~/QCaml/common/util.org::*Gamma][Gamma:2]] *) external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] -(* Gamma:2 ends here *) -(* [[file:~/QCaml/common/util.org::*Popcnt][Popcnt:2]] *) + +(** Bit manipulation functions *) val popcnt : int64 -> int -(* Popcnt:2 ends here *) - -(* [[file:~/QCaml/common/util.org::*Trailz][Trailz:2]] *) val trailz : int64 -> int -(* Trailz:2 ends here *) - -(* [[file:~/QCaml/common/util.org::*Leadz][Leadz:2]] *) val leadz : int64 -> int -(* Leadz:2 ends here *) - -(* General functions *) -(* [[file:~/QCaml/common/util.org::*General functions][General functions:1]] *) +(** General functions *) + val fact : int -> float -(* @raise Invalid_argument for negative arguments or arguments >100. *) -val binom : int -> int -> int -val binom_float : int -> int -> float +(** Factorial function. + * @raise Invalid_argument for negative arguments or arguments >100. *) +val binom : int -> int -> int +(** Binomial coefficient. ~binom n k~ = $C_n^k$ *) + +val binom_float : int -> int -> float +(** Binomial coefficient (float). ~binom n k~ = $C_n^k$ *) + val chop : float -> (unit -> float) -> float +(** In ~chop a f~, evaluate ~f~ only if the absolute value of ~a~ is larger than ~Constants.epsilon~, and return ~a *. f ()~. *) + val pow : float -> int -> float +(** Fast implementation of the power function for small integer powers *) + val float_of_int_fast : int -> float +(** Faster implementation of float_of_int for small positive ints *) val of_some : 'a option -> 'a +(** Extracts the value of an option *) exception Not_implemented of string val not_implemented : string -> 'a -(* @raise Not_implemented. *) -(* General functions:1 ends here *) - -(* Functions related to the Boys function *) +(** Fails if some functionality is not implemented + * @raise Not_implemented. *) + + +(** Functions related to the Boys function *) -(* [[file:~/QCaml/common/util.org::*Functions related to the Boys function][Functions related to the Boys function:1]] *) val incomplete_gamma : alpha:float -> float -> float -(* @raise Failure when the calculation doesn't converge. *) -(* Functions related to the Boys function:1 ends here *) +(** The lower [[https://en.wikipedia.org/wiki/Incomplete_gamma_function][Incomplete Gamma function]] is implemented : + * \[ + * \gamma(\alpha,x) = \int_0^x e^{-t} t^{\alpha-1} dt + * \] + * + * p: $\frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ + * + * q: $\frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ + * + * Reference : Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten + * (New Algorithm handbook in C language) (Gijyutsu hyouron sha, + * Tokyo, 1991) p.227 [in Japanese]. + * + * @raise Failure when the calculation doesn't converge. *) -(* [[file:~/QCaml/common/util.org::*Functions related to the Boys function][Functions related to the Boys function:3]] *) val boys_function : maxm:int -> float -> float array -(* Functions related to the Boys function:3 ends here *) +(** Functions related to the Boys function:3 ends here -(* List functions *) + The [[https://link.springer.com/article/10.1007/s10910-005-9023-3][Generalized Boys function]] is implemented, + ~maxm~ is the maximum total angular momentum. + + \[ + F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}} + \] + where $\gamma$ is the incomplete gamma function. + + - $F_0(0.) = 1$ + - $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ + - $F_m(0.) = \frac{1}{2m+1}$ + - $F_m(t) = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}}$ + - $F_m(t) = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ +*) + + +(** List functions *) - -(* [[file:~/QCaml/common/util.org::*List functions][List functions:1]] *) val list_some : 'a option list -> 'a list +(** Filters out all ~None~ elements of the list, and returns the elements without the ~Some~ *) + val list_range : int -> int -> int list +(** Creates a list of consecutive integers *) + val list_pack : int -> 'a list -> 'a list list -(* List functions:1 ends here *) +(** ~list_pack n l~ Creates a list of ~n~-elements lists *) -(* Array functions *) +(** Array functions *) -(* [[file:~/QCaml/common/util.org::*Array functions][Array functions:1]] *) val array_range : int -> int -> int array -val array_sum : float array -> float -val array_product : float array -> float -(* Array functions:1 ends here *) +(** Creates an array of consecutive integers *) -(* Seq functions *) +val array_sum : float array -> float +(** Returns the sum of all the elements of the array *) + +val array_product : float array -> float +(** Returns the product of all the elements of the array *) -(* [[file:~/QCaml/common/util.org::*Seq functions][Seq functions:1]] *) +(** Seq functions *) + val seq_range : int -> int -> int Seq.t +(** Creates a sequence returning consecutive integers *) + val seq_to_list : 'a Seq.t -> 'a list +(** Read a sequence and put items in a list *) + val seq_fold : ('a -> 'b -> 'a) -> 'a -> 'b Seq.t -> 'a -(* Seq functions:1 ends here *) - -(* Printers *) +(** Apply a fold to the elements of the sequence *) -(* [[file:~/QCaml/common/util.org::*Printers][Printers:1]] *) +(** Printers *) + +(** + | ~pp_float_array~ | Printer for float arrays | + | ~pp_float_array_size~ | Printer for float arrays with size | + | ~pp_float_2darray~ | Printer for matrices | + | ~pp_float_2darray_size~ | Printer for matrices with size | + | ~pp_bitstring~ | Printer for bit strings (used by ~Bitstring~ module) | + + Example: + #+begin_example +pp_float_array_size: +[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + +pp_float_array: +[ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + +pp_float_2darray_size +[ + 2:[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + [ 4: 1.000000 2.000000 3.000000 4.000000 ] ] + +pp_float_2darray: +[ [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + [ 1.000000 2.000000 3.000000 4.000000 ] ] + +pp_bitstring 14: ++++++------+-- + #+end_example +*) + val pp_float_array_size : Format.formatter -> float array -> unit val pp_float_array : Format.formatter -> float array -> unit val pp_float_2darray_size : Format.formatter -> float array array -> unit val pp_float_2darray : Format.formatter -> float array array -> unit val pp_bitstring : int -> Format.formatter -> Z.t -> unit -(* Printers:1 ends here *) diff --git a/common/util.org b/common/util.org deleted file mode 100644 index 19b922a..0000000 --- a/common/util.org +++ /dev/null @@ -1,672 +0,0 @@ -#+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 c (concat lib name ".c")) -(setq test-ml (concat testdir name ".ml")) -(org-babel-tangle) -#+end_src - -* Util - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - Utility functions. - - -** Test header :noexport: - - #+begin_src ocaml :tangle (eval test-ml) :exports none -open Common.Util -open Alcotest - #+end_src - -** External C functions - - | ~erf_float~ | Error function ~erf~ from =libm= | - | ~erfc_float~ | Complementary error function ~erfc~ from =libm= | - | ~gamma_float~ | Gamma function ~gamma~ from =libm= | - | ~popcnt~ | ~popcnt~ instruction | - | ~trailz~ | ~ctz~ instruction | - | ~leadz~ | ~bsf~ instruction | - - #+begin_src c :tangle (eval c) :exports none -#include -#include -#include - #+end_src - -*** Erf - - #+begin_src c :tangle (eval c) :exports none -CAMLprim value erf_float_bytecode(value x) { - return caml_copy_double(erf(Double_val(x))); -} - -CAMLprim double erf_float(double x) { - return erf(x); -} - #+end_src - - #+begin_src ocaml :tangle (eval mli) -external erf_float : float -> float - = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -external erf_float : float -> float - = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] - #+end_src - -*** Erfc - - #+begin_src c :tangle (eval c) :exports none -CAMLprim value erfc_float_bytecode(value x) { - return caml_copy_double(erfc(Double_val(x))); -} - -CAMLprim double erfc_float(double x) { - return erfc(x); -} - #+end_src - - #+begin_src ocaml :tangle (eval mli) -external erfc_float : float -> float - = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] - #+end_src - -*** Gamma - - #+begin_src c :tangle (eval c) :exports none -CAMLprim value gamma_float_bytecode(value x) { - return caml_copy_double(tgamma(Double_val(x))); -} - - -CAMLprim double gamma_float(double x) { - return tgamma(x); -} - #+end_src - - #+begin_src ocaml :tangle (eval mli) -external gamma_float : float -> float - = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -external gamma_float : float -> float - = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] - #+end_src - -*** Popcnt - - #+begin_src c :tangle (eval c) :exports none -CAMLprim int32_t popcnt(int64_t i) { - return __builtin_popcountll (i); -} - - -CAMLprim value popcnt_bytecode(value i) { - return caml_copy_int32(__builtin_popcountll (Int64_val(i))); -} - #+end_src - - #+begin_src ocaml :tangle (eval mli) -val popcnt : int64 -> int - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -external popcnt : int64 -> int32 = "popcnt_bytecode" "popcnt" -[@@unboxed] [@@noalloc] - -let popcnt i = (popcnt [@inlined] ) i |> Int32.to_int - #+end_src - -*** Trailz - - #+begin_src c :tangle (eval c) :exports none -CAMLprim int32_t trailz(int64_t i) { - return i == 0L ? 64 : __builtin_ctzll (i); -} - - -CAMLprim value trailz_bytecode(value i) { - return caml_copy_int32(i == 0L ? 64 : __builtin_ctzll (Int64_val(i))); -} - #+end_src - - #+begin_src ocaml :tangle (eval mli) -val trailz : int64 -> int - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -external trailz : int64 -> int32 = "trailz_bytecode" "trailz" "int" -[@@unboxed] [@@noalloc] - -let trailz i = trailz i |> Int32.to_int - #+end_src - -*** Leadz - - #+begin_src c :tangle (eval c) :exports none -CAMLprim int32_t leadz(int64_t i) { - return i == 0L ? 64 : __builtin_clzll(i); -} - - -CAMLprim value leadz_bytecode(value i) { - return caml_copy_int32(i == 0L ? 64 : __builtin_clzll (Int64_val(i))); -} - #+end_src - - #+begin_src ocaml :tangle (eval mli) -val leadz : int64 -> int - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -external leadz : int64 -> int32 = "leadz_bytecode" "leadz" "int" -[@@unboxed] [@@noalloc] - -let leadz i = leadz i |> Int32.to_int - #+end_src - -*** Test - - #+begin_src ocaml :tangle (eval test-ml) :exports none -let test_external () = - check (float 1.e-15) "erf" 0.842700792949715 (erf_float 1.0); - check (float 1.e-15) "erf" 0.112462916018285 (erf_float 0.1); - check (float 1.e-15) "erf" (-0.112462916018285) (erf_float (-0.1)); - check (float 1.e-15) "erfc" 0.157299207050285 (erfc_float 1.0); - check (float 1.e-15) "erfc" 0.887537083981715 (erfc_float 0.1); - check (float 1.e-15) "erfc" (1.112462916018285) (erfc_float (-0.1)); - check (float 1.e-14) "gamma" (1.77245385090552) (gamma_float 0.5); - check (float 1.e-14) "gamma" (9.51350769866873) (gamma_float (0.1)); - check (float 1.e-14) "gamma" (-3.54490770181103) (gamma_float (-0.5)); - check int "popcnt" 6 (popcnt @@ Int64.of_int 63); - check int "popcnt" 8 (popcnt @@ Int64.of_int 299605); - check int "popcnt" 1 (popcnt @@ Int64.of_int 65536); - check int "popcnt" 0 (popcnt @@ Int64.of_int 0); - check int "trailz" 3 (trailz @@ Int64.of_int 8); - check int "trailz" 2 (trailz @@ Int64.of_int 12); - check int "trailz" 0 (trailz @@ Int64.of_int 1); - check int "trailz" 64 (trailz @@ Int64.of_int 0); - check int "leadz" 60 (leadz @@ Int64.of_int 8); - check int "leadz" 60 (leadz @@ Int64.of_int 12); - check int "leadz" 63 (leadz @@ Int64.of_int 1); - check int "leadz" 64 (leadz @@ Int64.of_int 0); - () - #+end_src - -** General functions - - #+begin_src ocaml :tangle (eval mli) -val fact : int -> float -(* @raise Invalid_argument for negative arguments or arguments >100. *) -val binom : int -> int -> int -val binom_float : int -> int -> float - -val chop : float -> (unit -> float) -> float -val pow : float -> int -> float -val float_of_int_fast : int -> float - -val of_some : 'a option -> 'a - -exception Not_implemented of string -val not_implemented : string -> 'a -(* @raise Not_implemented. *) - #+end_src - - | ~fact~ | Factorial function. | - | ~binom~ | Binomial coefficient. ~binom n k~ = $C_n^k$ | - | ~binom_float~ | float variant of ~binom~ | - | ~pow~ | Fast implementation of the power function for small integer powers | - | ~chop~ | In ~chop a f~, evaluate ~f~ only if the absolute value of ~a~ is larger than ~Constants.epsilon~, and return ~a *. f ()~. | - | ~float_of_int_fast~ | Faster implementation of float_of_int for small positive ints | - | ~not_implemented~ | Fails if some functionality is not implemented | - | ~of_some~ | Extracts the value of an option | - - #+begin_src ocaml :tangle (eval ml) :exports none -let memo_float_of_int = - Array.init 64 float_of_int - -let float_of_int_fast i = - if Int.logand i 63 = i then - memo_float_of_int.(i) - else - float_of_int i - - -let factmax = 150 - -let fact_memo = - let rec aux accu_l accu = function - | 0 -> (aux [@tailcall]) [1.] 1. 1 - | i when (i = factmax) -> - let x = (float_of_int factmax) *. accu in - List.rev (x::accu_l) - | i -> let x = (float_of_int i) *. accu in - (aux [@tailcall]) (x::accu_l) x (i+1) - in - aux [] 0. 0 - |> Array.of_list - -let fact = function - | i when (i < 0) -> - raise (Invalid_argument "Argument of factorial should be non-negative") - | i when (i > 150) -> - raise (Invalid_argument "Result of factorial is infinite") - | i -> fact_memo.(i) - - -let binom = - let memo = - let m = Array.make_matrix 64 64 0 in - for n=0 to Array.length m - 1 do - m.(n).(0) <- 1; - m.(n).(n) <- 1; - for k=1 to (n - 1) do - m.(n).(k) <- m.(n-1).(k-1) + m.(n-1).(k) - done - done; - m - in - let rec f n k = - assert (k >= 0); - assert (n >= k); - if k = 0 || k = n then - 1 - else if n < 64 then - memo.(n).(k) - else - f (n-1) (k-1) + f (n-1) k - in f - - -let binom_float n k = - binom n k - |> float_of_int_fast - - -let rec pow a = function - | 0 -> 1. - | 1 -> a - | 2 -> a *. a - | 3 -> a *. a *. a - | -1 -> 1. /. a - | n when n > 0 -> - let b = pow a (n / 2) in - b *. b *. (if n mod 2 = 0 then 1. else a) - | n when n < 0 -> (pow [@tailcall]) (1./.a) (-n) - | _ -> assert false - - -let chop f g = - if (abs_float f) < Constants.epsilon then 0. - else f *. (g ()) - - -exception Not_implemented of string -let not_implemented string = - raise (Not_implemented string) - - -let of_some = function - | Some a -> a - | None -> assert false - - #+end_src - - - #+begin_src ocaml :tangle (eval test-ml) :exports none -let test_general () = - check int "of_some_of_int_fast" 1 (of_some (Some 1)) ; - check int "binom" 35 (binom 7 4); - check (float 1.e-15) "fact" 5040. (fact 7); - check (float 1.e-15) "binom_float" 35.0 (binom_float 7 4); - check (float 1.e-15) "pow" 729.0 (pow 3.0 6); - check (float 1.e-15) "float_of_int_fast" 10.0 (float_of_int_fast 10); - () - #+end_src - -** Functions related to the Boys function - - #+begin_src ocaml :tangle (eval mli) -val incomplete_gamma : alpha:float -> float -> float -(* @raise Failure when the calculation doesn't converge. *) - #+end_src - - The lower [[https://en.wikipedia.org/wiki/Incomplete_gamma_function][Incomplete Gamma function]] is implemented : - \[ - \gamma(\alpha,x) = \int_0^x e^{-t} t^{\alpha-1} dt - \] - - p: $\frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ - - q: $\frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ - - Reference : Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten - (New Algorithm handbook in C language) (Gijyutsu hyouron sha, - Tokyo, 1991) p.227 [in Japanese] - - - #+begin_src ocaml :tangle (eval ml) :exports none -let incomplete_gamma ~alpha x = - assert (alpha >= 0.); - assert (x >= 0.); - let a = alpha in - let a_inv = 1./. a in - let gf = gamma_float alpha in - let loggamma_a = log gf in - let rec p_gamma x = - if x >= 1. +. a then 1. -. q_gamma x - else if x = 0. then 0. - else - let rec pg_loop prev res term k = - if k > 1000. then failwith "p_gamma did not converge." - else if prev = res then res - else - let term = term *. x /. (a +. k) in - (pg_loop [@tailcall]) res (res +. term) term (k +. 1.) - in - let r0 = exp (a *. log x -. x -. loggamma_a) *. a_inv in - pg_loop min_float r0 r0 1. - - and q_gamma x = - if x < 1. +. a then 1. -. p_gamma x - else - let rec qg_loop prev res la lb w k = - if k > 1000. then failwith "q_gamma did not converge." - else if prev = res then res - else - let k_inv = 1. /. k in - let kma = (k -. 1. -. a) *. k_inv in - let la, lb = - lb, kma *. (lb -. la) +. (k +. x) *. lb *. k_inv - in - let w = w *. kma in - let prev, res = res, res +. w /. (la *. lb) in - (qg_loop [@tailcall]) prev res la lb w (k +. 1.) - in - let w = exp (a *. log x -. x -. loggamma_a) in - let lb = (1. +. x -. a) in - qg_loop min_float (w /. lb) 1. lb w 2.0 - in - gf *. p_gamma x - #+end_src - - #+begin_src ocaml :tangle (eval mli) -val boys_function : maxm:int -> float -> float array - #+end_src - - The [[https://link.springer.com/article/10.1007/s10910-005-9023-3][Generalized Boys function]] is implemented, - ~maxm~ is the maximum total angular momentum. - - \[ - F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}} - \] - where $\gamma$ is the incomplete gamma function. - - - $F_0(0.) = 1$ - - $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ - - $F_m(0.) = \frac{1}{2m+1}$ - - $F_m(t) = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}}$ - - $F_m(t) = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ - - #+begin_src ocaml :tangle (eval ml) :exports none -let boys_function ~maxm t = - assert (t >= 0.); - match maxm with - | 0 -> - begin - if t = 0. then [| 1. |] else - let sq_t = sqrt t in - [| (Constants.sq_pi_over_two /. sq_t) *. erf_float sq_t |] - end - | _ -> - begin - assert (maxm > 0); - let result = - Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) - in - let power_t_inv = (maxm+maxm+1) in - try - let fmax = - let t_inv = sqrt (1. /. t) in - let n = float_of_int maxm in - let dm = 0.5 +. n in - let f = (pow t_inv power_t_inv ) in - match classify_float f with - | FP_normal -> (incomplete_gamma ~alpha:dm t) *. 0.5 *. f - | FP_zero - | FP_subnormal -> 0. - | _ -> raise Exit - in - let emt = exp (-. t) in - result.(maxm) <- fmax; - for n=maxm-1 downto 0 do - result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n) - done; - result - with Exit -> result - end - #+end_src - - #+begin_src ocaml :tangle (eval test-ml) :exports none -let test_boys () = - check (float 1.e-15) "incomplete_gamma" 0.0 (incomplete_gamma ~alpha:0.5 0.); - check (float 1.e-15) "incomplete_gamma" 1.114707979049507 (incomplete_gamma ~alpha:0.5 0.4); - check (float 1.e-15) "incomplete_gamma" 1.4936482656248544 (incomplete_gamma ~alpha:0.5 1.); - check (float 1.e-15) "incomplete_gamma" 1.7724401246392805 (incomplete_gamma ~alpha:0.5 10.); - check (float 1.e-15) "incomplete_gamma" 1.7724538509055159 (incomplete_gamma ~alpha:0.5 100.); - - check (float 1.e-15) "boys" 1.0 (boys_function ~maxm:0 0.).(0); - check (float 1.e-15) "boys" 0.2 (boys_function ~maxm:2 0.).(2); - check (float 1.e-15) "boys" (1./.3.) (boys_function ~maxm:2 0.).(1); - check (float 1.e-15) "boys" 0.8556243918921488 (boys_function ~maxm:0 0.5).(0); - check (float 1.e-15) "boys" 0.14075053682591263 (boys_function ~maxm:2 0.5).(2); - check (float 1.e-15) "boys" 0.00012711171070276764 (boys_function ~maxm:3 15.).(3); - () - #+end_src - -** List functions - - #+begin_src ocaml :tangle (eval mli) -val list_some : 'a option list -> 'a list -val list_range : int -> int -> int list -val list_pack : int -> 'a list -> 'a list list - #+end_src - - | ~list_some~ | Filters out all ~None~ elements of the list, and returns the elements without the ~Some~ | - | ~list_range~ | Creates a list of consecutive integers | - | ~list_pack~ | ~list_pack n l~ Creates a list of ~n~-elements lists | - - #+begin_src ocaml :tangle (eval ml) :exports none -let list_some l = - List.filter (function None -> false | _ -> true) l - |> List.rev_map (function Some x -> x | _ -> assert false) - |> List.rev - - -let list_range first last = - if last < first then [] else - let rec aux accu = function - | 0 -> first :: accu - | i -> (aux [@tailcall]) ( (first+i)::accu ) (i-1) - in - aux [] (last-first) - - -let list_pack n l = - assert (n>=0); - let rec aux i accu1 accu2 = function - | [] -> if accu1 = [] then - List.rev accu2 - else - List.rev ((List.rev accu1) :: accu2) - | a :: rest -> - match i with - | 0 -> (aux [@tailcall]) (n-1) [] ((List.rev (a::accu1)) :: accu2) rest - | _ -> (aux [@tailcall]) (i-1) (a::accu1) accu2 rest - in - aux (n-1) [] [] l - - #+end_src - - #+begin_src ocaml :tangle (eval test-ml) :exports none -let test_list () = - check bool "list_range" true ([ 2; 3; 4 ] = list_range 2 4); - check bool "list_some" true ([ 2; 3; 4 ] = - list_some ([ None ; Some 2 ; None ; Some 3 ; None ; None ; Some 4]) ); - check bool "list_pack" true (list_pack 3 (list_range 1 20) = - [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]; [10; 11; 12]; [13; 14; 15]; - [16; 17; 18]; [19; 20]]); - () - #+end_src - -** Array functions - - #+begin_src ocaml :tangle (eval mli) -val array_range : int -> int -> int array -val array_sum : float array -> float -val array_product : float array -> float - #+end_src - - | ~array_range~ | Creates an array of consecutive integers | - | ~array_sum~ | Returns the sum of all the elements of the array | - | ~array_product~ | Returns the product of all the elements of the array | - - #+begin_src ocaml :tangle (eval ml) :exports none -let array_range first last = - if last < first then [| |] else - Array.init (last-first+1) (fun i -> i+first) - - -let array_sum a = - Array.fold_left ( +. ) 0. a - - -let array_product a = - Array.fold_left ( *. ) 1. a - #+end_src - - #+begin_src ocaml :tangle (eval test-ml) :exports none -let test_array () = - check bool "array_range" true ([| 2; 3; 4 |] = array_range 2 4); - check (float 1.e-15) "array_sum" 9. (array_sum [| 2.; 3.; 4. |]); - check (float 1.e-15) "array_product" 24. (array_product [| 2.; 3.; 4. |]); - () - #+end_src - -** Seq functions - - #+begin_src ocaml :tangle (eval mli) -val seq_range : int -> int -> int Seq.t -val seq_to_list : 'a Seq.t -> 'a list -val seq_fold : ('a -> 'b -> 'a) -> 'a -> 'b Seq.t -> 'a - #+end_src - - | ~seq_range~ | Creates a sequence returning consecutive integers | - | ~seq_to_list~ | Read a sequence and put items in a list | - | ~seq_fold~ | Apply a fold to the elements of the sequence | - - #+begin_src ocaml :tangle (eval ml) :exports none -let seq_range first last = - Seq.init (last-first) (fun i -> i+first) - -let seq_to_list seq = - let rec aux accu xs = - match Seq.uncons xs with - | Some (x, xs) -> aux (x::accu) xs - | None -> List.rev accu - in - aux [] seq - - -let seq_fold f init seq = - Seq.fold_left f init seq - #+end_src - -** Printers - - #+begin_src ocaml :tangle (eval mli) -val pp_float_array_size : Format.formatter -> float array -> unit -val pp_float_array : Format.formatter -> float array -> unit -val pp_float_2darray_size : Format.formatter -> float array array -> unit -val pp_float_2darray : Format.formatter -> float array array -> unit -val pp_bitstring : int -> Format.formatter -> Z.t -> unit - #+end_src - - | ~pp_float_array~ | Printer for float arrays | - | ~pp_float_array_size~ | Printer for float arrays with size | - | ~pp_float_2darray~ | Printer for matrices | - | ~pp_float_2darray_size~ | Printer for matrices with size | - | ~pp_bitstring~ | Printer for bit strings (used by ~Bitstring~ module) | - - Example: - #+begin_example -pp_float_array_size: -[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - -pp_float_array: -[ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - -pp_float_2darray_size -[ - 2:[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - [ 4: 1.000000 2.000000 3.000000 4.000000 ] ] - -pp_float_2darray: -[ [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - [ 1.000000 2.000000 3.000000 4.000000 ] ] - -pp_bitstring 14: -+++++------+-- - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -let pp_float_array ppf a = - Format.fprintf ppf "@[<2>[@ "; - Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; - Format.fprintf ppf "]@]" - -let pp_float_array_size ppf a = - Format.fprintf ppf "@[<2>@[ %d:@[<2>" (Array.length a); - Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; - Format.fprintf ppf "]@]@]" - -let pp_float_2darray ppf a = - Format.fprintf ppf "@[<2>[@ "; - Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array f) a; - Format.fprintf ppf "]@]" - -let pp_float_2darray_size ppf a = - Format.fprintf ppf "@[<2>@[ %d:@[" (Array.length a); - Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array_size f) a; - Format.fprintf ppf "]@]@]" - -let pp_bitstring n ppf bs = - String.init n (fun i -> if (Z.testbit bs i) then '+' else '-') - |> Format.fprintf ppf "@[%s@]" - #+end_src - - -** Test footer :noexport: - - #+begin_src ocaml :tangle (eval test-ml) :exports none -let tests = [ - "External", `Quick, test_external; - "General" , `Quick, test_general; - "Boys" , `Quick, test_boys; - "List" , `Quick, test_list; - "Array" , `Quick, test_array; -] - #+end_src diff --git a/gaussian/lib/contracted_shell.ml b/gaussian/lib/contracted_shell.ml index b3846f1..9f29158 100644 --- a/gaussian/lib/contracted_shell.ml +++ b/gaussian/lib/contracted_shell.ml @@ -1,3 +1,5 @@ +[@@@landmark "auto-off"] + open Common type t = { @@ -16,7 +18,7 @@ module Co = Coordinate module Ps = Primitive_shell -let make ?(index=0) lc = +let[@landmark] make ?(index=0) lc = assert (Array.length lc > 0); let coef = Array.map fst lc @@ -30,7 +32,7 @@ let make ?(index=0) lc = 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 @@ -38,7 +40,7 @@ let make ?(index=0) lc = 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 = @@ -76,14 +78,14 @@ let primitives x = x.prim let zkey_array x = Ps.zkey_array x.prim.(0) -let values t point = +let[@landmark] 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 = + | i -> let new_accu = t.norm_coef.(i) *. t.coef.(i) *. exp(-. t.expo.(i) *. r2) +. accu in aux new_accu (i-1) in @@ -95,7 +97,7 @@ let values t point = 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 = + let fill arr v = arr.(0) <- 1.; for i=1 to n do arr.(i) <- arr.(i-1) *. v @@ -107,10 +109,10 @@ let values t point = 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 + t.norm_coef_scale.(i) *. x.(p.(0)) *. y.(p.(1)) *. z.(p.(2)) *. radial ) powers - + (** {2 Printers} *) diff --git a/gaussian/lib/contracted_shell_pair.ml b/gaussian/lib/contracted_shell_pair.ml index 438ec05..c4e171b 100644 --- a/gaussian/lib/contracted_shell_pair.ml +++ b/gaussian/lib/contracted_shell_pair.ml @@ -1,5 +1,6 @@ +[@@@landmark "auto-off"] open Common - + type t = { coefs_and_shell_pairs : (float * Primitive_shell_pair.t) list; @@ -18,11 +19,11 @@ module Psp = Primitive_shell_pair A contracted shell with N functions combined with a contracted shell with M functions generates a NxM array of shell pairs. *) -let make ?(cutoff=Constants.epsilon) s_a s_b = +let[@landmark] make ?(cutoff=Constants.epsilon) s_a s_b = let make = Psp.create_make_of (Cs.primitives s_a).(0) (Cs.primitives s_b).(0) in - let coefs_and_shell_pairs = + let coefs_and_shell_pairs = Array.mapi (fun i p_a -> let c_a = (Cs.coefficients s_a).(i) in let make = make p_a in @@ -49,15 +50,15 @@ let shell_a x = x.shell_a let shell_b x = x.shell_b let coefs_and_shell_pairs x = x.coefs_and_shell_pairs -let shell_pairs x = +let[@landmark] shell_pairs x = List.map snd x.coefs_and_shell_pairs |> Array.of_list -let coefficients x = +let[@landmark] coefficients x = List.map fst x.coefs_and_shell_pairs |> Array.of_list -let exponents_inv x = +let[@landmark] exponents_inv x = List.map (fun (_,sp) -> Psp.exponent_inv sp) x.coefs_and_shell_pairs |> Array.of_list @@ -66,17 +67,17 @@ let a_minus_b x = | [] -> assert false | (_,sp)::_ -> Psp.a_minus_b sp -let a_minus_b_sq x = +let a_minus_b_sq x = match x.coefs_and_shell_pairs with | [] -> assert false | (_,sp)::_ -> Psp.a_minus_b_sq sp -let ang_mom x = +let ang_mom x = match x.coefs_and_shell_pairs with | [] -> assert false | (_,sp)::_ -> Psp.ang_mom sp -let norm_scales x = +let norm_scales x = match x.coefs_and_shell_pairs with | [] -> assert false | (_,sp)::_ -> Psp.norm_scales sp @@ -89,18 +90,18 @@ let monocentric x = (** Returns an integer characteristic of a contracted shell pair *) -let hash a = +let[@landmark] hash a = List.rev_map Hashtbl.hash a |> Array.of_list (** Comparison function, used for sorting *) -let compare t t' = +let[@landmark] compare t t' = let a = hash t.coefs_and_shell_pairs in let b = hash t'.coefs_and_shell_pairs in if a = b then 0 else if (Array.length a < Array.length b) then -1 else if (Array.length a > Array.length b) then 1 - else + else let out = ref 0 in begin try @@ -118,11 +119,11 @@ let compare t t' = (** The array of all shell pairs with their correspondance in the list of contracted shells. *) -let of_contracted_shell_array ?(cutoff=Constants.epsilon) basis = +let[@landmark] of_contracted_shell_array ?(cutoff=Constants.epsilon) basis = let rec loop accu = function | [] -> accu | (s_a :: rest) as l -> - let new_accu = + let new_accu = (List.map (fun s_b -> make ~cutoff s_a s_b) l) :: accu in (loop [@tailcall]) new_accu rest in @@ -145,15 +146,15 @@ let equivalent x y = (** A list of unique shell pairs *) let unique sp = - let sp = + let sp = Array.to_list sp |> Array.concat |> Array.to_list in let rec aux accu = function | [] -> accu - | x::rest -> - let newaccu = + | x::rest -> + let newaccu = try ignore @@ List.find (fun y -> equivalent x y) accu; accu @@ -165,8 +166,8 @@ let unique sp = *) -let zkey_array x = - Am.zkey_array (Am.Doublet +let[@landmark] zkey_array x = + Am.zkey_array (Am.Doublet Cs.(ang_mom x.shell_a, ang_mom x.shell_b) ) diff --git a/gaussian/lib/contracted_shell_pair_couple.ml b/gaussian/lib/contracted_shell_pair_couple.ml index 7df7c6a..1a1ae54 100644 --- a/gaussian/lib/contracted_shell_pair_couple.ml +++ b/gaussian/lib/contracted_shell_pair_couple.ml @@ -1,6 +1,8 @@ +[@@@landmark "auto-off"] + open Common -type t = +type t = { coefs_and_shell_pair_couples : (float * Primitive_shell_pair_couple.t) list ; shell_pair_p: Contracted_shell_pair.t ; @@ -18,7 +20,7 @@ module Cs = Contracted_shell module Csp = Contracted_shell_pair module Pspc = Primitive_shell_pair_couple -let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q = +let[@landmark] make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q = let ang_mom = Am.(Csp.ang_mom shell_pair_p + Csp.ang_mom shell_pair_q) in @@ -29,8 +31,8 @@ let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q = in let cutoff = 1.e-3 *. cutoff in let coefs_and_shell_pair_couples = - List.concat_map (fun (c_ab, sp_ab) -> - List.map (fun (c_cd, sp_cd) -> + List.concat_map (fun (c_ab, sp_ab) -> + List.map (fun (c_cd, sp_cd) -> let coef_prod = c_ab *. c_cd in if abs_float coef_prod < cutoff then None else Some (coef_prod, Pspc.make sp_ab sp_cd) @@ -39,7 +41,7 @@ let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q = |> Util.list_some in match coefs_and_shell_pair_couples with - | [] -> None + | [] -> None | _ -> Some { shell_pair_p ; shell_pair_q ; ang_mom ; shell_a ; shell_b ; shell_c ; shell_d ; coefs_and_shell_pair_couples ; @@ -68,7 +70,7 @@ let zkey_array t = | (_,f)::_ -> Pspc.zkey_array f | _ -> invalid_arg "ContractedShellPairCouple.zkey_array" -let norm_scales t = +let norm_scales t = match t.coefs_and_shell_pair_couples with | (_,f)::_ -> Pspc.norm_scales f | _ -> invalid_arg "ContractedShellPairCouple.norm_scales" diff --git a/gaussian/lib/primitive_shell_pair.ml b/gaussian/lib/primitive_shell_pair.ml index 83c08c8..9f34029 100644 --- a/gaussian/lib/primitive_shell_pair.ml +++ b/gaussian/lib/primitive_shell_pair.ml @@ -1,3 +1,5 @@ +[@@@landmark "auto-off"] + open Common open Constants @@ -12,7 +14,7 @@ type t = { center : Coordinate.t; (* {% $P = (\alpha A + \beta B)/(\alpha+\beta)$ %} *) center_minus_a : Coordinate.t; (* {% $P - A$ %} *) a_minus_b : Coordinate.t; (* {% $A - B$ %} *) - ang_mom : Angular_momentum.t; + ang_mom : Angular_momentum.t; shell_a : Primitive_shell.t; shell_b : Primitive_shell.t; } @@ -21,124 +23,6 @@ module Am = Angular_momentum module Co = Coordinate module Ps = Primitive_shell - -let hash a = - Hashtbl.hash a - - -let equivalent a b = - a.exponent = b.exponent && - a.ang_mom = b.ang_mom && - a.normalization = b.normalization && - a.center = b.center && - a.center_minus_a = b.center_minus_a && - a.a_minus_b = b.a_minus_b - - -let cmp a b = - hash a - hash b - - -let compute_norm_scales p_a p_b = - Array.map (fun v1 -> - Array.map (fun v2 -> v1 *. v2) (Ps.norm_scales p_b) - ) (Ps.norm_scales p_a) - |> Array.to_list - |> Array.concat - -let create_make_of p_a p_b = - - let a_minus_b = - Co.( Ps.center p_a |- Ps.center p_b ) - in - - let a_minus_b_sq = - Co.dot a_minus_b a_minus_b - in - - let norm_scales = - lazy (compute_norm_scales p_a p_b) - in - - let ang_mom = - Am.( Ps.ang_mom p_a + Ps.ang_mom p_b ) - in - - function p_a -> - - let norm_coef_a = - Ps.normalization p_a - in - - let alfa_a = - Co.( Ps.exponent p_a |. Ps.center p_a ) - in - - function p_b -> - - let normalization = - norm_coef_a *. Ps.normalization p_b - in - - let exponent = - Ps.exponent p_a +. Ps.exponent p_b - in - - let exponent_inv = 1. /. exponent in - - let normalization = - let argexpo = - Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv - in - normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo) - in - - function cutoff -> - - if abs_float normalization > cutoff then ( - - let beta_b = - Co.( Ps.exponent p_b |. Ps.center p_b ) - in - - let center = - Co.(exponent_inv |. (alfa_a |+ beta_b)) - in - - let center_minus_a = - Co.(center |- Ps.center p_a) - in - - Some { - ang_mom ; - exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ; - a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a; - shell_b = p_b } - - ) - else None - - -let make p_a p_b = - let f = - create_make_of p_a p_b - in - match f p_a p_b 0. with - | Some result -> result - | None -> assert false - - -let norm_scales x = - try - Lazy.force x.norm_scales - with Lazy.Undefined -> compute_norm_scales x.shell_a x.shell_b - -let exponent_inv x = x.exponent_inv - -let monocentric x = - Ps.center x.shell_a = Ps.center x.shell_b - - let ang_mom x = x.ang_mom let a_minus_b x = x.a_minus_b @@ -158,8 +42,126 @@ let shell_a x = x.shell_a let shell_b x = x.shell_b -let zkey_array x = - Am.zkey_array (Am.Doublet +let hash a = + Hashtbl.hash a + + +let equivalent a b = + a.exponent = b.exponent && + a.ang_mom = b.ang_mom && + a.normalization = b.normalization && + a.center = b.center && + a.center_minus_a = b.center_minus_a && + a.a_minus_b = b.a_minus_b + + +let cmp a b = + hash a - hash b + + +let[@landmark] compute_norm_scales p_a p_b = + Array.map (fun v1 -> + Array.map (fun v2 -> v1 *. v2) (Ps.norm_scales p_b) + ) (Ps.norm_scales p_a) + |> Array.to_list + |> Array.concat + +let[@landmark] create_make_of p_a p_b = + + let a_minus_b = + Co.( Ps.center p_a |- Ps.center p_b ) + in + + let a_minus_b_sq = + Co.dot a_minus_b a_minus_b + in + + let norm_scales = + lazy (compute_norm_scales p_a p_b) + in + + let ang_mom = + Am.( Ps.ang_mom p_a + Ps.ang_mom p_b ) + in + + function p_a -> + + let norm_coef_a = + Ps.normalization p_a + in + + let alfa_a = + Co.( Ps.exponent p_a |. Ps.center p_a ) + in + + function p_b -> + + let normalization = + norm_coef_a *. Ps.normalization p_b + in + + let exponent = + Ps.exponent p_a +. Ps.exponent p_b + in + + let exponent_inv = 1. /. exponent in + + let normalization = + let argexpo = + Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv + in + normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo) + in + + function cutoff -> + + if abs_float normalization > cutoff then ( + + let beta_b = + Co.( Ps.exponent p_b |. Ps.center p_b ) + in + + let center = + Co.(exponent_inv |. (alfa_a |+ beta_b)) + in + + let center_minus_a = + Co.(center |- Ps.center p_a) + in + + Some { + ang_mom ; + exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ; + a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a; + shell_b = p_b } + + ) + else None + + +let[@landmark] make p_a p_b = + let f = + create_make_of p_a p_b + in + match f p_a p_b 0. with + | Some result -> result + | None -> assert false + + +let[@landmark] norm_scales x = + try + Lazy.force x.norm_scales + with Lazy.Undefined -> compute_norm_scales x.shell_a x.shell_b + +let exponent_inv x = x.exponent_inv + +let monocentric x = + Ps.center x.shell_a = Ps.center x.shell_b + + + +let[@landmark] zkey_array x = + Am.zkey_array (Am.Doublet Ps.(ang_mom x.shell_a, ang_mom x.shell_b) ) diff --git a/gaussian/lib/primitive_shell_pair_couple.ml b/gaussian/lib/primitive_shell_pair_couple.ml index 5508fb2..d605a9a 100644 --- a/gaussian/lib/primitive_shell_pair_couple.ml +++ b/gaussian/lib/primitive_shell_pair_couple.ml @@ -1,6 +1,8 @@ +[@@@landmark "auto-off"] + open Common - -type t = + +type t = { zkey_array : Zkey.t array lazy_t; shell_pair_p: Primitive_shell_pair.t ; @@ -17,7 +19,7 @@ module Co = Coordinate module Ps = Primitive_shell module Psp = Primitive_shell_pair -let make shell_pair_p shell_pair_q = +let[@landmark] make shell_pair_p shell_pair_q = let ang_mom = Am.(Psp.ang_mom shell_pair_p + Psp.ang_mom shell_pair_q) in @@ -27,9 +29,10 @@ let make shell_pair_p shell_pair_q = and shell_d = Psp.shell_b shell_pair_q in let zkey_array = lazy ( + let open Ps in Am.zkey_array (Am.Quartet - Ps.(ang_mom shell_a, ang_mom shell_b, - ang_mom shell_c, ang_mom shell_d) + (ang_mom shell_a, ang_mom shell_b, + ang_mom shell_c, ang_mom shell_d) ) ) in @@ -55,14 +58,14 @@ let shell_b t = t.shell_b let shell_c t = t.shell_c let shell_d t = t.shell_d -let p_minus_q t = +let p_minus_q t = let p = Psp.center t.shell_pair_p and q = Psp.center t.shell_pair_q in Co.(p |- q) let zkey_array t = Lazy.force t.zkey_array -let norm_scales t = +let[@landmark] norm_scales t = let norm_coef_scale_p_list = Array.to_list (Psp.norm_scales t.shell_pair_p) in let norm_coef_scale_q = Psp.norm_scales t.shell_pair_q in List.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) diff --git a/gaussian_integrals/lib/eri.ml b/gaussian_integrals/lib/eri.ml index a98e072..7637cbc 100644 --- a/gaussian_integrals/lib/eri.ml +++ b/gaussian_integrals/lib/eri.ml @@ -1,5 +1,7 @@ (** Electron-electron repulsion integrals *) +[@@@landmark "auto-off"] + open Common open Gaussian @@ -12,7 +14,16 @@ module T = struct open Zero_m_parameters - let zero_m z = + let rec aux_zero_m result expo_pq accu k = function + | 0 -> result.(k) <- result.(k) *. accu + | l -> + begin + result.(k) <- result.(k) *. accu; + let new_accu = -. accu *. expo_pq in + (aux_zero_m [@tailcall]) result expo_pq new_accu (k+1) (l-1) + end + + let[@landmark] zero_m z = let expo_pq_inv = z.expo_p_inv +. z.expo_q_inv in assert (expo_pq_inv <> 0.); let expo_pq = 1. /. expo_pq_inv in @@ -23,23 +34,14 @@ module T = struct in let maxm = z.maxm in let result = Util.boys_function ~maxm t in - let rec aux accu k = function - | 0 -> result.(k) <- result.(k) *. accu - | l -> - begin - result.(k) <- result.(k) *. accu; - let new_accu = -. accu *. expo_pq in - (aux [@tailcall]) new_accu (k+1) (l-1) - end - in let f = Constants.two_over_sq_pi *. (sqrt expo_pq) in - aux f 0 maxm; + aux_zero_m result expo_pq f 0 maxm; result - let class_of_contracted_shell_pair_couple ?operator shell_pair_couple = + let[@landmark] class_of_contracted_shell_pair_couple ?operator shell_pair_couple = assert (operator = None); let shell_p = Cspc.shell_pair_p shell_pair_couple - and shell_q = Cspc.shell_pair_q shell_pair_couple + and shell_q = Cspc.shell_pair_q shell_pair_couple in if Array.length (Csp.shell_pairs shell_p) + (Array.length (Csp.shell_pairs shell_q)) < 4 then @@ -51,6 +53,6 @@ module T = struct end -module M = Two_electron_integrals.Make(T) +module M = Two_electron_integrals.Make(T) include M diff --git a/gaussian_integrals/lib/two_electron_rr_vectorized.ml b/gaussian_integrals/lib/two_electron_rr_vectorized.ml index 216e817..c63e2d0 100644 --- a/gaussian_integrals/lib/two_electron_rr_vectorized.ml +++ b/gaussian_integrals/lib/two_electron_rr_vectorized.ml @@ -15,7 +15,7 @@ module Zp = Zero_m_parameters exception NullQuartet exception Found -let cutoff = Constants.integrals_cutoff +let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff let empty = Zmap.create 0 @@ -47,16 +47,16 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let expo_p_inv = abcd.expo_p_inv and expo_q_inv = abcd.expo_q_inv - and center_ab = abcd.center_ab - and center_cd = abcd.center_cd - and center_pq = abcd.center_pq + and center_ab = abcd.center_ab + and center_cd = abcd.center_cd + and center_pq = abcd.center_pq in let zero_m_array = abcd.zero_m_array in let maxm = Array.length zero_m_array - 1 in - let get_xyz angMom = + let get_xyz angMom = match angMom with | { Po.y=0 ; z=0 ; _ } -> Co.X | { z=0 ; _ } -> Co.Y @@ -64,7 +64,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in (* Vertical recurrence relations *) - let rec vrr0_v angMom_a = + let rec vrr0_v angMom_a = match angMom_a.Po.tot with | 0 -> zero_m_array | _ -> @@ -72,8 +72,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in try Zmap.find map_1d key with - | Not_found -> - let result = + | Not_found -> + let result = let xyz = get_xyz angMom_a in let am = Po.decr xyz angMom_a in let cab = Co.get xyz center_ab in @@ -83,11 +83,11 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) begin if abs_float cab >= cutoff then let expo_b = abcd.expo_b in - Array.iteri (fun m result_m -> + Array.iteri (fun m result_m -> let v0 = v_am.(m) in Array.iteri (fun l result_ml -> let f0 = -. expo_b.(l) *. expo_p_inv.(l) *. cab - and v0_l = v0.(l) + and v0_l = v0.(l) in Array.iteri (fun k v0_lk -> result_ml.(k) <- v0_lk *. f0) v0_l @@ -96,9 +96,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) end; let amxyz = Po.get xyz am in if amxyz < 1 then + let center_pq_xyz = center_pq xyz in Array.iteri (fun l expo_inv_p_l -> - let center_pq_xyz_l = (center_pq xyz).(l) in - Array.iteri (fun m result_m -> + let center_pq_xyz_l = center_pq_xyz.(l) in + Array.iteri (fun m result_m -> let result_ml = result_m.(l) in let p0 = v_am.(m+1) in let p0_l = p0.(l) @@ -114,10 +115,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let amm = Po.decr xyz am in let amxyz = Util.float_of_int_fast amxyz in let v_amm = vrr0_v amm in + let center_pq_xyz = center_pq xyz in Array.iteri (fun l expo_inv_p_l -> - let f = amxyz *. expo_p_inv.(l) *. 0.5 - and center_pq_xyz_l = (center_pq xyz).(l) - in + let f = amxyz *. expo_p_inv.(l) *. 0.5 + and center_pq_xyz_l = center_pq_xyz.(l) in Array.iteri (fun m result_m -> let v1 = v_amm.(m) in let v1_l = v1.(l) in @@ -129,12 +130,12 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) Array.iteri (fun k p0_lk -> result_ml.(k) <- result_ml.(k) +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +. - f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) + f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) ) p0.(l) ) result ) expo_p_inv end; - result + result in Zmap.add map_1d key result; result @@ -148,8 +149,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let key = Zkey.of_powers_six angMom_a angMom_c in try Zmap.find map_2d.(m) key with - | Not_found -> - let result = + | Not_found -> + let result = begin let xyz = get_xyz angMom_c in let cm = Po.decr xyz angMom_c in @@ -170,7 +171,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) if (!do_compute) then match vrr_v m angMom_a cm with | None -> None - | Some v1 -> + | Some v1 -> begin Some (Array.init np (fun l -> let v1_l = v1.(l) in @@ -182,8 +183,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let v2 = let f2 = + let center_pq_xyz = center_pq xyz in Array.init np (fun l -> - let cpq_l = (center_pq xyz).(l) in + let cpq_l = center_pq_xyz.(l) in Array.init nq (fun k -> let x = expo_q_inv.(k) *. cpq_l.(k) in if (!do_compute) then x @@ -193,12 +195,11 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) if (!do_compute) then match vrr_v (m+1) angMom_a cm with | None -> None - | Some v2 -> + | Some v2 -> begin for l=0 to np-1 do let f2_l = f2.(l) - and v2_l = v2.(l) - in + and v2_l = v2.(l) in for k=0 to nq-1 do f2_l.(k) <- -. v2_l.(k) *. f2_l.(k) done @@ -208,7 +209,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) else None in - let p1 = + let p1 = match v1, v2 with | None, None -> None | None, Some v2 -> Some v2 @@ -217,8 +218,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) begin for l=0 to np-1 do let v1_l = v1.(l) - and v2_l = v2.(l) - in + and v2_l = v2.(l) in for k=0 to nq-1 do v2_l.(k) <- v2_l.(k) +. v1_l.(k) done @@ -228,7 +228,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in let cxyz = Po.get xyz angMom_c in - let p2 = + let p2 = if cxyz < 2 then p1 else let cmm = Po.decr xyz cm in let fcm = (Util.float_of_int_fast (cxyz-1)) *. 0.5 in @@ -239,17 +239,16 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) else (if abs_float x > cutoff then do_compute := true ; x) ) in - let v1 = + let v1 = if (!do_compute) then match vrr_v m angMom_a cmm with | None -> None - | Some v1 -> + | Some v1 -> begin let result = Array.make_matrix np nq 0. in for l=0 to np-1 do let v1_l = v1.(l) - and result_l = result.(l) - in + and result_l = result.(l) in for k=0 to nq-1 do result_l.(k) <- v1_l.(k) *. f1.(k) done; @@ -259,9 +258,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) else None in - let v3 = - let f2 = - Array.init nq (fun k -> + let v3 = + let f2 = + Array.init nq (fun k -> let x = expo_q_inv.(k) *. f1.(k) in if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) @@ -270,7 +269,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) if (!do_compute) then match vrr_v (m+1) angMom_a cmm with | None -> None - | Some v3 -> + | Some v3 -> begin let result = Array.make_matrix np nq 0. in for l=0 to np-1 do @@ -295,8 +294,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) for l=0 to np-1 do let v3_l = v3.(l) and v1_l = v1.(l) - and p1_l = p1.(l) - in + and p1_l = p1.(l) in for k=0 to nq-1 do v3_l.(k) <- p1_l.(k) +. v1_l.(k) +. v3_l.(k) done @@ -307,8 +305,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) begin for l=0 to np-1 do let v1_l = v1.(l) - and p1_l = p1.(l) - in + and p1_l = p1.(l) in for k=0 to nq-1 do p1_l.(k) <- v1_l.(k) +. p1_l.(k) done @@ -319,8 +316,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) begin for l=0 to np-1 do let v3_l = v3.(l) - and p1_l = p1.(l) - in + and p1_l = p1.(l) in for k=0 to nq-1 do p1_l.(k) <- p1_l.(k) +. v3_l.(k) done @@ -331,8 +327,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) begin for l=0 to np-1 do let v3_l = v3.(l) - and v1_l = v1.(l) - in + and v1_l = v1.(l) in for k=0 to nq-1 do v3_l.(k) <- v1_l.(k) +. v3_l.(k) done @@ -342,7 +337,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in if (axyz < 1) || (cxyz < 1) then p2 else let am = Po.decr xyz angMom_a in - let v = + let v = vrr_v (m+1) am cm in match (p2, v) with @@ -380,7 +375,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let key = Zkey.of_powers_six angMom_a angMom_c in try Zmap.find map_2d.(0) key with - | Not_found -> + | Not_found -> let xyz = get_xyz angMom_c in let axyz = Po.get xyz angMom_a in let cm = Po.decr xyz angMom_c in @@ -426,22 +421,22 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) and result_l = result.(l) and expo_inv_q_over_p_l = expo_inv_q_over_p.(l) in - Array.iteri (fun k v1_lk -> + Array.iteri (fun k v1_lk -> let cqc = (center_qc xyz).(k) in result_l.(k) <- result_l.(k) +. - (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk + (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk ) v1_l ) v1 ; Some result) | None, None -> None | None, Some v1 -> Some (Array.init np (fun l -> - let v1_l = v1.(l) + let v1_l = v1.(l) and cpa = (center_pa xyz).(l) and expo_inv_q_over_p_l = expo_inv_q_over_p.(l) in Array.mapi (fun k v1_lk -> let cqc = (center_qc xyz).(k) in - (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk + (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk ) v1_l ) ) end @@ -457,7 +452,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result_l = result.(l) in Array.iteri (fun k v4_lk -> let expo_inv_q_over_p_l = expo_inv_q_over_p.(l) in - result_l.(k) <- result_l.(k) + result_l.(k) <- result_l.(k) -. expo_inv_q_over_p_l.(k) *. v4_lk) v4_l ) v4 ; Some result) | None, None -> None @@ -502,8 +497,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix in - let vrr_v a c = - let v = + let vrr_v a c = + let v = (* if c.Po.tot <> 0 then vrr_v 0 a c @@ -511,7 +506,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) *) vrr_v 0 a c in - match v with + match v with | Some matrix -> sum matrix | None -> 0. in @@ -553,20 +548,20 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | (_,0) -> if angMom_b.Po.tot = 0 then vrr_v angMom_a angMom_c else - hrr0_v angMom_a angMom_b angMom_c + hrr0_v angMom_a angMom_b angMom_c | (_,_) -> let xyz = get_xyz angMom_d in let cp = Po.incr xyz angMom_c in - let dm = Po.decr xyz angMom_d in - let h1 = - hrr_v angMom_a angMom_b cp dm + let dm = Po.decr xyz angMom_d in + let h1 = + hrr_v angMom_a angMom_b cp dm in let f = Co.get xyz center_cd in if abs_float f < cutoff then h1 else let h2 = - hrr_v angMom_a angMom_b angMom_c dm + hrr_v angMom_a angMom_b angMom_c dm in h1 +. f *. h2 in hrr_v angMom_a angMom_b angMom_c angMom_d @@ -585,7 +580,7 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell and cq = Csp.coefficients shell_q in - let np, nq = + let np, nq = Array.length sp, Array.length sq in @@ -593,7 +588,7 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell try match Cspc.make ~cutoff shell_p shell_q with | None -> raise NullQuartet - | Some shell_pair_couple -> + | Some shell_pair_couple -> let shell_a = Cspc.shell_a shell_pair_couple and shell_c = Cspc.shell_c shell_pair_couple @@ -605,7 +600,7 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell (* Pre-computation of integral class indices *) let class_indices = Cspc.zkey_array shell_pair_couple in - let contracted_class = + let contracted_class = Array.make (Array.length class_indices) 0.; in @@ -621,9 +616,9 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell contracted_class.(0) <- begin try - let expo_p_inv = - Vector.init np (fun ab -> Psp.exponent_inv sp.(ab-1)) - and expo_q_inv = + let expo_p_inv = + Vector.init np (fun ab -> Psp.exponent_inv sp.(ab-1)) + and expo_q_inv = Vector.init nq (fun cd -> Psp.exponent_inv sq.(cd-1)) in @@ -637,16 +632,14 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell if (abs_float coefx.{j,i} ) < 1.e-3*.cutoff then raise NullQuartet; - let expo_p_inv, expo_q_inv = - (expo_p_inv%.(i)), - (expo_q_inv%.(j)) - in + let expo_p_inv = expo_p_inv%.(i) in + let expo_q_inv = expo_q_inv%.(j) in let center_pq = Co.(Psp.center sp.(i-1) |- Psp.center sq.(j-1)) - and center_pa = + and center_pa = Co.(Psp.center sp.(i-1) |- Cs.center shell_a) - and center_qc = + and center_qc = Co.(Psp.center sq.(i-1) |- Cs.center shell_c) in let norm_pq_sq = @@ -654,32 +647,35 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell in let zero = Zp.zero ?operator zero_m in - let zero_m_array = + let zero_m_array = zero_m {zero with expo_p_inv ; expo_q_inv ; norm_pq_sq ; - center_pq ; center_pa ; center_qc ; + center_pq ; center_pa ; center_qc ; } in zero_m_array.(0) with NullQuartet -> 0. ) in - Matrix.gemm_trace zm_array coef + Matrix.gemm_trace zm_array coef with (Invalid_argument _) -> 0. end - | _ -> + | _ -> let coef = - Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) + Array.init np (fun l -> + let cpl = cp.(l) in + Array.map (fun cqk -> cqk *. cpl) cq + ) in let norm = Cspc.norm_scales shell_pair_couple in - let expo_p_inv = + let expo_p_inv = Array.map (fun shell_ab -> Psp.exponent_inv shell_ab) sp - and expo_q_inv = - Array.map (fun shell_cd -> Psp.exponent_inv shell_cd) sq + and expo_q_inv = + Array.map (fun shell_cd -> Psp.exponent_inv shell_cd) sq in let expo_b = @@ -688,99 +684,98 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell Array.map (fun shell_cd -> Ps.exponent (Psp.shell_b shell_cd) ) sq in - let center_pq = - let result = - Array.init 3 (fun xyz -> - Array.init np (fun ab -> - let shell_ab = sp.(ab) in - Array.init nq (fun cd -> - let shell_cd = sq.(cd) - in - let cpq = - Co.(Psp.center shell_ab |- Psp.center shell_cd) - in - match xyz with - | 0 -> Co.get X cpq; - | 1 -> Co.get Y cpq; - | _ -> Co.get Z cpq; - ) - ) - ) - in function - | Co.X -> result.(0) - | Co.Y -> result.(1) - | Co.Z -> result.(2) - in - let center_pa = - let result = - Array.init 3 (fun xyz -> - Array.init np (fun ab -> - let shell_ab = sp.(ab) in - let cpa = - Co.(Psp.center shell_ab |- Cs.center shell_a) - in - match xyz with - | 0 -> Co.(get X cpa); - | 1 -> Co.(get Y cpa); - | _ -> Co.(get Z cpa); - ) - ) - in function - | Co.X -> result.(0) - | Co.Y -> result.(1) - | Co.Z -> result.(2) - in - let center_qc = - let result = + let center_pq_x, center_pq_y, center_pq_z = + let result = Array.init 3 (fun xyz -> + let xyz = match xyz with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in + Array.init np (fun ab -> + let shell_ab = sp.(ab) in Array.init nq (fun cd -> - let shell_cd = sq.(cd) in - let cqc = - Co.(Psp.center shell_cd |- Cs.center shell_c) - in - match xyz with - | 0 -> Co.(get X cqc); - | 1 -> Co.(get Y cqc); - | _ -> Co.(get Z cqc); - ) + let shell_cd = sq.(cd) in + let cpq = + Co.(Psp.center shell_ab |- Psp.center shell_cd) + in + Co.get xyz cpq; + ) ) - in function - | Co.X -> result.(0) - | Co.Y -> result.(1) - | Co.Z -> result.(2) + ) + in + result.(0), result.(1), result.(2) + in + let center_pa_x, center_pa_y, center_pa_z = + let result = + Array.init 3 (fun xyz -> + let xyz = match xyz with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in + Array.init np (fun ab -> + let shell_ab = sp.(ab) in + let cpa = + Co.(Psp.center shell_ab |- Cs.center shell_a) + in + Co.get xyz cpa; + ) + ) + in result.(0), result.(1), result.(2) + in + let center_qc_x, center_qc_y, center_qc_z = + let result = + Array.init 3 (fun xyz -> + let xyz = match xyz with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in + Array.init nq (fun cd -> + let shell_cd = sq.(cd) in + let cqc = + Co.(Psp.center shell_cd |- Cs.center shell_c) + in + Co.get xyz cqc; + ) + ) + in result.(0), result.(1), result.(2) in let zero_m_array = - let result = + let result = Array.init (maxm+1) (fun _ -> - Array.init np (fun _ -> Array.make nq 0. ) ) + Array.make_matrix np nq 0. ) in let empty = Array.make (maxm+1) 0. in let center_qc_tmp = Array.init nq (fun cd -> Coordinate.make { Coordinate. - x = (center_qc Co.X).(cd) ; - y = (center_qc Co.Y).(cd) ; - z = (center_qc Co.Z).(cd) ; + x = center_qc_x.(cd) ; + y = center_qc_y.(cd) ; + z = center_qc_z.(cd) ; }) in Array.iteri (fun ab _shell_ab -> let center_pa = Coordinate.make { Coordinate. - x = (center_pa Co.X).(ab) ; - y = (center_pa Co.Y).(ab) ; - z = (center_pa Co.Z).(ab) ; + x = center_pa_x.(ab) ; + y = center_pa_y.(ab) ; + z = center_pa_z.(ab) ; } in + let coef_ab = coef.(ab) in + let expo_p_inv = expo_p_inv.(ab) in let zero_m_array_tmp = + let xab = center_pq_x.(ab) + and yab = center_pq_y.(ab) + and zab = center_pq_z.(ab) in Array.mapi (fun cd _shell_cd -> - if (abs_float coef.(ab).(cd) < cutoff) then + if (abs_float coef_ab.(cd) < cutoff) then empty else - let expo_p_inv, expo_q_inv = - expo_p_inv.(ab), expo_q_inv.(cd) - in - let x = (center_pq X).(ab).(cd) - and y = (center_pq Y).(ab).(cd) - and z = (center_pq Z).(ab).(cd) - in + let expo_q_inv = expo_q_inv.(cd) in + let x = xab.(cd) + and y = yab.(cd) + and z = zab.(cd) in let norm_pq_sq = x *. x +. y *. y +. z *. z in @@ -795,8 +790,7 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell (* Transpose result *) let coef_ab = coef.(ab) in for m=0 to maxm do - let result_m_ab = result.(m).(ab) - in + let result_m_ab = result.(m).(ab) in for cd=0 to nq-1 do result_m_ab.(cd) <- zero_m_array_tmp.(cd).(m) *. coef_ab.(cd) done @@ -828,18 +822,18 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell (* Schwartz screening *) if (np+nq> 24) then ( - let schwartz_p = + let schwartz_p = let key = Zkey.of_powers_twelve - angMom_a angMom_b angMom_a angMom_b + angMom_a angMom_b angMom_a angMom_b in match schwartz_p with | None -> 1. | Some schwartz_p -> Zmap.find schwartz_p key in if schwartz_p < cutoff then raise NullQuartet; - let schwartz_q = + let schwartz_q = let key = Zkey.of_powers_twelve - angMom_c angMom_d angMom_c angMom_d + angMom_c angMom_d angMom_c angMom_d in match schwartz_q with | None -> 1. @@ -848,7 +842,22 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; ); - let abcd = + let center_pq = function + | Co.X -> center_pq_x + | Co.Y -> center_pq_y + | Co.Z -> center_pq_z + in + let center_pa = function + | Co.X -> center_pa_x + | Co.Y -> center_pa_y + | Co.Z -> center_pa_z + in + let center_qc = function + | Co.X -> center_qc_x + | Co.Y -> center_qc_y + | Co.Z -> center_qc_z + in + let abcd = { expo_b ; expo_d ; expo_p_inv ; expo_q_inv ; center_ab = Csp.a_minus_b shell_p; center_cd = Csp.a_minus_b shell_q ; @@ -856,7 +865,7 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell center_qc ; zero_m_array } in - let integral = + let integral = hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) abcd map_1d map_2d np nq in @@ -866,7 +875,7 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell end; - let result = + let result = Zmap.create (Array.length contracted_class) in Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; diff --git a/gaussian_integrals/lib/zero_m_parameters.ml b/gaussian_integrals/lib/zero_m_parameters.ml index a94efaa..d02ea81 100644 --- a/gaussian_integrals/lib/zero_m_parameters.ml +++ b/gaussian_integrals/lib/zero_m_parameters.ml @@ -1,6 +1,6 @@ open Common open Operators - + type t = { expo_p_inv : float ; @@ -14,11 +14,11 @@ type t = operator : Operator.t option; } -let zero ?operator zero_m_func = +let zero ?operator zero_m_func = { zero_m_func ; operator ; - maxm=0 ; + maxm=0 ; expo_p_inv = 0.; expo_q_inv = 0.; norm_pq_sq = 0.;