diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 252b96b..e7e8eaf 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -1,41 +1,20 @@ -(** General basis set read from a file *) -type primitive = { - exponent: float ; - coefficient: float - } +type t -type element_name = string - -type general_contracted_shell = Angular_momentum.t * (primitive array) - -type basis_set = element_name * (general_contracted_shell array) +let of_nuclei_and_general_basis n b = + Array.map (fun (e, center) -> + List.assoc e b + |> Array.map (fun (totAngMom, shell) -> + let expo = Array.map (fun General_basis.{exponent ; coefficient} -> + exponent) shell + and coef = Array.map (fun General_basis.{exponent ; coefficient} -> + coefficient) shell + in + Contracted_shell.create ~expo ~coef ~totAngMom ~center) + ) n -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" - (Angular_momentum.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 +let to_string b = + Array.map (fun i -> Contracted_shell.to_string i) b |> 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/Contracted_shell.ml b/Basis/Contracted_shell.ml index 5bc414c..2c4e195 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -17,6 +17,14 @@ let totAngMom a = a.totAngMom let norm_coef a i = a.norm_coef.(i) +let to_string s = + let open Printf in + [ sprintf "center: %s" (Coordinate.to_string s.center) ; + sprintf "angular momentum: %s" (Angular_momentum.to_string s.totAngMom) ] + @ (Array.map2 (fun e c -> sprintf "expo: %e coeff: %e" e c) s.expo s.coef + |> Array.to_list) @ ["\n"] + |> String.concat "\n" + (** 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 diff --git a/Basis/ERI.ml b/Basis/ERI.ml index d27432e..f88e5f0 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -210,8 +210,8 @@ let erint_contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = let class_indices = Angular_momentum.zkey_array (Angular_momentum.Kind_4 - (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, - Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)) + Contracted_shell.(totAngMom shell_a, totAngMom shell_b, + totAngMom shell_c, totAngMom shell_d)) in let contracted_class = diff --git a/Basis/Gamess_parser.mly b/Basis/Gamess_parser.mly index ae84ac0..e15978a 100644 --- a/Basis/Gamess_parser.mly +++ b/Basis/Gamess_parser.mly @@ -13,7 +13,7 @@ %token EOF %start input -%type input +%type input %% /* Grammar rules and actions follow */ @@ -43,7 +43,7 @@ primitive_list: | primitive_list primitive { $2 :: $1 } primitive: - | INTEGER FLOAT FLOAT EOL { Basis.{exponent=$2 ; coefficient=$3 } } + | INTEGER FLOAT FLOAT EOL { General_basis.{exponent=$2 ; coefficient=$3 } } diff --git a/Basis/General_basis.ml b/Basis/General_basis.ml new file mode 100644 index 0000000..7ea7642 --- /dev/null +++ b/Basis/General_basis.ml @@ -0,0 +1,40 @@ +(** General basis set read from a file *) +type primitive = { + exponent: float ; + coefficient: float + } + +type element_name = string + +type general_contracted_shell = Angular_momentum.t * (primitive array) + +type t = element_name * (general_contracted_shell array) + + +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" + (Angular_momentum.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/Shell_pair.ml b/Basis/Shell_pair.ml index 87d97c6..af46c61 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/Shell_pair.ml @@ -21,15 +21,15 @@ let create_array ?(cutoff=0.) p_a p_b = else -. (log cutoff) in - let center_ab = - Coordinate.(Contracted_shell.center p_a |- Contracted_shell.center p_b) + let center_ab = Coordinate.( + Contracted_shell.center p_a |- Contracted_shell.center p_b ) in let norm_sq = Coordinate.dot center_ab center_ab in Array.init (Contracted_shell.size p_a) (fun i -> - let p_a_expo_center = - Coordinate.(Contracted_shell.expo p_a i |. Contracted_shell.center p_a) + let p_a_expo_center = Coordinate.( + Contracted_shell.expo p_a i |. Contracted_shell.center p_a ) in let f1 = Contracted_shell.norm_coef p_a i @@ -50,13 +50,13 @@ let create_array ?(cutoff=0.) p_a p_b = in if (norm < cutoff) then raise Null_contribution; - let p_b_expo_center = - Coordinate.(Contracted_shell.expo p_b j |. Contracted_shell.center p_b) + let p_b_expo_center = Coordinate.( + Contracted_shell.expo p_b j |. Contracted_shell.center p_b ) in let expo = Contracted_shell.(expo p_a i +. expo p_b j) in let expo_inv = 1. /. expo in - let center = - Coordinate.( expo_inv |. (p_a_expo_center |+ p_b_expo_center ) ) + let center = Coordinate.( + expo_inv |. (p_a_expo_center |+ p_b_expo_center ) ) in let argexpo = Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv diff --git a/Makefile b/Makefile index e06bb59..7b5d1c0 100644 --- a/Makefile +++ b/Makefile @@ -9,12 +9,11 @@ MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) MLYFILES=$(wildcard */*.mly) $(wildcard *.mly) MLFILES= $(wildcard */*.ml) $(wildcard *.ml) MLIFILES=$(wildcard */*.mli) $(wildcard *.mli) -ALL_EXE=$(patsubst %.ml,%.native,$(wildcard qp_*.ml)) +ALL_EXE=$(patsubst %.ml,%.native,$(wildcard run_*.ml)) .PHONY: default -#default: $(ALL_EXE) -default: test.byte +default: $(ALL_EXE) tests: $(ALL_TESTS) diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index a4b528c..f0834c1 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -1,4 +1,4 @@ -type coordinates = ( (Element.t * (float array)) array) +type coordinates = ( (Element.t * Coordinate.t) array) let of_xyz_file ~filename = let lexbuf = @@ -21,5 +21,5 @@ let of_zmt_file ~filename = in aux [] |> Zmatrix.of_string |> Zmatrix.to_xyz - |> Array.map (fun (e,x,y,z) -> (e, [|x;y;z|])) + |> Array.map (fun (e,x,y,z) -> (e, Coordinate.of_3_floats x y z )) diff --git a/Nuclei/Xyz_parser.mly b/Nuclei/Xyz_parser.mly index 19c9c4b..4785570 100644 --- a/Nuclei/Xyz_parser.mly +++ b/Nuclei/Xyz_parser.mly @@ -12,7 +12,7 @@ exception InputError of string %token EOF %start input -%type <(Element.t * (float array)) array> input +%type <(Element.t * Coordinate.t) array> input %% /* Grammar rules and actions follow */ @@ -49,10 +49,10 @@ atoms_xyz: 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 } - | 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 } + | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_string $2, Coordinate.of_3_floats $4 $6 $8 ) :: $1 } + | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_string $2, Coordinate.of_3_floats $4 $6 $8 ) :: $1 } + | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_int $2, Coordinate.of_3_floats $4 $6 $8 ) :: $1 } + | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_int $2, Coordinate.of_3_floats $4 $6 $8 ) :: $1 } diff --git a/run_integrals.ml b/run_integrals.ml new file mode 100644 index 0000000..7328191 --- /dev/null +++ b/run_integrals.ml @@ -0,0 +1,37 @@ +let basis_file : string option ref = ref None +let coord_file : string option ref = ref None + + +let speclist = [ +( "-b" , Arg.String (fun x -> basis_file := Some x) , + "File containing the atomic basis set.") ; +( "-c" , Arg.String (fun x -> coord_file := Some x) , + "File containing the nuclear coordinates.") ; +] + +let run ~coord ~basis = + let coord_file = + match coord with + | None -> raise (Invalid_argument "Coordinate file should be specified") + | Some x -> x + and basis_file = + match basis with + | None -> raise (Invalid_argument "Basis set file should be specified") + | Some x -> x + in + let nuclei = + Nuclei.of_xyz_file ~filename:coord_file + and general_basis = + Gamess_reader.read ~filename:basis_file + in + let basis = + Basis.of_nuclei_and_general_basis nuclei general_basis + in + Array.map Basis.to_string basis + |> Array.iter print_endline + + +let () = + let usage_msg = "Quantum Chemistry Package. Options available:" in + Arg.parse speclist (fun _ -> ()) usage_msg; + run ~coord:!coord_file ~basis:!basis_file