LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

Remove org-mode

This commit is contained in:
Anthony Scemama 2024-01-17 10:30:24 +01:00
parent d53c0e6ea3
commit c4b022aeee
24 changed files with 710 additions and 2167 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ao/lib/dune
View File

@ -12,6 +12,7 @@
qcaml.gaussian_integrals
qcaml.operators
)
(instrumentation (backend landmarks))
(modules_without_implementation ao_dim)
)

96
common/charge.org
View File

@ -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

354
common/command_line.org
View File

@ -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 "<string>";
doc="Name of the file containing the basis set"; } ;
{ short='m' ; long="multiplicity" ; opt=Optional;
arg=With_arg "<int>";
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
- <<<Short option>>>: in the command line, a dash with a single character
(ex: =ls -l=)
- <<<Long option>>>: 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[=<n>]~)
#+begin_src ocaml :tangle (eval ml) :exports none
<<type>>
#+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 "@[<v 0>";
begin
match Str.split (Str.regexp "\n") t with
| x :: [] ->
Format.printf "@[<hov 0>";
Str.split (Str.regexp " ") x
|> List.iter (fun y -> Format.printf "@[%s@]@ " y) ;
Format.printf "@]"
| t ->
List.iter (fun x ->
Format.printf "@[<hov 0>";
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 <string> [-h] [-u <float>] -x <string> [--]
#+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=<string> 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 "@[<v 0>";
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 "@[<v>@[<v 2>Usage:@,@,@[<hov 4>@[%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 "@[<v 2>Arguments:@,";
Format.printf "@[<v 0>" ;
List.iter (fun x -> Format.printf "@ "; output_long max_width x) anon;
Format.printf "@]@,@]@,";
(* Print options and doc *)
Format.printf "@[<v 2>Options:@,";
Format.printf "@[<v 0>" ;
List.iter (fun x -> Format.printf "@ "; output_long max_width x) options;
Format.printf "@]@,@]@,";
(* Print footer *)
if !description_doc <> "" then
begin
Format.printf "@[<v 2>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

79
common/constants.org
View File

@ -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

218
common/coordinate.org
View File

@ -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

23
common/lib/charge.ml
View File

@ -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 *)

27
common/lib/charge.mli
View File

@ -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 *)

71
common/lib/command_line.ml
View File

@ -1,17 +1,4 @@
(* - <<<Short option>>>: in the command line, a dash with a single character
* (ex: =ls -l=)
* - <<<Long option>>>: 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[=<n>]~) *)
(* [[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 "@[<v 0>";
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 <string> [-h] [-u <float>] -x <string> [--]
* #+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=<string> 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 *)

98
common/lib/command_line.mli
View File

@ -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 "<string>";
* doc="Name of the file containing the basis set"; } ;
*
* { short='m' ; long="multiplicity" ; opt=Optional;
* arg=With_arg "<int>";
* 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
* ...
*)
(* - <<<Short option>>>: in the command line, a dash with a single character
* (ex: =ls -l=)
* - <<<Long option>>>: 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[=<n>]`) *)
(** 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. *)

29
common/lib/constants.ml
View File

@ -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 *)

36
common/lib/constants.mli
View File

@ -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~ *)

80
common/lib/coordinate.ml
View File

@ -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 *)

91
common/lib/coordinate.mli
View File

@ -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 *)

152
common/lib/util.ml
View File

@ -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 "@[<h>%s@]"
(* Printers:2 ends here *)

151
common/lib/util.mli
View File

@ -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 *)

672
common/util.org
View File

@ -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 <math.h>
#include <caml/mlvalues.h>
#include <caml/alloc.h>
#+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 "@[<h>%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

18
gaussian/lib/contracted_shell.ml
View File

@ -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} *)

39
gaussian/lib/contracted_shell_pair.ml
View File

@ -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)
)

14
gaussian/lib/contracted_shell_pair_couple.ml
View File

@ -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"

244
gaussian/lib/primitive_shell_pair.ml
View File

@ -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)
)

17
gaussian/lib/primitive_shell_pair_couple.ml
View File

@ -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)

30
gaussian_integrals/lib/eri.ml
View File

@ -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

331
gaussian_integrals/lib/two_electron_rr_vectorized.ml
View File

@ -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;

6
gaussian_integrals/lib/zero_m_parameters.ml
View File

@ -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.;