diff --git a/qcaml-tools/dune b/qcaml-tools/dune index a2a77f8..7bc46c7 100644 --- a/qcaml-tools/dune +++ b/qcaml-tools/dune @@ -4,7 +4,6 @@ quack_integrals ) (libraries - qcaml.common qcaml ) ) diff --git a/qcaml-tools/quack_input.ml b/qcaml-tools/quack_input.ml index bee8f62..931e17e 100644 --- a/qcaml-tools/quack_input.ml +++ b/qcaml-tools/quack_input.ml @@ -5,84 +5,71 @@ let quack_dir = let quack_basis_filename = quack_dir ^ "/input/basis" let quack_molecule_filename = quack_dir ^ "/input/molecule" +module Command_line = Qcaml.Common.Command_line +module Util = Qcaml.Common.Util -let basis_file : string option ref = ref None -let nuclei_file : string option ref = ref None -let charge : int option ref = ref None -let multiplicity : int option ref = ref None +let () = + let open Command_line in + begin + set_header_doc (Sys.argv.(0) ^ " - QuAcK command"); + set_description_doc "Prepares the input data for QuAcK. +Writes $QUACK_DIR/input/basis and $QUACK_DIR/input/molecule. +If $QUACK_DIR is not set, $QUACK_DIR is replaces by the current +directory."; + set_specs + [ { short='b' ; long="basis" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the basis set"; } ; + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format"; } ; -let speclist = [ - ( "-b" , Arg.String (fun x -> basis_file := Some x), - "File containing the atomic basis set") ; - ( "-c" , Arg.Int (fun x -> charge := Some x), - "Total charge of the system") ; - ( "-m" , Arg.Int (fun x -> multiplicity := Some x), - "Multiplicity of the system") ; - ( "-x" , Arg.String (fun x -> nuclei_file := Some x), - "File containing the nuclear coordinates") ; -] + { short='m' ; long="multiplicity" ; opt=Optional; + arg=With_arg ""; + doc="Spin multiplicity (2S+1). Default is singlet"; } ; + { short='c' ; long="charge" ; opt=Optional; + arg=With_arg ""; + doc="Total charge of the molecule. Specify negative charges with 'm' instead of the minus sign, for example m1 instead of -1. Default is 0"; } ; -let print_basis nuclei basis = - let oc = open_out quack_basis_filename in - let ocf = Format.formatter_of_out_channel oc in - let g_basis = Qcaml.Gaussian.Basis.general_basis basis in - let dict = - Array.to_list nuclei - |> List.mapi (fun i (e, _) -> - (i+1), List.assoc e g_basis - ) + { short='f' ; long="frozen-core" ; opt=Optional; + arg=Without_arg ; + doc="Freeze core MOs. Default is false"; } ; + + { short='r' ; long="rydberg" ; opt=Optional; + arg=With_arg "" ; + doc="Number of Rydberg electrons. Default is 0"; } ; + ] + end; + + (* Handle options *) + let basis_file = Util.of_some @@ Command_line.get "basis" in + let nuclei_file = Util.of_some @@ Command_line.get "xyz" in + let frozen_core = Command_line.get_bool "frozen-core" in + + let charge = + match Command_line.get "charge" with + | Some x -> ( if x.[0] = 'm' then + ~- (int_of_string (String.sub x 1 (String.length x - 1))) + else + int_of_string x ) + | None -> 0 in - List.iter (fun (i,b) -> - Format.fprintf ocf "%d %d\n" i (Array.length b); - Array.iter (fun x -> - Format.fprintf ocf "%a" Qcaml.Gaussian.General_basis.pp_gcs x) b - ) dict; - close_out oc - - -let print_molecule nuclei electrons = - let oc = open_out quack_molecule_filename in - let nat = Array.length nuclei in - let nela = Qcaml.Particles.Electrons.n_alfa electrons in - let nelb = Qcaml.Particles.Electrons.n_beta electrons in - let ncore = Qcaml.Particles.Nuclei.small_core nuclei in - let nryd = 0 in - Printf.fprintf oc "# nAt nEla nElb nCore nRyd\n"; - Printf.fprintf oc " %4d %4d %4d %5d %4d\n" nat nela nelb ncore nryd; - Printf.fprintf oc "# Znuc x y z\n"; - let open Qcaml.Common.Coordinate in - Array.iter (fun (element, coord) -> - Printf.fprintf oc "%3s %16.10f %16.10f %16.10f\n" - (Qcaml.Particles.Element.to_string element) - coord.x coord.y coord.z - ) nuclei; - close_out oc - - - - -let run () = - let basis_file = - match !basis_file with - | None -> raise (Invalid_argument "Basis set file should be specified with -b") - | Some x -> x - and nuclei_file = - match !nuclei_file with - | None -> raise (Invalid_argument "Coordinate file should be specified with -x") - | Some x -> x - and charge = - match !charge with - | None -> 0 - | Some c -> c - and multiplicity = - match !multiplicity with + let multiplicity = + match Command_line.get "multiplicity" with + | Some x -> int_of_string x | None -> 1 - | Some m -> m in + let rydberg = + match Command_line.get "rydberg" with + | Some x -> int_of_string x + | None -> 0 + in + + let nuclei = Qcaml.Particles.Nuclei.of_xyz_file nuclei_file in @@ -95,14 +82,52 @@ let run () = Qcaml.Gaussian.Basis.of_nuclei_and_basis_filename ~nuclei basis_file in - (* Print basis *) - print_molecule nuclei electrons; + + let print_basis nuclei basis = + let oc = open_out quack_basis_filename in + let ocf = Format.formatter_of_out_channel oc in + let g_basis = Qcaml.Gaussian.Basis.general_basis basis in + let dict = + Array.to_list nuclei + |> List.mapi (fun i (e, _) -> + (i+1), List.assoc e g_basis + ) + in + List.iter (fun (i,b) -> + Format.fprintf ocf "%d %d\n" i (Array.length b); + Array.iter (fun x -> + Format.fprintf ocf "%a" Qcaml.Gaussian.General_basis.pp_gcs x) b + ) dict; + close_out oc + in print_basis nuclei basis; + + + let print_molecule nuclei electrons = + let oc = open_out quack_molecule_filename in + let nat = Array.length nuclei in + let nela = Qcaml.Particles.Electrons.n_alfa electrons in + let nelb = Qcaml.Particles.Electrons.n_beta electrons in + let ncore = + if frozen_core then + Qcaml.Particles.Nuclei.small_core nuclei + else 0 + in + let nryd = rydberg in + Printf.fprintf oc "# nAt nEla nElb nCore nRyd\n"; + Printf.fprintf oc " %4d %4d %4d %5d %4d\n" nat nela nelb ncore nryd; + Printf.fprintf oc "# Znuc x y z\n"; + let open Qcaml.Common.Coordinate in + Array.iter (fun (element, coord) -> + Printf.fprintf oc "%3s %16.10f %16.10f %16.10f\n" + (Qcaml.Particles.Element.to_string element) + coord.x coord.y coord.z + ) nuclei; + close_out oc + in + print_molecule nuclei electrons; + () -let () = - let usage_msg = "Available options:" in - Arg.parse speclist (fun _ -> ()) usage_msg; - run () diff --git a/qcaml-tools/quack_integrals.ml b/qcaml-tools/quack_integrals.ml index 0aae58d..37bcf9f 100644 --- a/qcaml-tools/quack_integrals.ml +++ b/qcaml-tools/quack_integrals.ml @@ -1,30 +1,32 @@ -let out_file : string option ref = ref None -let basis_file : string option ref = ref None -let nuclei_file : string option ref = ref None -let charge : int option ref = ref None -let multiplicity : int option ref = ref None -let range_separation : float option ref = ref None +module Command_line = Qcaml.Common.Command_line +module Util = Qcaml.Common.Util +let () = + let open Command_line in + begin + set_header_doc (Sys.argv.(0) ^ " - QuAcK command"); + set_description_doc "Computes the one- and two-electron integrals on the Gaussian atomic basis set."; + set_specs + [ { short='b' ; long="basis" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the basis set"; } ; -let speclist = [ - ( "-b" , Arg.String (fun x -> basis_file := Some x), - "File containing the atomic basis set") ; - ( "-x" , Arg.String (fun x -> nuclei_file := Some x), - "File containing the nuclear coordinates") ; - ( "-u" , Arg.Float (fun x -> range_separation := Some x), - "Value of mu, the range separation factor") ; -] + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format"; } ; -let run () = - let basis_file = - match !basis_file with - | None -> raise (Invalid_argument "Basis set file should be specified with -b") - | Some x -> x - and nuclei_file = - match !nuclei_file with - | None -> raise (Invalid_argument "Coordinate file should be specified with -x") - | Some x -> x - and range_separation = !range_separation + { short='u' ; long="range-separation" ; opt=Optional; + arg=With_arg ""; + doc="Range-separation parameter."; } ; + ] + end; + + let basis_file = Util.of_some @@ Command_line.get "basis" in + let nuclei_file = Util.of_some @@ Command_line.get "xyz" in + let range_separation = + match Command_line.get "range-separation" with + | None -> None + | Some mu -> Some (float_of_string mu) in let nuclei = @@ -56,8 +58,4 @@ let run () = | None -> () -let () = - let usage_msg = "Available options:" in - Arg.parse speclist (fun _ -> ()) usage_msg; - run ()