diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 861fa29..6e8e131 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -1,11 +1,15 @@ +(** General basis set read from a file *) + type primitive = { exponent: float ; coefficient: float } -type contracted_shell = Angular_momentum.t * (primitive array) +type element_name = string -type gaussian_basis_set = string * (contracted_shell array) +type general_contracted_shell = Angular_momentum.t * (primitive array) + +type basis_set = element_name * (general_contracted_shell array) let string_of_primitive ?id prim = @@ -30,6 +34,9 @@ let string_of_contracted_shell_array 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/Basis_lexer.mll b/Basis/Basis_lexer.mll index 5ac3885..2540494 100644 --- a/Basis/Basis_lexer.mll +++ b/Basis/Basis_lexer.mll @@ -1,7 +1,7 @@ { exception SyntaxError of string -open Basis_parser +open Gamess_parser } diff --git a/Basis/Basis_parser.mly b/Basis/Gamess_parser.mly similarity index 90% rename from Basis/Basis_parser.mly rename to Basis/Gamess_parser.mly index 61d4f00..ae84ac0 100644 --- a/Basis/Basis_parser.mly +++ b/Basis/Gamess_parser.mly @@ -1,4 +1,4 @@ -/* Parses basis sets in GAMESS format */ +/* Parses basis sets GAMESS format */ %{ @@ -13,7 +13,7 @@ %token EOF %start input -%type input +%type input %% /* Grammar rules and actions follow */ @@ -46,3 +46,4 @@ primitive: | INTEGER FLOAT FLOAT EOL { Basis.{exponent=$2 ; coefficient=$3 } } + diff --git a/Basis/Gamess_reader.ml b/Basis/Gamess_reader.ml new file mode 100644 index 0000000..f6f69c5 --- /dev/null +++ b/Basis/Gamess_reader.ml @@ -0,0 +1,21 @@ +(** 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 element, basis = + Gamess_parser.input Basis_lexer.read_all lexbuf + in + let key = + Element.of_string element + in + aux ((key, basis)::accu) + with + | Parsing.Parse_error -> List.rev accu + in + aux [] + diff --git a/Makefile b/Makefile index 50fbdc4..e06bb59 100644 --- a/Makefile +++ b/Makefile @@ -13,14 +13,16 @@ ALL_EXE=$(patsubst %.ml,%.native,$(wildcard qp_*.ml)) .PHONY: default -default: $(ALL_EXE) +#default: $(ALL_EXE) +default: test.byte + tests: $(ALL_TESTS) qpackage.odocl: $(MLIFILES) ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl doc: qpackage.odocl - $(OCAMLBUILD) qpackage.docdir/index.html -use-ocamlfind $(PKGS) + $(OCAMLBUILD) qpackage.docdir/index.html -use-ocamlfind $(PKGS) %.inferred.mli: $(MLFILES) $(OCAMLBUILD) $*.inferred.mli -use-ocamlfind $(PKGS) @@ -28,7 +30,7 @@ doc: qpackage.odocl %.byte: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) rm -f -- $* - $(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS) + $(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS) ln -s $*.byte $* %.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) diff --git a/Nuclei/Element.ml b/Nuclei/Element.ml index 248731c..82bb2f6 100644 --- a/Nuclei/Element.ml +++ b/Nuclei/Element.ml @@ -79,44 +79,48 @@ let to_long_string = function | Sb -> "Antimony" | Te -> "Tellurium" | I -> "Iodine" | Xe -> "Xenon" | Pt -> "Platinum" + +let to_int = function +| X -> 0 | H -> 1 | He -> 2 | Li -> 3 +| Be -> 4 | B -> 5 | C -> 6 | N -> 7 +| O -> 8 | F -> 9 | Ne -> 10 | Na -> 11 +| Mg -> 12 | Al -> 13 | Si -> 14 | P -> 15 +| S -> 16 | Cl -> 17 | Ar -> 18 | K -> 19 +| Ca -> 20 | Sc -> 21 | Ti -> 22 | V -> 23 +| Cr -> 24 | Mn -> 25 | Fe -> 26 | Co -> 27 +| Ni -> 28 | Cu -> 29 | Zn -> 30 | Ga -> 31 +| Ge -> 32 | As -> 33 | Se -> 34 | Br -> 35 +| Kr -> 36 | Rb -> 37 | Sr -> 38 | Y -> 39 +| Zr -> 40 | Nb -> 41 | Mo -> 42 | Tc -> 43 +| Ru -> 44 | Rh -> 45 | Pd -> 46 | Ag -> 47 +| Cd -> 48 | In -> 49 | Sn -> 50 | Sb -> 51 +| Te -> 52 | I -> 53 | Xe -> 54 | Pt -> 78 + + let to_charge c = - begin - match c with - | X -> 0 | H -> 1 | He -> 2 | Li -> 3 - | Be -> 4 | B -> 5 | C -> 6 | N -> 7 - | O -> 8 | F -> 9 | Ne -> 10 | Na -> 11 - | Mg -> 12 | Al -> 13 | Si -> 14 | P -> 15 - | S -> 16 | Cl -> 17 | Ar -> 18 | K -> 19 - | Ca -> 20 | Sc -> 21 | Ti -> 22 | V -> 23 - | Cr -> 24 | Mn -> 25 | Fe -> 26 | Co -> 27 - | Ni -> 28 | Cu -> 29 | Zn -> 30 | Ga -> 31 - | Ge -> 32 | As -> 33 | Se -> 34 | Br -> 35 - | Kr -> 36 | Rb -> 37 | Sr -> 38 | Y -> 39 - | Zr -> 40 | Nb -> 41 | Mo -> 42 | Tc -> 43 - | Ru -> 44 | Rh -> 45 | Pd -> 46 | Ag -> 47 - | Cd -> 48 | In -> 49 | Sn -> 50 | Sb -> 51 - | Te -> 52 | I -> 53 | Xe -> 54 | Pt -> 78 - end - |> Charge.of_int + to_int c |> Charge.of_int + + +let of_int = function +| 0 -> X | 1 -> H | 2 -> He | 3 -> Li +| 4 -> Be | 5 -> B | 6 -> C | 7 -> N +| 8 -> O | 9 -> F | 10 -> Ne | 11 -> Na +| 12 -> Mg | 13 -> Al | 14 -> Si | 15 -> P +| 16 -> S | 17 -> Cl | 18 -> Ar | 19 -> K +| 20 -> Ca | 21 -> Sc | 22 -> Ti | 23 -> V +| 24 -> Cr | 25 -> Mn | 26 -> Fe | 27 -> Co +| 28 -> Ni | 29 -> Cu | 30 -> Zn | 31 -> Ga +| 32 -> Ge | 33 -> As | 34 -> Se | 35 -> Br +| 36 -> Kr | 37 -> Rb | 38 -> Sr | 39 -> Y +| 40 -> Zr | 41 -> Nb | 42 -> Mo | 43 -> Tc +| 44 -> Ru | 45 -> Rh | 46 -> Pd | 47 -> Ag +| 48 -> Cd | 49 -> In | 50 -> Sn | 51 -> Sb +| 52 -> Te | 53 -> I | 54 -> Xe | 78 -> Pt +| x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown")) let of_charge c = - match Charge.to_int c with - | 0 -> X | 1 -> H | 2 -> He | 3 -> Li - | 4 -> Be | 5 -> B | 6 -> C | 7 -> N - | 8 -> O | 9 -> F | 10 -> Ne | 11 -> Na - | 12 -> Mg | 13 -> Al | 14 -> Si | 15 -> P - | 16 -> S | 17 -> Cl | 18 -> Ar | 19 -> K - | 20 -> Ca | 21 -> Sc | 22 -> Ti | 23 -> V - | 24 -> Cr | 25 -> Mn | 26 -> Fe | 27 -> Co - | 28 -> Ni | 29 -> Cu | 30 -> Zn | 31 -> Ga - | 32 -> Ge | 33 -> As | 34 -> Se | 35 -> Br - | 36 -> Kr | 37 -> Rb | 38 -> Sr | 39 -> Y - | 40 -> Zr | 41 -> Nb | 42 -> Mo | 43 -> Tc - | 44 -> Ru | 45 -> Rh | 46 -> Pd | 47 -> Ag - | 48 -> Cd | 49 -> In | 50 -> Sn | 51 -> Sb - | 52 -> Te | 53 -> I | 54 -> Xe | 78 -> Pt - | x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown")) + Charge.to_int c |> of_int let covalent_radius x = diff --git a/Nuclei/Element.mli b/Nuclei/Element.mli index 81e7acf..3c78084 100644 --- a/Nuclei/Element.mli +++ b/Nuclei/Element.mli @@ -1,6 +1,14 @@ exception ElementError of string -type t +type t = +|X +|H |He +|Li|Be |B |C |N |O |F |Ne +|Na|Mg |Al|Si|P |S |Cl|Ar +|K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr +|Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe + |Pt + (** String conversion functions *) val of_string : string -> t @@ -8,6 +16,8 @@ val to_string : t -> string val to_long_string : t -> string (** Properties *) +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 diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index cf030de..a4b528c 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -1 +1,25 @@ -type xyz_input = string * ( (Element.t * (float array)) array) +type coordinates = ( (Element.t * (float array)) array) + +let of_xyz_file ~filename = + let lexbuf = + let ic = open_in filename in + Lexing.from_channel ic + in + Xyz_parser.input Nuclei_lexer.read_all lexbuf + + +let of_zmt_file ~filename = + let ic = open_in filename in + let rec aux accu = + try + let line = input_line ic in + aux (line::accu) + with End_of_file -> + close_in ic; + List.rev accu + |> String.concat "\n" + in aux [] + |> Zmatrix.of_string + |> Zmatrix.to_xyz + |> Array.map (fun (e,x,y,z) -> (e, [|x;y;z|])) + diff --git a/Nuclei/Nuclei_lexer.mll b/Nuclei/Nuclei_lexer.mll index 86efdb3..33b68bb 100644 --- a/Nuclei/Nuclei_lexer.mll +++ b/Nuclei/Nuclei_lexer.mll @@ -1,5 +1,5 @@ { -open Nuclei_parser +open Xyz_parser } let eol = ['\n'] diff --git a/Nuclei/Nuclei_parser.mly b/Nuclei/Xyz_parser.mly similarity index 69% rename from Nuclei/Nuclei_parser.mly rename to Nuclei/Xyz_parser.mly index b43857e..19c9c4b 100644 --- a/Nuclei/Nuclei_parser.mly +++ b/Nuclei/Xyz_parser.mly @@ -1,4 +1,4 @@ -/* Parses basis sets in GAMESS format */ +/* Parses nuclear coordinates in xyz format */ %{ exception InputError of string @@ -11,19 +11,19 @@ exception InputError of string %token FLOAT %token EOF -%start input_xyz -%type input_xyz +%start input +%type <(Element.t * (float array)) array> input %% /* Grammar rules and actions follow */ -input_xyz: +input: | integer title atoms_xyz { let len = List.length $3 in if len <> $1 then let error_msg = Printf.sprintf "%d atoms entered, expected %d" len $1 in raise (InputError error_msg) else - ($2,Array.of_list $3) + Array.of_list $3 } integer: @@ -51,21 +51,8 @@ atoms_list: | { [] } | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_string $2, [| $4;$6;$8 |]) :: $1 } | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_string $2, [| $4;$6;$8 |]) :: $1 } - -label: - | FLOAT { Zmatrix.Value $1 } - | WORD { Zmatrix.Label $1 } - -first_line: - | WORD { $1 } - -second_line: - | WORD INTEGER label { ($1,$2,$3) } - -third_line: - | WORD INTEGER label INTEGER label { ($1,$2,$3,$4,$5) } - -nth_line: - | WORD INTEGER label INTEGER label INTEGER label { ($1,$2,$3,$4,$5,$6,$7) } + | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_int $2, [| $4;$6;$8 |]) :: $1 } + | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_int $2, [| $4;$6;$8 |]) :: $1 } + diff --git a/Nuclei/Zmatrix.ml b/Nuclei/Zmatrix.ml index 3eb4577..cba38d6 100644 --- a/Nuclei/Zmatrix.ml +++ b/Nuclei/Zmatrix.ml @@ -200,20 +200,6 @@ let rotation_matrix axis angle = let apply_rotation_matrix rot u = (dot rot.(0) u, dot rot.(1) u, dot rot.(2) u) -let center_of_mass l = -let (x,y,z) = - let sum_mass, com = - Array.fold_left (fun (s,com) (e,x,y,z) -> - let mass = - Positive_float.to_float @@ Element.mass e - in - (s +. mass, ( mass |. (x,y,z) ) |+ com) ) - (0., (0.,0.,0.)) l - in - (1. /. sum_mass) |. com -in -Printf.printf "%f %f %f\n" x y z ; (x,y,z) - let to_xyz (z,map) = let result = Array.make (Array.length z) None @@ -300,14 +286,16 @@ let to_xyz (z,map) = | Some x -> x | None -> failwith "Some atoms were not defined" ) result in - Array.to_list result + result let to_xyz_string (l,map) = String.concat "\n" ( to_xyz (l,map) - |> List.map (fun (e,x,y,z) -> - Printf.sprintf "%s %f %f %f\n" (Element.to_string e) x y z) ) + |> Array.map (fun (e,x,y,z) -> + Printf.sprintf "%s %f %f %f\n" (Element.to_string e) x y z) + |> Array.to_list + ) diff --git a/Nuclei/Radius.ml b/Utils/Radius.ml similarity index 100% rename from Nuclei/Radius.ml rename to Utils/Radius.ml diff --git a/Utils/Units.ml b/Utils/Units.ml index ae8b21c..0ec8030 100644 --- a/Utils/Units.ml +++ b/Utils/Units.ml @@ -1,7 +1,22 @@ 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 index eff3fb2..6c01628 100644 --- a/Utils/Units.mli +++ b/Utils/Units.mli @@ -1,7 +1,8 @@ -type units = Bohr | Angstrom +type units = Bohr | Angstrom +type angle_units = Degree | Radian + +val to_radian : float -> float +val to_degree : float -> float -(** Conversion functions *) val angstrom_to_bohr : float val bohr_to_angstrom : float - - diff --git a/Utils/Util.ml b/Utils/Util.ml new file mode 100644 index 0000000..6fb0b6f --- /dev/null +++ b/Utils/Util.ml @@ -0,0 +1,73 @@ +(** Constants *) +let pi = acos (-1.) +let pi_inv = 1. /. pi +let two_over_sq_pi = 2. /. (sqrt pi) +let factmax = 150 + + +(** Generalized Boys function. Uses GSL's incomplete Gamma function. + maxm : Maximum total angular momentum +*) +let boys_function ~maxm t = + begin + if t = 0. then + Array.init (maxm+1) (fun m -> 1. /. float_of_int (m+m+1)) + else + let incomplete_gamma ~alpha x = + Gsl.Sf.gamma alpha *. ( Gsl.Sf.gamma_inc_P alpha x ) + in + let t_inv = + 1. /. t + in + let factor = + Array.make (maxm+1) (0.5, sqrt t_inv); + in + for i=1 to maxm + do + let (dm, f) = factor.(i-1) in + factor.(i) <- (dm +. 1., f *. t_inv); + done; + Array.map (fun (dm, f) -> (incomplete_gamma dm t ) *. 0.5 *. f) factor + end + + + +let fact_memo = + let rec aux accu_l accu = function + | 0 -> aux [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 (x::accu_l) x (i+1) + in + aux [] 0. 0 + |> Array.of_list + + + +(** 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") + | i when (i > 150) -> + raise (Invalid_argument "Result of factorial is infinite") + | 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 -> + let b = pow a (n / 2) in + b *. b *. (if n mod 2 = 0 then 1. else a) +;; + + diff --git a/_tags b/_tags index 633df42..e3819e0 100644 --- a/_tags +++ b/_tags @@ -1 +1 @@ -true: package(str) +true: package(str,gsl,zarith)