From e1da54cd67c1e352d4c7bd4f9b10c300d9066f99 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Feb 2018 18:41:30 +0100 Subject: [PATCH] Working on documenation --- Basis/Basis.ml | 11 +- Basis/Basis.mli | 33 +++-- Basis/Basis_lexer.mll | 33 ----- Basis/ContractedShellPair.ml | 8 +- Basis/Contracted_shell.ml | 129 ------------------- Basis/Contracted_shell.mli | 27 ---- Basis/ERI.ml | 2 +- Basis/Gamess_parser.mly | 49 -------- Basis/Gamess_reader.ml | 18 --- Basis/General_basis.ml | 40 ------ Basis/KinInt.ml | 4 +- Basis/OneElectronRR.ml | 4 +- Basis/Overlap.ml | 4 +- Basis/ShellPair.ml | 6 +- Basis/TwoElectronRR.ml | 4 +- Basis/TwoElectronRRVectorized.ml | 206 ++++++++++++++++--------------- Nuclei/Element.ml | 4 +- Nuclei/Mass.ml | 2 +- Utils/Angular_momentum.ml | 134 -------------------- Utils/Positive_float.ml | 14 --- Utils/Positive_float.mli | 5 - Utils/Radius.ml | 2 +- Utils/Units.ml | 25 ---- Utils/Units.mli | 8 -- 24 files changed, 151 insertions(+), 621 deletions(-) delete mode 100644 Basis/Basis_lexer.mll delete mode 100644 Basis/Contracted_shell.ml delete mode 100644 Basis/Contracted_shell.mli delete mode 100644 Basis/Gamess_parser.mly delete mode 100644 Basis/Gamess_reader.ml delete mode 100644 Basis/General_basis.ml delete mode 100644 Utils/Angular_momentum.ml delete mode 100644 Utils/Positive_float.ml delete mode 100644 Utils/Positive_float.mli delete mode 100644 Utils/Units.ml delete mode 100644 Utils/Units.mli diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 5259dc7..3795ba1 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -1,11 +1,12 @@ type t = { size : int; - contracted_shells : Contracted_shell.t array; + contracted_shells : ContractedShell.t array; } -module Cs = Contracted_shell -module Gb = General_basis +module Cs = ContractedShell +module Gb = GeneralBasis + (** Returns an array of the basis set per atom *) let of_nuclei_and_general_basis n b = @@ -62,9 +63,11 @@ let to_string b = ) ^ line + + let of_nuclei_and_basis_filename ~nuclei ~filename = let general_basis = - Gamess_reader.read ~filename + GamessReader.read ~filename in of_nuclei_and_general_basis nuclei general_basis diff --git a/Basis/Basis.mli b/Basis/Basis.mli index 5aa65cb..3b9dc2c 100644 --- a/Basis/Basis.mli +++ b/Basis/Basis.mli @@ -1,22 +1,29 @@ +(** The atomic basis set is represented as an array of {!ContractedShell.t}. *) + type t = private { - (** Number of contracted Gaussians *) - size : int; - - (** Array of contracted shells *) - contracted_shells : Contracted_shell.t array; + size : int ; (** Number of contracted Gaussians *) + contracted_shells : + ContractedShell.t array ; (** Contracted shells *) } - -(** Returns an array of the basis set per atom *) -val of_nuclei_and_general_basis : Nuclei.t -> General_basis.t list -> t - - -(** Pretty prints the basis set in a string *) val to_string : t -> string +(** Pretty prints the basis set in a string. *) + + +val of_nuclei_and_general_basis : Nuclei.t -> GeneralBasis.t list -> t +(** Takes an array of {!Nuclei.t}, and a {!GeneralBasis.t} (such as cc-pVDZ + for instance) and creates the corresponding atomic basis set. + All the {!Element.t}s of the array of {!Nuclei.t} are searched in + the {!GeneralBasis.t}, and the basis is built by creating + {!ContractedShell.t}s centered on the nuclei with the exponents + and contraction coefficients given by the {!GeneralBasis.t}. +*) + -(** Create a basis using the coordinates of Nuclei and a the filename of - the general basis set *) val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> filename:string -> t +(** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t} + from a file. + *) diff --git a/Basis/Basis_lexer.mll b/Basis/Basis_lexer.mll deleted file mode 100644 index 2540494..0000000 --- a/Basis/Basis_lexer.mll +++ /dev/null @@ -1,33 +0,0 @@ -{ -exception SyntaxError of string - -open Gamess_parser - -} - -let eol = ['\n'] -let white = [' ' '\t']+ -let element = ['A'-'Z' 'a'-'z']+ white? eol -let ang_mom = ['S' 'P' 'D' 'F' 'G' 'H' 'I' 'J' 'K' 'L' 'M' 'N' 'O' - 's' 'p' 'd' 'f' 'g' 'h' 'i' 'j' 'k' 'l' 'm' 'n' 'o' ] - white -let integer = ['0'-'9']+ -let real = '-'? integer '.' integer (['e' 'E'] ('+'|'-')? integer)? - - -rule read_all_rule = parse - | eol { EOL } - | white { SPACE } - | ang_mom as a { ANG_MOM (a.[0]) } - | element as e { ELEMENT (String.trim e) } - | integer as i { INTEGER (int_of_string i) } - | real as f { FLOAT (float_of_string f) } - | eof { EOF } - - -{ - let rec read_all lexbuf = - match read_all_rule lexbuf with - | SPACE -> read_all_rule lexbuf - | x -> x -} diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 9fb51f7..f7f43a8 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -5,8 +5,8 @@ exception Null_contribution type t = { - shell_a : Contracted_shell.t; - shell_b : Contracted_shell.t; + shell_a : ContractedShell.t; + shell_b : ContractedShell.t; shell_pairs : ShellPair.t array; coef : float array; expo_inv : float array; @@ -18,9 +18,9 @@ type t = } -module Am = Angular_momentum +module Am = AngularMomentum module Co = Coordinate -module Cs = Contracted_shell +module Cs = ContractedShell module Sp = ShellPair (** Creates an contracted shell pair : an array of pairs of primitive shells. diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml deleted file mode 100644 index d47af6a..0000000 --- a/Basis/Contracted_shell.ml +++ /dev/null @@ -1,129 +0,0 @@ -open Util -open Constants -open Coordinate - -type shell_contracted = { - expo : float array; (* Gaussian exponents *) - coef : float array; (* Contraction coefficients *) - center : Coordinate.t; (* Center of all the Gaussians *) - totAngMom : Angular_momentum.t; (* Total angular momentum *) - size : int; (* Number of contracted Gaussians *) - norm_coef : float array; (* Normalization coefficient of the class - corresponding to the i-th contraction *) - norm_coef_scale : float array; (* Inside a class, the norm is the norm - of the function with (totAngMom,0,0) *. - this scaling factor *) - index : int; (* Index in the array of contracted shells *) - powers : Zkey.t array; (* Array of Zkeys corresponding to the - powers of (x,y,z) in the class *) -} - -type t = shell_contracted - -module Am = Angular_momentum - -(** Normalization coefficient of contracted function i, which depends on the - exponent and the angular momentum. Two conventions can be chosen : a single - normalisation factor for all functions of the class, or a coefficient which - depends on the powers of x,y and z. - Returns, for each contracted function, an array of functions taking as - argument the [|x;y;z|] powers. -*) -let compute_norm_coef expo totAngMom = - let atot = - Am.to_int totAngMom - in - let factor int_array = - let dfa = Array.map (fun j -> - ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) - ) int_array - in - sqrt (dfa.(0) *.dfa.(1) *. dfa.(2)) - in - let expo = - if atot mod 2 = 0 then - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2)) - ) expo - else - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot) - ) expo - in - Array.map (fun x -> let f a = x *. (factor a) in f) expo - - -let make ~index ~expo ~coef ~center ~totAngMom = - assert (Array.length expo = Array.length coef); - assert (Array.length expo > 0); - let norm_coef_func = - compute_norm_coef expo totAngMom - in - let powers = - Am.zkey_array (Am.Singlet totAngMom) - in - let norm_coef = - Array.map (fun f -> f [| Am.to_int totAngMom ; 0 ; 0 |]) norm_coef_func - in - let norm_coef_scale = - Array.map (fun a -> - (norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0) - ) powers - in - { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ; - norm_coef_scale ; powers } - - -let with_index a i = - { a with index = i } - - -let to_string s = - let coord = s.center in - let open Printf in - (match s.totAngMom with - | Am.S -> sprintf "%3d " (s.index+1) - | _ -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.powers)) - ) ^ - ( sprintf "%1s %8.3f %8.3f %8.3f " (Am.to_string s.totAngMom) - (get X coord) (get Y coord) (get Z coord) ) ^ - (Array.map2 (fun e c -> sprintf "%16.8e %16.8e" e c) s.expo s.coef - |> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) - - -(** Normalization coefficient of contracted function i, which depends on the - exponent and the angular momentum. Two conventions can be chosen : a single - normalisation factor for all functions of the class, or a coefficient which - depends on the powers of x,y and z. - Returns, for each contracted function, an array of functions taking as - argument the [|x;y;z|] powers. -*) -let compute_norm_coef expo totAngMom = - let atot = - Am.to_int totAngMom - in - let factor int_array = - let dfa = Array.map (fun j -> - (float_of_int (1 lsl j) *. fact j) /. fact (j+j) - ) int_array - in - sqrt (dfa.(0) *.dfa.(1) *. dfa.(2)) - in - let expo = - if atot mod 2 = 0 then - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2)) - ) expo - else - Array.map (fun alpha -> - let alpha_2 = alpha +. alpha in - (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot) - ) expo - in - Array.map (fun x -> let f a = x *. factor a in f) expo - - - diff --git a/Basis/Contracted_shell.mli b/Basis/Contracted_shell.mli deleted file mode 100644 index 6f1a73b..0000000 --- a/Basis/Contracted_shell.mli +++ /dev/null @@ -1,27 +0,0 @@ -type shell_contracted = private { - expo : float array; - coef : float array; - center : Coordinate.t; - totAngMom : Angular_momentum.t; - size : int; - norm_coef : float array; - norm_coef_scale : float array; - index : int; - powers : Zkey.t array; -} - -type t = shell_contracted - - -(** Pretty-printing of the contracted shell in a string *) -val to_string : t -> string - -(** Creates a contracted shell *) -val make : - index:int -> - expo:float array -> - coef:float array -> - center:Coordinate.t -> totAngMom:Angular_momentum.t -> t - -(** Returns a copy of the contracted shell with a modified index *) -val with_index : t -> int -> t diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 044aca5..b707cfa 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -7,7 +7,7 @@ open Bigarray type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t module Bs = Basis -module Cs = Contracted_shell +module Cs = ContractedShell module Csp = ContractedShellPair (** (00|00)^m : Fundamental electron repulsion integral diff --git a/Basis/Gamess_parser.mly b/Basis/Gamess_parser.mly deleted file mode 100644 index 082e4f8..0000000 --- a/Basis/Gamess_parser.mly +++ /dev/null @@ -1,49 +0,0 @@ -/* Parses basis sets GAMESS format */ - -%{ - -%} - -%token ELEMENT -%token ANG_MOM -%token INTEGER -%token FLOAT -%token SPACE -%token EOL -%token EOF - -%start input -%type input - -%% /* Grammar rules and actions follow */ - -input: - | basis { $1 } - | EOL input { $2 } - -basis: - | element shell_array EOL { ($1, $2) } - | element shell_array EOF { ($1, $2) } - -element: - | ELEMENT { Element.of_string $1 } - -shell_array: - | shell_list { Array.of_list @@ List.rev $1 } - -shell_list: - | { [] } - | shell_list shell { $2 :: $1 } - -shell: - | ANG_MOM INTEGER EOL primitive_list { (Angular_momentum.of_char $1, Array.of_list @@ List.rev $4 ) } - -primitive_list: - | { [] } - | primitive_list primitive { $2 :: $1 } - -primitive: - | INTEGER FLOAT FLOAT EOL { General_basis.{exponent=$2 ; coefficient=$3 } } - - - diff --git a/Basis/Gamess_reader.ml b/Basis/Gamess_reader.ml deleted file mode 100644 index e2a2b7e..0000000 --- a/Basis/Gamess_reader.ml +++ /dev/null @@ -1,18 +0,0 @@ -(** Read a basis set file in GAMESS format and return an association list where the key is an - Element.t and the value is the parsed basis set. *) -let read ~filename = - let lexbuf = - let ic = open_in filename in - Lexing.from_channel ic - in - let rec aux accu = - try - let key, basis = - Gamess_parser.input Basis_lexer.read_all lexbuf - in - aux ((key, basis)::accu) - with - | Parsing.Parse_error -> List.rev accu - in - aux [] - diff --git a/Basis/General_basis.ml b/Basis/General_basis.ml deleted file mode 100644 index 752d9b1..0000000 --- a/Basis/General_basis.ml +++ /dev/null @@ -1,40 +0,0 @@ -(** General basis set read from a file *) -type primitive = { - exponent: float ; - coefficient: float - } - -type general_contracted_shell = Angular_momentum.t * (primitive array) - -type t = Element.t * (general_contracted_shell array) - - -module Am = Angular_momentum - -let string_of_primitive ?id prim = - match id with - | None -> (string_of_float prim.exponent)^" "^(string_of_float prim.coefficient) - | Some i -> (string_of_int i)^" "^(string_of_float prim.exponent)^" "^(string_of_float prim.coefficient) - - -let string_of_contracted_shell (angular_momentum, prim_array) = - let n = - Array.length prim_array - in - Printf.sprintf "%s %d\n%s" - (Am.to_string angular_momentum) n - (Array.init n (fun i -> string_of_primitive ~id:(i+1) prim_array.(i)) - |> Array.to_list - |> String.concat "\n") - - -let string_of_contracted_shell_array a = - Array.map string_of_contracted_shell a - |> Array.to_list - |> String.concat "\n" - - -let to_string (name, contracted_shell_array) = - Printf.sprintf "%s\n%s" name (string_of_contracted_shell_array contracted_shell_array) - - diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 78c9325..e82b3a0 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -2,10 +2,10 @@ open Util open Constants open Lacaml.D -module Am = Angular_momentum +module Am = AngularMomentum module Bs = Basis module Co = Coordinate -module Cs = Contracted_shell +module Cs = ContractedShell module Csp = ContractedShellPair module Sp = ShellPair diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 065240b..4b76d03 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -3,9 +3,9 @@ open Constants exception NullPair -module Am = Angular_momentum +module Am = AngularMomentum module Co = Coordinate -module Cs = Contracted_shell +module Cs = ContractedShell module Csp = ContractedShellPair module Po = Powers module Sp = ShellPair diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index b312db3..e2c6d71 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -4,10 +4,10 @@ open Lacaml.D type t = Mat.t -module Am = Angular_momentum +module Am = AngularMomentum module Bs = Basis module Co = Coordinate -module Cs = Contracted_shell +module Cs = ContractedShell module Csp = ContractedShellPair module Sp = ShellPair diff --git a/Basis/ShellPair.ml b/Basis/ShellPair.ml index 98d5a99..9660f0f 100644 --- a/Basis/ShellPair.ml +++ b/Basis/ShellPair.ml @@ -14,8 +14,8 @@ type t = { totAngMomInt : int ; i : int; j : int; - shell_a : Contracted_shell.t; - shell_b : Contracted_shell.t; + shell_a : ContractedShell.t; + shell_b : ContractedShell.t; monocentric : bool } @@ -27,7 +27,7 @@ let hash a = let equivalent a b = a = b (* - Hashtbl.hash (a.expo, a.center_a, a.center_ab, a.coef, Contracted_shell.totAngMom a.shell_a, Contracted_shell.totAngMom a.shell_b) + Hashtbl.hash (a.expo, a.center_a, a.center_ab, a.coef, ContractedShell.totAngMom a.shell_a, ContractedShell.totAngMom a.shell_b) *) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 3fdd502..9684f04 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -1,9 +1,9 @@ open Util open Constants -module Am = Angular_momentum +module Am = AngularMomentum module Co = Coordinate -module Cs = Contracted_shell +module Cs = ContractedShell module Csp = ContractedShellPair module Sp = ShellPair module Po = Powers diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index f3225d0..903403e 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -1,16 +1,19 @@ open Util open Lacaml.D open Bigarray -open Powers -open Coordinate -open Contracted_shell_type -let cutoff = Constants.cutoff -let cutoff2 = cutoff *. cutoff +module Am = AngularMomentum +module Co = Coordinate +module Csp = ContractedShellPair +module Sp = ShellPair +module Po = Powers exception NullQuartet exception Found +let cutoff = Constants.cutoff +let cutoff2 = cutoff *. cutoff + let at_least_one_valid arr = try Array.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false @@ -31,14 +34,14 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let get_xyz angMom = match angMom with - | { y=0 ; z=0 ; _ } -> X - | { z=0 ; _ } -> Y - | _ -> Z + | { Po.y=0 ; z=0 ; _ } -> Co.X + | { z=0 ; _ } -> Co.Y + | _ -> Co.Z in (** Vertical recurrence relations *) let rec vrr0_v angMom_a = - match angMom_a.tot with + match angMom_a.Po.tot with | 0 -> zero_m_array | _ -> let key = Zkey.of_powers_three angMom_a @@ -48,9 +51,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | Not_found -> let result = let xyz = get_xyz angMom_a in - let am = Powers.decr xyz angMom_a in - let cab = Coordinate.get xyz center_ab in - let result = Array.init (maxm+1-angMom_a.tot) (fun _ -> Array.make_matrix np nq 0.) in + let am = Po.decr xyz angMom_a in + let cab = Co.get xyz center_ab in + let result = Array.init (maxm+1-angMom_a.Po.tot) (fun _ -> Array.make_matrix np nq 0.) in let v_am= vrr0_v am in begin @@ -66,7 +69,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) ) result_m ) result end; - let amxyz = Powers.get xyz am in + let amxyz = Po.get xyz am in if amxyz < 1 then Array.iteri (fun l expo_inv_p_l -> let center_pq_xyz_l = (center_pq xyz).(l) in @@ -83,7 +86,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) ) expo_inv_p else begin - let amm = Powers.decr xyz am in + let amm = Po.decr xyz am in let amxyz = float_of_int amxyz in let v_amm = vrr0_v amm in Array.iteri (fun l expo_inv_p_l -> @@ -113,7 +116,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) and vrr_v m angMom_a angMom_c = - match (angMom_a.tot, angMom_c.tot) with + match (angMom_a.Po.tot, angMom_c.Po.tot) with | (i,0) -> Some (vrr0_v angMom_a).(m) | (_,_) -> @@ -124,12 +127,12 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = begin let xyz = get_xyz angMom_c in - let cm = Powers.decr xyz angMom_c in - let axyz = Powers.get xyz angMom_a in + let cm = Po.decr xyz angMom_c in + let axyz = Po.get xyz angMom_a in let do_compute = ref false in let v1 = - let f = -. (Coordinate.get xyz center_cd) in + let f = -. (Co.get xyz center_cd) in let f1 = Array.init nq (fun k -> @@ -198,10 +201,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) end in - let cxyz = Powers.get xyz angMom_c in + let cxyz = Po.get xyz angMom_c in let p2 = if cxyz < 2 then p1 else - let cmm = Powers.decr xyz cm in + let cmm = Po.decr xyz cm in let fcm = (float_of_int (cxyz-1)) *. 0.5 in let f1 = Array.init nq (fun k -> @@ -312,7 +315,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) end in if (axyz < 1) || (cxyz < 1) then p2 else - let am = Powers.decr xyz angMom_a in + let am = Po.decr xyz angMom_a in let v = vrr_v (m+1) am cm in @@ -344,7 +347,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (* and trr_v angMom_a angMom_c = - match (angMom_a.tot, angMom_c.tot) with + match (angMom_a.Po.tot, angMom_c.Po.tot) with | (i,0) -> Some (vrr0_v angMom_a).(0) | (_,_) -> @@ -353,9 +356,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) try Zmap.find map_2d.(0) key with | Not_found -> let xyz = get_xyz angMom_c in - let axyz = Powers.get xyz angMom_a in - let cm = Powers.decr xyz angMom_c in - let cmxyz = Powers.get xyz cm in + let axyz = Po.get xyz angMom_a in + let cm = Po.decr xyz angMom_c in + let cmxyz = Po.get xyz cm in let expo_inv_q_over_p = Array.mapi (fun l expo_inv_p_l -> let expo_p_l = 1./.expo_inv_p_l in @@ -368,7 +371,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) if cmxyz < 1 then result else begin let f = 0.5 *. (float_of_int cmxyz) in - let cmm = Powers.decr xyz cm in + let cmm = Po.decr xyz cm in match result, trr_v angMom_a cmm with | None, None -> None | None, Some v3 -> @@ -420,7 +423,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = if cmxyz < 0 then result else begin - let ap = Powers.incr xyz angMom_a in + let ap = Po.incr xyz angMom_a in match result, trr_v ap cm with | Some result, None -> Some result | Some result, Some v4 -> @@ -445,7 +448,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) if axyz < 1 then result else begin let f = 0.5 *. (float_of_int axyz) in - let am = Powers.decr xyz angMom_a in + let am = Po.decr xyz angMom_a in match result, trr_v am cm with | Some result, None -> Some result | Some result, Some v2 -> @@ -476,7 +479,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let vrr_v a c = let v = (* - if c.tot <> 0 then + if c.Po.tot <> 0 then vrr_v 0 a c else trr_v a c *) @@ -491,48 +494,48 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (** Horizontal recurrence relations *) let rec hrr0_v angMom_a angMom_b angMom_c = - match angMom_b.tot with + match angMom_b.Po.tot with | 0 -> begin - match (angMom_a.tot, angMom_c.tot) with + match (angMom_a.Po.tot, angMom_c.Po.tot) with | (0,0) -> sum zero_m_array.(0) | (_,_) -> vrr_v angMom_a angMom_c end | 1 -> let xyz = get_xyz angMom_b in - let ap = Powers.incr xyz angMom_a in - let f = Coordinate.get xyz center_ab in + let ap = Po.incr xyz angMom_a in + let f = Co.get xyz center_ab in let v1 = vrr_v ap angMom_c in if (abs_float f < cutoff) then v1 else let v2 = vrr_v angMom_a angMom_c in v1 +. v2 *. f | _ -> let xyz = get_xyz angMom_b in - let bxyz = Powers.get xyz angMom_b in + let bxyz = Po.get xyz angMom_b in if (bxyz < 0) then 0. else - let ap = Powers.incr xyz angMom_a in - let bm = Powers.decr xyz angMom_b in + let ap = Po.incr xyz angMom_a in + let bm = Po.decr xyz angMom_b in let h1 = hrr0_v ap bm angMom_c in - let f = Coordinate.get xyz center_ab in + let f = Co.get xyz center_ab in if abs_float f < cutoff then h1 else let h2 = hrr0_v angMom_a bm angMom_c in h1 +. h2 *. f and hrr_v angMom_a angMom_b angMom_c angMom_d = - match (angMom_b.tot, angMom_d.tot) with - | (_,0) -> if angMom_b.tot = 0 then + match (angMom_b.Po.tot, angMom_d.Po.tot) with + | (_,0) -> if angMom_b.Po.tot = 0 then vrr_v angMom_a angMom_c else hrr0_v angMom_a angMom_b angMom_c | (_,_) -> let xyz = get_xyz angMom_d in - let cp = Powers.incr xyz angMom_c in - let dm = Powers.decr 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 in - let f = Coordinate.get xyz center_cd in + let f = Co.get xyz center_cd in if abs_float f < cutoff then h1 else @@ -550,24 +553,23 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - let shell_a = shell_p.ContractedShellPair.shell_a - and shell_b = shell_p.ContractedShellPair.shell_b - and shell_c = shell_q.ContractedShellPair.shell_a - and shell_d = shell_q.ContractedShellPair.shell_b - and sp = shell_p.ContractedShellPair.shell_pairs - and sq = shell_q.ContractedShellPair.shell_pairs + let shell_a = shell_p.Csp.shell_a + and shell_b = shell_p.Csp.shell_b + and shell_c = shell_q.Csp.shell_a + and shell_d = shell_q.Csp.shell_b + and sp = shell_p.Csp.shell_pairs + and sq = shell_q.Csp.shell_pairs in let maxm = - shell_p.ContractedShellPair.totAngMomInt + - shell_q.ContractedShellPair.totAngMomInt + shell_p.Csp.totAngMomInt + + shell_q.Csp.totAngMomInt in (* Pre-computation of integral class indices *) let class_indices = - Angular_momentum.zkey_array - (Angular_momentum.Quartet - (shell_a.totAngMom, shell_b.totAngMom, - shell_c.totAngMom, shell_d.totAngMom)) + Am.zkey_array (Am.Quartet + (shell_a.totAngMom, shell_b.totAngMom, + shell_c.totAngMom, shell_d.totAngMom)) in let contracted_class = @@ -575,21 +577,21 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let monocentric = - shell_p.ContractedShellPair.monocentric && - shell_q.ContractedShellPair.monocentric && - shell_p.ContractedShellPair.shell_a.center = - shell_q.ContractedShellPair.shell_a.center + shell_p.Csp.monocentric && + shell_q.Csp.monocentric && + shell_p.Csp.shell_a.center = + shell_q.Csp.shell_a.center in (** Screening on the product of coefficients *) let coef_max_p = Array.fold_left (fun accu x -> if (abs_float x) > accu then (abs_float x) else accu) - 0. shell_p.ContractedShellPair.coef + 0. shell_p.Csp.coef and coef_max_q = Array.fold_left (fun accu x -> if (abs_float x) > accu then (abs_float x) else accu) - 0. shell_q.ContractedShellPair.coef + 0. shell_q.Csp.coef in let rec build_list cutoff vec accu = function @@ -599,10 +601,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q else accu ) (k-1) in let p_list = - let vec = shell_p.ContractedShellPair.coef in + let vec = shell_p.Csp.coef in build_list (cutoff /. coef_max_q) vec [] (Array.length vec - 1) and q_list = - let vec = shell_q.ContractedShellPair.coef in + let vec = shell_q.Csp.coef in build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1) in @@ -623,21 +625,21 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q begin match (shell_a.totAngMom, shell_b.totAngMom, shell_c.totAngMom, shell_d.totAngMom) with - | Angular_momentum.(S,S,S,S) -> + | Am.(S,S,S,S) -> contracted_class.(0) <- begin try let expo_inv_p = - Vec.init np (fun ab -> sp.(ab-1).ShellPair.expo_inv) + Vec.init np (fun ab -> sp.(ab-1).Sp.expo_inv) and expo_inv_q = - Vec.init nq (fun cd -> sq.(cd-1).ShellPair.expo_inv) + Vec.init nq (fun cd -> sq.(cd-1).Sp.expo_inv) in let coef = let result = Mat.make0 nq np in Lacaml.D.ger - (Vec.of_array @@ filter_q shell_q.ContractedShellPair.coef) - (Vec.of_array @@ filter_p shell_p.ContractedShellPair.coef) + (Vec.of_array @@ filter_q shell_q.Csp.coef) + (Vec.of_array @@ filter_p shell_p.Csp.coef) result; result in @@ -651,10 +653,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let center_pq = - sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center + Co.(sp.(i-1).Sp.center |- sq.(j-1).Sp.center) in let norm_pq_sq = - Coordinate.dot center_pq center_pq + Co.dot center_pq center_pq in let zero_m_array = @@ -669,24 +671,24 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q | _ -> let coef = - let cp = filter_p shell_p.ContractedShellPair.coef - and cq = filter_q shell_q.ContractedShellPair.coef + let cp = filter_p shell_p.Csp.coef + and cq = filter_q shell_q.Csp.coef in Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) in let expo_inv_p = - Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp + Array.map (fun shell_ab -> shell_ab.Sp.expo_inv) sp and expo_inv_q = - Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq + Array.map (fun shell_cd -> shell_cd.Sp.expo_inv) sq in let expo_b = - Array.map (fun shell_ab -> shell_b.expo.(shell_ab.ShellPair.j)) sp + Array.map (fun shell_ab -> shell_b.expo.(shell_ab.Sp.j)) sp and expo_d = - Array.map (fun shell_cd -> shell_d.expo.(shell_cd.ShellPair.j)) sq + Array.map (fun shell_cd -> shell_d.expo.(shell_cd.Sp.j)) sq in - let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in + let norm_coef_scale_p = shell_p.Csp.norm_coef_scale in let center_pq = let result = @@ -697,19 +699,19 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let shell_cd = sq.(cd) in let cpq = - shell_ab.ShellPair.center |- shell_cd.ShellPair.center + Co.(shell_ab.Sp.center |- shell_cd.Sp.center) in match xyz with - | 0 -> Coordinate.get X cpq; - | 1 -> Coordinate.get Y cpq; - | _ -> Coordinate.get Z cpq; + | 0 -> Co.get X cpq; + | 1 -> Co.get Y cpq; + | _ -> Co.get Z cpq; ) ) ) in function - | X -> result.(0) - | Y -> result.(1) - | Z -> result.(2) + | Co.X -> result.(0) + | Co.Y -> result.(1) + | Co.Z -> result.(2) in let center_pa = let result = @@ -717,18 +719,18 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.init np (fun ab -> let shell_ab = sp.(ab) in let cpa = - shell_ab.ShellPair.center |- shell_a.center + Co.(shell_ab.Sp.center |- shell_a.center) in match xyz with - | 0 -> Coordinate.get X cpa; - | 1 -> Coordinate.get Y cpa; - | _ -> Coordinate.get Z cpa; + | 0 -> Co.(get X cpa); + | 1 -> Co.(get Y cpa); + | _ -> Co.(get Z cpa); ) ) in function - | X -> result.(0) - | Y -> result.(1) - | Z -> result.(2) + | Co.X -> result.(0) + | Co.Y -> result.(1) + | Co.Z -> result.(2) in let center_qc = let result = @@ -736,18 +738,18 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.init nq (fun cd -> let shell_cd = sq.(cd) in let cqc = - shell_cd.ShellPair.center |- shell_c.center + Co.(shell_cd.Sp.center |- shell_c.center) in match xyz with - | 0 -> Coordinate.get X cqc; - | 1 -> Coordinate.get Y cqc; - | _ -> Coordinate.get Z cqc; + | 0 -> Co.(get X cqc); + | 1 -> Co.(get Y cqc); + | _ -> Co.(get Z cqc); ) ) in function - | X -> result.(0) - | Y -> result.(1) - | Z -> result.(2) + | Co.X -> result.(0) + | Co.Y -> result.(1) + | Co.Z -> result.(2) in let zero_m_array = let result = @@ -787,7 +789,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm = let norm_coef_scale_q = - shell_q.ContractedShellPair.norm_coef_scale + shell_q.Csp.norm_coef_scale in Array.to_list norm_coef_scale_p |> List.map (fun v1 -> @@ -843,8 +845,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q zero_m_array (expo_b, expo_d) (expo_inv_p, expo_inv_q) - (shell_p.ContractedShellPair.center_ab, - shell_q.ContractedShellPair.center_ab, center_pq) + (shell_p.Csp.center_ab, + shell_q.Csp.center_ab, center_pq) (center_pa, center_qc) map_1d map_2d np nq in @@ -865,8 +867,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (** Computes all the two-electron integrals of the contracted shell quartet *) let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = - let shell_p = ContractedShellPair.create ~cutoff shell_a shell_b - and shell_q = ContractedShellPair.create ~cutoff shell_c shell_d + let shell_p = Csp.create ~cutoff shell_a shell_b + and shell_q = Csp.create ~cutoff shell_c shell_d in contracted_class_shell_pairs ~zero_m shell_p shell_q diff --git a/Nuclei/Element.ml b/Nuclei/Element.ml index 82bb2f6..159671b 100644 --- a/Nuclei/Element.ml +++ b/Nuclei/Element.ml @@ -140,7 +140,7 @@ let covalent_radius x = | Cd -> 1.44 | In -> 1.42 | Sn -> 1.39 | Sb -> 1.39 | Te -> 1.38 | I -> 1.39 | Xe -> 1.40 | Pt -> 1.30 in - Units.angstrom_to_bohr *. (result x) + Constants.a0 *. (result x) |> Radius.of_float @@ -161,7 +161,7 @@ let vdw_radius x = | Cd -> 1.58 | In -> 1.93 | Sn -> 2.17 | Sb -> 2.06 | Te -> 2.06 | I -> 1.98 | Xe -> 2.16 | Pt -> 1.75 in - Units.angstrom_to_bohr *. (result x) + Constants.a0 *. (result x) |> Radius.of_float diff --git a/Nuclei/Mass.ml b/Nuclei/Mass.ml index 0c9e0cc..035c570 100644 --- a/Nuclei/Mass.ml +++ b/Nuclei/Mass.ml @@ -1 +1 @@ -include Positive_float +include PositiveFloat diff --git a/Utils/Angular_momentum.ml b/Utils/Angular_momentum.ml deleted file mode 100644 index d7877a2..0000000 --- a/Utils/Angular_momentum.ml +++ /dev/null @@ -1,134 +0,0 @@ -open Powers - -exception AngularMomentumError of string - -type t = - | S | P | D | F | G | H | I | J | K | L | M | N | O - -let of_char = function - | 's' | 'S' -> S | 'p' | 'P' -> P - | 'd' | 'D' -> D | 'f' | 'F' -> F - | 'g' | 'G' -> G | 'h' | 'H' -> H - | 'i' | 'I' -> I | 'j' | 'J' -> J - | 'k' | 'K' -> K | 'l' | 'L' -> L - | 'm' | 'M' -> M | 'n' | 'N' -> N - | 'o' | 'O' -> O - | c -> raise (AngularMomentumError (String.make 1 c)) - -let to_string = function - | S -> "S" | P -> "P" - | D -> "D" | F -> "F" - | G -> "G" | H -> "H" - | I -> "I" | J -> "J" - | K -> "K" | L -> "L" - | M -> "M" | N -> "N" - | O -> "O" - -let to_char = function - | S -> 'S' | P -> 'P' - | D -> 'D' | F -> 'F' - | G -> 'G' | H -> 'H' - | I -> 'I' | J -> 'J' - | K -> 'K' | L -> 'L' - | M -> 'M' | N -> 'N' - | O -> 'O' - -let to_int = function - | S -> 0 | P -> 1 - | D -> 2 | F -> 3 - | G -> 4 | H -> 5 - | I -> 6 | J -> 7 - | K -> 8 | L -> 9 - | M -> 10 | N -> 11 - | O -> 12 - -let of_int = function - | 0 -> S | 1 -> P - | 2 -> D | 3 -> F - | 4 -> G | 5 -> H - | 6 -> I | 7 -> J - | 8 -> K | 9 -> L - | 10 -> M | 11 -> N - | 12 -> O - | c -> raise (AngularMomentumError (string_of_int c)) - - - -type kind = -| Singlet of t -| Doublet of (t*t) -| Triplet of (t*t*t) -| Quartet of (t*t*t*t) - - -let n_functions a = - let a = - to_int a - in - (a*a + 3*a + 2)/2 - - -(** Returns an array of Zkeys corresponding to all possible angular momenta *) -let zkey_array a = - let keys_1d l = - let create_z { x ; y ; _ } = - Powers.of_int_tuple (x,y,l-(x+y)) - in - let rec create_y accu xyz = - let { x ; y ; z } = xyz in - match y with - | 0 -> (create_z xyz)::accu - | i -> let ynew = y-1 in - create_y ( (create_z xyz)::accu) (Powers.of_int_tuple (x,ynew,z)) - in - let rec create_x accu xyz = - let { x ; y ; z } = xyz in - match x with - | 0 -> (create_y [] xyz)@accu - | i -> let xnew = x-1 in - let ynew = l-xnew in - create_x ((create_y [] xyz)@accu) (Powers.of_int_tuple (xnew, ynew, z)) - in - create_x [] (Powers.of_int_tuple (l,0,0)) - |> List.rev - in - - begin - match a with - | Singlet l1 -> - List.map (fun x -> Zkey.of_powers (Zkey.Three x)) (keys_1d @@ to_int l1) - - | Doublet (l1, l2) -> - List.map (fun a -> - List.map (fun b -> - Zkey.of_powers (Zkey.Six (a,b))) (keys_1d @@ to_int l2) - ) (keys_1d @@ to_int l1) - |> List.concat - - | Triplet (l1, l2, l3) -> - - List.map (fun a -> - List.map (fun b -> - List.map (fun c -> - Zkey.of_powers (Zkey.Nine (a,b,c))) (keys_1d @@ to_int l3) - ) (keys_1d @@ to_int l2) - |> List.concat - ) (keys_1d @@ to_int l1) - |> List.concat - - | Quartet (l1, l2, l3, l4) -> - - List.map (fun a -> - List.map (fun b -> - List.map (fun c -> - List.map (fun d -> - Zkey.of_powers (Zkey.Twelve (a,b,c,d))) (keys_1d @@ to_int l4) - ) (keys_1d @@ to_int l3) - |> List.concat - ) (keys_1d @@ to_int l2) - |> List.concat - ) (keys_1d @@ to_int l1) - |> List.concat - end - |> Array.of_list - diff --git a/Utils/Positive_float.ml b/Utils/Positive_float.ml deleted file mode 100644 index fbf46c0..0000000 --- a/Utils/Positive_float.ml +++ /dev/null @@ -1,14 +0,0 @@ -type t = float - -let of_float x = - assert ( x >= 0. ); - x - -external to_float : t -> float = "%identity" - -let to_string x = - let f = to_float x in string_of_float f - -let of_string x = - let f = float_of_string x in of_float f - diff --git a/Utils/Positive_float.mli b/Utils/Positive_float.mli deleted file mode 100644 index 2d3d60b..0000000 --- a/Utils/Positive_float.mli +++ /dev/null @@ -1,5 +0,0 @@ -type t = private float -val of_float : float -> t -val to_float : t -> float -val to_string : t -> string -val of_string : string -> t diff --git a/Utils/Radius.ml b/Utils/Radius.ml index 0c9e0cc..035c570 100644 --- a/Utils/Radius.ml +++ b/Utils/Radius.ml @@ -1 +1 @@ -include Positive_float +include PositiveFloat diff --git a/Utils/Units.ml b/Utils/Units.ml deleted file mode 100644 index 0ec8030..0000000 --- a/Utils/Units.ml +++ /dev/null @@ -1,25 +0,0 @@ -type units = -| Bohr -| Angstrom - -type angle_units = -| Degree -| Radian - -let pi = acos (-1.) - -let to_degree x = - assert (x <= 2.*.pi); - assert (x >= -2.*.pi); - x *. 180. /. pi - -let to_radian x = - assert (x <= 360.); - assert (x >= -360.); - x *. pi /. 180. - -let angstrom_to_bohr = 1. /. 0.52917721092 -let bohr_to_angstrom = 0.52917721092 -;; - - diff --git a/Utils/Units.mli b/Utils/Units.mli deleted file mode 100644 index 6c01628..0000000 --- a/Utils/Units.mli +++ /dev/null @@ -1,8 +0,0 @@ -type units = Bohr | Angstrom -type angle_units = Degree | Radian - -val to_radian : float -> float -val to_degree : float -> float - -val angstrom_to_bohr : float -val bohr_to_angstrom : float