diff --git a/GoDuck b/GoDuck index 13eb269..06927ce 100755 --- a/GoDuck +++ b/GoDuck @@ -1,18 +1,68 @@ #! /bin/bash + +function help() { + cat << EOF +Usage: + ./qcaml-tools/quack_input -b [-c ] [-f] [-h] [-m ] + [-r ] -x [--] + +Arguments: + + +Options: + + -b --basis= Name of the file containing the basis set + -c --charge= Total charge of the molecule. Specify negative + charges with 'm' instead of the minus sign, + for example m1 instead of -1. Default is + 0 + -f --frozen-core Freeze core MOs. Default is false + -h --help Prints the help message. + -m --multiplicity= Spin multiplicity (2S+1). Default is singlet + -r --rydberg= Number of Rydberg electrons. Default is 0 + -x --xyz= Name of the file containing the nuclear + coordinates in xyz format + +Description: + + 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. + + +Usage: + + $(basename $0) [-h|--help] + +Options: + + -h --help Prints the help message + +EOF + exit 0 +} + set -e -if [ $# -ne 2 ] +if [ $# -ne 2 ] || [ $# -ne 4 ] then - echo "You need two arguments [Molecule] [Basis] !!" + echo "You need at least two arguments [Molecule] [Basis] [Charge] [Multiplicity]" fi if [ $# = 2 ] then - cp examples/molecule."$1" input/molecule - cp examples/basis."$1"."$2" input/basis - cp basis/"$2" input/basis.qcaml -# ./bin/IntPak + ./qcaml-tools/quack_integrals -x mol/${1}.xyz -b basis/${2} + ./qcaml-tools/quack_input -x mol/${1}.xyz -b basis/${2} ./bin/QuAcK fi +if [ $# = 4 ] +then + ./qcaml-tools/quack_integrals -x mol/${1}.xyz -b basis/${2} + ./qcaml-tools/quack_input -x mol/${1}.xyz -b basis/${2} -c ${3} -m ${4} + ./bin/QuAcK +fi + + diff --git a/qcaml-tools/GoDuck b/qcaml-tools/GoDuck new file mode 100755 index 0000000..776b1c7 Binary files /dev/null and b/qcaml-tools/GoDuck differ diff --git a/qcaml-tools/GoDuck.ml b/qcaml-tools/GoDuck.ml new file mode 100644 index 0000000..4a2b481 --- /dev/null +++ b/qcaml-tools/GoDuck.ml @@ -0,0 +1,185 @@ +let quack_dir = + try Sys.getenv "QUACK_ROOT" with + Not_found -> "." + +let quack_input = quack_dir ^ "/input/" +let quack_mol = quack_dir ^ "/mol/" +let quack_basis = quack_dir ^ "/basis/" +let quack_int = quack_dir ^ "/int/" + +let quack_basis_filename = quack_input ^ "basis" +let quack_molecule_filename = quack_input ^ "molecule" + +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 "Prepares the input data for QuAcK. +Writes $QUACK_ROOT/input/basis and $QUACK_ROOT/input/molecule. +If $QUACK_ROOT is not set, $QUACK_ROOT 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 in the $QUACK_ROOT/basis/ directory"; } ; + + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format in the $QUACK_ROOT/mol/ directory without the .xyz extension"; } ; + + { 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"; } ; + + { 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"; } ; + + { short='u' ; long="range-separation" ; opt=Optional; + arg=With_arg ""; + doc="Range-separation parameter."; } ; + ] + end; + + (* Handle options *) + let basis_file = + quack_basis ^ (Util.of_some @@ Command_line.get "basis") + in + let nuclei_file = + quack_mol ^ (Util.of_some @@ Command_line.get "xyz") ^ ".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 + + let multiplicity = + match Command_line.get "multiplicity" with + | Some x -> int_of_string x + | None -> 1 + in + + let rydberg = + match Command_line.get "rydberg" with + | Some x -> int_of_string x + | None -> 0 + in + + let range_separation = + match Command_line.get "range-separation" with + | None -> None + | Some mu -> Some (float_of_string mu) + in + + + let nuclei = + Qcaml.Particles.Nuclei.of_xyz_file nuclei_file + in + + let electrons = + Qcaml.Particles.Electrons.of_atoms ~multiplicity ~charge nuclei + in + + let basis = + Qcaml.Gaussian.Basis.of_nuclei_and_basis_filename ~nuclei basis_file + in + + + 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 operators = + match range_separation with + | None -> [] + | Some mu -> [ Qcaml.Operators.Operator.of_range_separation mu ] + in + + let ao_basis = + Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian + ~operators ~cartesian:true ~nuclei basis_file + in + + let overlap = Qcaml.Ao.Basis.overlap ao_basis in + let eN_ints = Qcaml.Ao.Basis.eN_ints ao_basis in + let kin_ints = Qcaml.Ao.Basis.kin_ints ao_basis in + let ee_ints = Qcaml.Ao.Basis.ee_ints ao_basis in + let multipole = Qcaml.Ao.Basis.multipole ao_basis in + let x_mat = Qcaml.Gaussian_integrals.Multipole.matrix_x multipole in + let y_mat = Qcaml.Gaussian_integrals.Multipole.matrix_y multipole in + let z_mat = Qcaml.Gaussian_integrals.Multipole.matrix_z multipole in + + Qcaml.Gaussian_integrals.Overlap.to_file ~filename:(quack_int ^ "Ov.dat") overlap; + Qcaml.Gaussian_integrals.Electron_nucleus.to_file ~filename:(quack_int ^ "Nuc.dat") eN_ints; + Qcaml.Gaussian_integrals.Kinetic.to_file ~filename:(quack_int ^ "Kin.dat") kin_ints; + Qcaml.Gaussian_integrals.Eri.to_file ~filename:(quack_int ^ "ERI.dat") ee_ints; + Qcaml.Gaussian_integrals.Multipole.to_file ~filename:(quack_int ^ "x.dat") x_mat; + Qcaml.Gaussian_integrals.Multipole.to_file ~filename:(quack_int ^ "y.dat") y_mat; + Qcaml.Gaussian_integrals.Multipole.to_file ~filename:(quack_int ^ "z.dat") z_mat; + + match range_separation with + | Some _mu -> + Qcaml.Gaussian_integrals.Eri_long_range.to_file ~filename:(quack_int ^ "ERI_lr.dat") (Qcaml.Ao.Basis.ee_lr_ints ao_basis) + | None -> () + + + + + diff --git a/qcaml-tools/Makefile b/qcaml-tools/Makefile index 9ee8806..f3dd1ee 100644 --- a/qcaml-tools/Makefile +++ b/qcaml-tools/Makefile @@ -2,7 +2,8 @@ .NOTPARALLEL: TARGETS=quack_input \ - quack_integrals + quack_integrals \ + GoDuck .PHONY: default build install uninstall test clean diff --git a/qcaml-tools/dune b/qcaml-tools/dune index 7bc46c7..c9c995e 100644 --- a/qcaml-tools/dune +++ b/qcaml-tools/dune @@ -2,6 +2,7 @@ (names quack_input quack_integrals + GoDuck ) (libraries qcaml diff --git a/qcaml-tools/quack_input.ml b/qcaml-tools/quack_input.ml index 931e17e..2af2b61 100644 --- a/qcaml-tools/quack_input.ml +++ b/qcaml-tools/quack_input.ml @@ -1,5 +1,5 @@ let quack_dir = - try Sys.getenv "QUACK_DIR" with + try Sys.getenv "QUACK_ROOT" with Not_found -> "." let quack_basis_filename = quack_dir ^ "/input/basis" @@ -13,8 +13,8 @@ let () = 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 +Writes $QUACK_ROOT/input/basis and $QUACK_ROOT/input/molecule. +If $QUACK_ROOT is not set, $QUACK_ROOT is replaces by the current directory."; set_specs [ { short='b' ; long="basis" ; opt=Mandatory; diff --git a/qcaml-tools/quack_integrals.ml b/qcaml-tools/quack_integrals.ml index 14a6827..71a39fa 100644 --- a/qcaml-tools/quack_integrals.ml +++ b/qcaml-tools/quack_integrals.ml @@ -1,5 +1,5 @@ let quack_dir = - try Sys.getenv "QUACK_DIR" with + try Sys.getenv "QUACK_ROOT" with Not_found -> "." module Command_line = Qcaml.Common.Command_line diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index a7a2f0f..4f0b700 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -245,7 +245,7 @@ program QuAcK else - call system('./GoQCaml') +! call system('./GoQCaml') call read_integrals(nBas,S,T,V,Hc,ERI_AO) call read_dipole_integrals(nBas,dipole_int)