From 8918f493029eee8572360b863fd983012b89cb57 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 24 Feb 2018 23:57:38 +0100 Subject: [PATCH] Documentation --- .ocamlinit | 28 ++++++ Basis/ERI.ml | 2 + Basis/KinInt.ml | 1 + Basis/NucInt.ml | 2 + Basis/OneElectronRR.ml | 9 +- Basis/Overlap.ml | 4 +- Basis/TwoElectronRR.ml | 4 +- Basis/TwoElectronRRVectorized.ml | 2 +- Makefile | 7 +- Nuclei/Element.ml | 4 +- Nuclei/Element.mli | 4 +- Nuclei/Mass.ml | 2 +- Utils/Angstrom.ml | 13 +-- Utils/Angstrom.mli | 3 +- Utils/AngularMomentum.mli | 87 +++++++++++++++++++ Utils/Bohr.ml | 3 +- Utils/Bohr.mli | 2 + Utils/Constants.ml | 5 +- Utils/Constants.mli | 20 ++++- Utils/Coordinate_type.ml | 12 ++- Utils/Coordinate_type.mli | 3 + Utils/Electrons.mli | 17 ++++ .../{PositiveFloat.ml => NonNegativeFloat.ml} | 1 + ...PositiveFloat.mli => NonNegativeFloat.mli} | 0 Utils/Radius.ml | 1 - Utils/Util.ml | 34 +++++--- Utils/Util.mli | 67 ++++++++++++++ 27 files changed, 276 insertions(+), 61 deletions(-) create mode 100644 .ocamlinit create mode 100644 Utils/AngularMomentum.mli create mode 100644 Utils/Electrons.mli rename Utils/{PositiveFloat.ml => NonNegativeFloat.ml} (81%) rename Utils/{PositiveFloat.mli => NonNegativeFloat.mli} (100%) delete mode 100644 Utils/Radius.ml create mode 100644 Utils/Util.mli diff --git a/.ocamlinit b/.ocamlinit new file mode 100644 index 0000000..55aa028 --- /dev/null +++ b/.ocamlinit @@ -0,0 +1,28 @@ +let () = + try Topdirs.dir_directory (Sys.getenv "OCAML_TOPLEVEL_PATH") + with Not_found -> () +;; + +#use "topfind";; +#require "lacaml";; +#directory "_build";; +#directory "_build/Utils";; +#directory "_build/Basis";; +#directory "_build/HartreeFock";; + +#load "Powers.cmo";; +#load "Zkey.cmo";; +#load "AngularMomentum.cmo";; +#load "Constants.cmo";; +#load "Util.cmo";; +#load "Coordinate.cmo";; +#load "Zmap.cmo";; +#load "Zkey.cmo";; +#load "ContractedShell.cmo";; +#load "ShellPair.cmo";; +#load "TwoElectronRR.cmo";; +#load "ERI.cmo";; + +#use "topfind";; + + diff --git a/Basis/ERI.ml b/Basis/ERI.ml index b707cfa..fa3385c 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -10,6 +10,8 @@ module Bs = Basis module Cs = ContractedShell module Csp = ContractedShellPair +let cutoff = integrals_cutoff + (** (00|00)^m : Fundamental electron repulsion integral $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $ diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index e82b3a0..3e6b02e 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -11,6 +11,7 @@ module Sp = ShellPair type t = Mat.t +let cutoff = integrals_cutoff (** Computes all the kinetic integrals of the contracted shell pair *) let contracted_class shell_a shell_b : float Zmap.t = diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index feacf29..8428032 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -33,6 +33,8 @@ let contracted_class_shell_pair shell_p geometry: float Zmap.t = OneElectronRR.contracted_class_shell_pair ~zero_m shell_p geometry +let cutoff = integrals_cutoff + let cutoff2 = cutoff *. cutoff exception NullIntegral diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 4b76d03..e6b1345 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -10,11 +10,6 @@ module Csp = ContractedShellPair module Po = Powers module Sp = ShellPair -(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *) -let chop f g = - if (abs_float f) < cutoff then 0. - else f *. (g ()) - (** Horizontal and Vertical Recurrence Relations (HVRR) *) @@ -84,7 +79,7 @@ let hvrr_one_e angMom_a angMom_b let bm = Po.decr xyz angMom_b in let h1 = hrr ap bm in let f2 = Co.get xyz center_ab in - if abs_float f2 < cutoff then h1 else + if abs_float f2 < integrals_cutoff then h1 else let h2 = hrr angMom_a bm in h1 +. f2 *. h2 @@ -130,7 +125,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let coef_prod = shell_p.Csp.coef.(ab) in (** Screening on the product of coefficients *) - if abs_float coef_prod < 1.e-3 *. cutoff then + if abs_float coef_prod < 1.e-3 *. integrals_cutoff then raise NullPair; diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index e2c6d71..c5985e2 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -11,6 +11,8 @@ module Cs = ContractedShell module Csp = ContractedShellPair module Sp = ShellPair +let cutoff = integrals_cutoff + (** Computes all the overlap integrals of the contracted shell pair *) let contracted_class shell_a shell_b : float Zmap.t = @@ -45,7 +47,7 @@ let contracted_class shell_a shell_b : float Zmap.t = shell_p.Csp.coef.(ab) in (** Screening on thr product of coefficients *) - if (abs_float coef_prod) > 1.e-4*.cutoff then + if (abs_float coef_prod) > 1.e-3*.cutoff then begin let expo_inv = shell_p.Csp.expo_inv.(ab) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 9684f04..99c48bd 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -1,5 +1,4 @@ open Util -open Constants module Am = AngularMomentum module Co = Coordinate @@ -8,8 +7,7 @@ module Csp = ContractedShellPair module Sp = ShellPair module Po = Powers -let debug=false - +let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff exception NullQuartet diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 903403e..3e40d37 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -11,7 +11,7 @@ module Po = Powers exception NullQuartet exception Found -let cutoff = Constants.cutoff +let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff let at_least_one_valid arr = diff --git a/Makefile b/Makefile index b7b083a..7449bba 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,9 @@ MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) Utils/math_functions.c MLYFILES=$(wildcard */*.mly) $(wildcard *.mly) MLFILES= $(wildcard */*.ml) $(wildcard *.ml) MLIFILES=$(wildcard */*.mli) $(wildcard *.mli) -ALL_EXE=$(patsubst %.ml,%.native,$(wildcard run_*.ml)) +ALL_NATIVE=$(patsubst %.ml,%.native,$(wildcard run_*.ml)) +ALL_BYTE=$(patsubst %.ml,%.byte,$(wildcard run_*.ml)) +ALL_EXE=$(ALL_BYTE) $(ALL_NATIVE) .PHONY: default @@ -32,7 +34,6 @@ doc: qpackage.odocl %.byte: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) rm -f -- $* $(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS) - ln -s $*.byte $* %.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) rm -f -- $* @@ -42,12 +43,10 @@ doc: qpackage.odocl %.p.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) rm -f -- $* $(OCAMLBUILD) $*.p.native -use-ocamlfind $(PKGS) - ln -s $*.p.native $* %.p.byte: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) rm -f -- $* $(OCAMLBUILD) -ocamlc ocamlcp $*.byte -use-ocamlfind $(PKGS) - ln -s $*.byte $* clean: rm -rf _build $(ALL_EXE) $(ALL_TESTS) *.native *.byte diff --git a/Nuclei/Element.ml b/Nuclei/Element.ml index 159671b..50275e1 100644 --- a/Nuclei/Element.ml +++ b/Nuclei/Element.ml @@ -141,7 +141,7 @@ let covalent_radius x = | Te -> 1.38 | I -> 1.39 | Xe -> 1.40 | Pt -> 1.30 in Constants.a0 *. (result x) - |> Radius.of_float + |> NonNegativeFloat.of_float let vdw_radius x = @@ -162,7 +162,7 @@ let vdw_radius x = | Te -> 2.06 | I -> 1.98 | Xe -> 2.16 | Pt -> 1.75 in Constants.a0 *. (result x) - |> Radius.of_float + |> NonNegativeFloat.of_float let mass c = diff --git a/Nuclei/Element.mli b/Nuclei/Element.mli index 3c78084..2f8cba7 100644 --- a/Nuclei/Element.mli +++ b/Nuclei/Element.mli @@ -20,6 +20,6 @@ val to_int : t -> int val of_int : int -> t val to_charge : t -> Charge.t val of_charge : Charge.t -> t -val covalent_radius : t -> Radius.t -val vdw_radius : t -> Radius.t +val covalent_radius : t -> NonNegativeFloat.t +val vdw_radius : t -> NonNegativeFloat.t val mass : t -> Mass.t diff --git a/Nuclei/Mass.ml b/Nuclei/Mass.ml index 035c570..98c4393 100644 --- a/Nuclei/Mass.ml +++ b/Nuclei/Mass.ml @@ -1 +1 @@ -include PositiveFloat +include NonNegativeFloat diff --git a/Utils/Angstrom.ml b/Utils/Angstrom.ml index b553427..0d3883b 100644 --- a/Utils/Angstrom.ml +++ b/Utils/Angstrom.ml @@ -1,12 +1 @@ -type t = { - x : float ; - y : float ; - z : float ; -} - - -let make { Point.x ; y ; z } = { x ; y ; z } -(* = %identity *) - - - +include Bohr diff --git a/Utils/Angstrom.mli b/Utils/Angstrom.mli index 77aa9e2..fa6fd98 100644 --- a/Utils/Angstrom.mli +++ b/Utils/Angstrom.mli @@ -1,3 +1,5 @@ +(** 3D point in cartesian coordinates, where units are Angstroms. *) + type t = private { x : float ; y : float ; @@ -7,4 +9,3 @@ type t = private { val make : Point.t -> t - diff --git a/Utils/AngularMomentum.mli b/Utils/AngularMomentum.mli new file mode 100644 index 0000000..dd0bea9 --- /dev/null +++ b/Utils/AngularMomentum.mli @@ -0,0 +1,87 @@ +(** Azimuthal quantum number, represented as s,p,d,... *) + +type t = S | P | D | F | G | H | I | J | K | L | M | N | O + +exception AngularMomentumError of string +(** Raised when the {!AngularMomentum.t} element can't be created. + *) + +val of_char : char -> t +(** Returns an {!AngularMomentum.t} when a shell is given as a character (case + insensitive). + + Example: + + [AngularMomentum.of_char 'p' -> AngularMomentum.P] + *) + +val to_string : t -> string +(** + [AngularMomentum.(to_string D) -> "D"] + *) + + +val to_char : t -> char +(** + [AngularMomentum.(to_char D) -> 'D'] + *) + + +val to_int : t -> int +(** + Returns the l{_max} value of the shell. + + Example: + + [AngularMomentum.of_int 3 -> AngularMomentum.F] + *) + +val of_int : int -> t +(** + Opposite of {!of_int}. + + Example: + + [AngularMomentum.of_int 3 -> AngularMomentum.F] + *) + +type kind = + Singlet of t + | Doublet of (t * t) + | Triplet of (t * t * t) + | Quartet of (t * t * t * t) + + +val n_functions : t -> int +(** Number of cartesian functions in shell. + + Example: + + [AngularMomentum.n_functions D -> 6] + *) + + +val zkey_array : kind -> Zkey.t array +(** Array of {!Zkey.t}, where each element is a a key associated with the + the powers of x,y,z. + + Example: + + {[ + AngularMomentum.( zkey_array Doublet (P,S) ) -> + [| {Zkey.left = 0; right = 1125899906842624} ; + {Zkey.left = 0; right = 1099511627776} ; + {Zkey.left = 0; right = 1073741824} |] + = + + let s,x,y,z = + Powers.( of_int_tuple (0,0,0), + of_int_tuple (1,0,0), + of_int_tuple (0,1,0), + of_int_tuple (0,0,1) ) + in + Array.map (fun (a,b) -> {!Zkey.of_powers_six} a b) + [| (x,s) ; (y,s) ; (z,s) |] + ]} + + *) diff --git a/Utils/Bohr.ml b/Utils/Bohr.ml index b553427..c51cda6 100644 --- a/Utils/Bohr.ml +++ b/Utils/Bohr.ml @@ -5,8 +5,7 @@ type t = { } -let make { Point.x ; y ; z } = { x ; y ; z } -(* = %identity *) +external make : Point.t -> t = "%identity" diff --git a/Utils/Bohr.mli b/Utils/Bohr.mli index 77aa9e2..dd46b40 100644 --- a/Utils/Bohr.mli +++ b/Utils/Bohr.mli @@ -1,3 +1,5 @@ +(** 3D point in cartesian coordinates, where units are atomic units (Bohr). *) + type t = private { x : float ; y : float ; diff --git a/Utils/Constants.ml b/Utils/Constants.ml index 9d8f9f9..d2b1b25 100644 --- a/Utils/Constants.ml +++ b/Utils/Constants.ml @@ -1,4 +1,5 @@ -let cutoff = 1.e-15 +let integrals_cutoff = 1.e-15 +let epsilon = 1.e-20 (** Constants *) let pi = acos (-1.) @@ -7,5 +8,5 @@ let sq_pi_over_two = sq_pi *. 0.5 let pi_inv = 1. /. pi let two_over_sq_pi = 2. /. (sqrt pi) -let a0 = 0.529_177_210_67 +let a0 = 0.529_177_210_671_2 let a0_inv = 1. /. a0 diff --git a/Utils/Constants.mli b/Utils/Constants.mli index 48541b6..73a1b08 100644 --- a/Utils/Constants.mli +++ b/Utils/Constants.mli @@ -1,8 +1,26 @@ -val cutoff : float +val epsilon : float +(** Value below which a float is considered null. Default is 10{^-20}. *) + +val integrals_cutoff : float +(** Cutoff value for integrals. Default is 10{^-15}. *) + val pi : float +(** pi = 3.141_592_653_589_793_12 *) + 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 +(** [ 2. /. (sqrt pi) ] *) + val a0 : float +(** Bohr radius : a{_0} = 0.529_177_210_671_2 Angstrom *) + val a0_inv : float +(** [ 1. /. a0 ] *) diff --git a/Utils/Coordinate_type.ml b/Utils/Coordinate_type.ml index b2af9a8..25f8504 100644 --- a/Utils/Coordinate_type.ml +++ b/Utils/Coordinate_type.ml @@ -1,12 +1,6 @@ -type point = Point.t - -type bohr = Bohr.t -type angstrom = Angstrom.t - -let a_to_b a = Constants.a0_inv *. a -let b_to_a b = Constants.a0 *. b let bohr_to_angstrom { Bohr.x ; y ; z } = + let b_to_a b = Constants.a0 *. b in Angstrom.make { Point. @@ -15,7 +9,10 @@ let bohr_to_angstrom { Bohr.x ; y ; z } = z = b_to_a z ; } + + let angstrom_to_bohr { Angstrom.x ; y ; z } = + let a_to_b a = Constants.a0_inv *. a in Bohr.make { Point. @@ -24,3 +21,4 @@ let angstrom_to_bohr { Angstrom.x ; y ; z } = z = a_to_b z ; } + diff --git a/Utils/Coordinate_type.mli b/Utils/Coordinate_type.mli index e3f1c00..5d2cbb2 100644 --- a/Utils/Coordinate_type.mli +++ b/Utils/Coordinate_type.mli @@ -1,3 +1,6 @@ +(** Coordinate_type can be Bohr or Angstrom. This module provides the conversion + functions from {!Bohr.t} to {!Angstrom.t} and from {!Angstrom.t} to {!Bohr.t} *) + type bohr = Bohr.t type angstrom = Angstrom.t diff --git a/Utils/Electrons.mli b/Utils/Electrons.mli new file mode 100644 index 0000000..c5f750c --- /dev/null +++ b/Utils/Electrons.mli @@ -0,0 +1,17 @@ +(** Information related to electrons. *) + +type t = { + n_alpha : int ; (** Number of alpha electrons *) + n_beta : int ; (** Number of beta electrons *) + multiplicity : int; (** Spin multiplicity: 2S+1 *) +} + + +val make : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t +(** Creates the data relative to electrons for a molecular system + described by {!Nuclei.t} for a given total charge and spin multiplicity. + @param multiplicity default is 1 + @param charge default is 0 + @raise Invalid_argument if the spin multiplicity is not compatible with + the molecule and the total charge. + *) diff --git a/Utils/PositiveFloat.ml b/Utils/NonNegativeFloat.ml similarity index 81% rename from Utils/PositiveFloat.ml rename to Utils/NonNegativeFloat.ml index fbf46c0..fe2bc8b 100644 --- a/Utils/PositiveFloat.ml +++ b/Utils/NonNegativeFloat.ml @@ -5,6 +5,7 @@ let of_float x = x external to_float : t -> float = "%identity" +external unsafe_of_float : float -> t = "%identity" let to_string x = let f = to_float x in string_of_float f diff --git a/Utils/PositiveFloat.mli b/Utils/NonNegativeFloat.mli similarity index 100% rename from Utils/PositiveFloat.mli rename to Utils/NonNegativeFloat.mli diff --git a/Utils/Radius.ml b/Utils/Radius.ml deleted file mode 100644 index 035c570..0000000 --- a/Utils/Radius.ml +++ /dev/null @@ -1 +0,0 @@ -include PositiveFloat diff --git a/Utils/Util.ml b/Utils/Util.ml index 5feb872..68819b2 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -1,13 +1,22 @@ -(** Functions from libm *) -external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] -external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] -external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] +(** All utilities which should be included in all source files are defined here *) + +(** {1 Functions from libm} *) open Constants open Lacaml.D let factmax = 150 +external erf_float : float -> float = "erf_float_bytecode" "erf_float" + [@@unboxed] [@@noalloc] + +external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" + [@@unboxed] [@@noalloc] + +external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" + [@@unboxed] [@@noalloc] + + @@ -79,8 +88,6 @@ let fact_memo = -(** Factorial function. - @raise Invalid_argument for negative or arguments >100. *) let fact = function | i when (i < 0) -> raise (Invalid_argument "Argument of factorial should be non-negative") @@ -89,24 +96,23 @@ let fact = function | i -> fact_memo.(i) - -(** Integer powers of floats *) let rec pow a = function | 0 -> 1. | 1 -> a | 2 -> a *. a | 3 -> a *. a *. a | -1 -> 1. /. a - | n when (n<0) -> pow (1./.a) (-n) - | n -> + | 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 (1./.a) (-n) + | _ -> assert false + + -(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *) let chop f g = - if (abs_float f) < cutoff then 0. + if (abs_float f) < Constants.epsilon then 0. else f *. (g ()) @@ -162,7 +168,7 @@ let diagonalize_symm h = in v, result -let xt_o_x o x = +let xt_o_x ~o ~x = gemm o x |> gemm ~transa:`T x diff --git a/Utils/Util.mli b/Utils/Util.mli new file mode 100644 index 0000000..4560f6c --- /dev/null +++ b/Utils/Util.mli @@ -0,0 +1,67 @@ +(** All utilities which should be included in all source files are defined here *) + +(** {1 Functions from libm} *) + +external erf_float : float -> float = "erf_float_bytecode" "erf_float" + [@@unboxed] [@@noalloc] + (** Error function [erf] from [libm] *) + +external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" + [@@unboxed] [@@noalloc] + (** Complementary error function [erfc] from [libm] *) + +external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" + [@@unboxed] [@@noalloc] + (** Gamma function [gamma] from [libm] *) + + + +(** {1 General functions} *) + +val fact : int -> float +(** Factorial function. + @raise Invalid_argument for negative arguments or arguments >100. + *) + +val pow : float -> int -> float +(** Fast implementation of the power function for small integer powers *) + +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 ()]. + *) + + + +(** {1 Functions related to the Boys function} *) + +val incomplete_gamma : alpha:float -> float -> float +(** {{:https://en.wikipedia.org/wiki/Incomplete_gamma_function} + Lower incomplete gamma function} + @raise Failure when the calculation doesn't converge. + *) + +val boys_function : maxm:int -> float -> float array +(** {{:https://link.springer.com/article/10.1007/s10910-005-9023-3} + Generalized Boys function}. + @param maxm Maximum total angular momentum. + *) + + +(** {1 Extension of the Array module} *) + +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 *) + + +(** {1 Linear algebra } *) + +val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec +(** Diagonalize a symmetric matrix. Returns the eigenvectors and the eigenvalues. *) + +val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat +(** Computes X{^T}.O.X *) +