From 3890394e5f1f7cb7157c32b492e80eef3f37f597 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Sat, 30 Nov 2024 21:10:21 +0100 Subject: [PATCH] Add support for counter-poise calculations (#116) --- README.md | 7 + app/cli.f90 | 152 ++- app/driver.f90 | 245 +++- app/help.f90 | 26 +- doc/_static/references.bib | 87 ++ doc/api/c.rst | 50 +- doc/index.rst | 3 +- fpm.toml | 5 + include/s-dftd3.h | 42 +- src/dftd3/CMakeLists.txt | 2 + src/dftd3/api.f90 | 145 +++ src/dftd3/cutoff.f90 | 12 + src/dftd3/gcp.f90 | 1431 +++++++++++++++++++++++ src/dftd3/gcp/CMakeLists.txt | 24 + src/dftd3/gcp/meson.build | 19 + src/dftd3/gcp/param.f90 | 1084 +++++++++++++++++ src/dftd3/meson.build | 2 + src/dftd3/output.f90 | 80 +- test/api/api-test.c | 214 +++- test/unit/CMakeLists.txt | 1 + test/unit/main.f90 | 2 + test/unit/meson.build | 1 + test/unit/test_gcp.f90 | 1082 +++++++++++++++++ test/validation/21-ad3.coord | 24 + test/validation/21-gcp-pbeh3c.json | 83 ++ test/validation/21-gcp-pbeh3c.resp | 6 + test/validation/22-benI-acetone.xyz | 24 + test/validation/22-gcp-b973c.json | 13 + test/validation/22-gcp-b973c.resp | 6 + test/validation/23-444.sdf | 43 + test/validation/23-gcp-hf3c.json | 4 + test/validation/23-gcp-hf3c.resp | 4 + test/validation/24-gcp-hf-deftzvp.json | 13 + test/validation/24-gcp-hf-deftzvp.resp | 7 + test/validation/24-h2o_6.qchem | 21 + test/validation/25-g3.gen | 26 + test/validation/25-gcp-b3lyp-631gd.json | 89 ++ test/validation/25-gcp-b3lyp-631gd.resp | 8 + test/validation/meson.build | 5 + test/validation/tester.py | 2 +- 40 files changed, 5024 insertions(+), 70 deletions(-) create mode 100644 src/dftd3/gcp.f90 create mode 100644 src/dftd3/gcp/CMakeLists.txt create mode 100644 src/dftd3/gcp/meson.build create mode 100644 src/dftd3/gcp/param.f90 create mode 100644 test/unit/test_gcp.f90 create mode 100644 test/validation/21-ad3.coord create mode 100644 test/validation/21-gcp-pbeh3c.json create mode 100644 test/validation/21-gcp-pbeh3c.resp create mode 100644 test/validation/22-benI-acetone.xyz create mode 100644 test/validation/22-gcp-b973c.json create mode 100644 test/validation/22-gcp-b973c.resp create mode 100644 test/validation/23-444.sdf create mode 100644 test/validation/23-gcp-hf3c.json create mode 100644 test/validation/23-gcp-hf3c.resp create mode 100644 test/validation/24-gcp-hf-deftzvp.json create mode 100644 test/validation/24-gcp-hf-deftzvp.resp create mode 100644 test/validation/24-h2o_6.qchem create mode 100644 test/validation/25-g3.gen create mode 100644 test/validation/25-gcp-b3lyp-631gd.json create mode 100644 test/validation/25-gcp-b3lyp-631gd.resp diff --git a/README.md b/README.md index 36d5289e..d34053b3 100644 --- a/README.md +++ b/README.md @@ -15,6 +15,8 @@ Usable via the command line interface of the [`s-dftd3` executable](man/s-dftd3. in Fortran via the [`dftd3` module](https://dftd3.readthedocs.io/en/latest/api/fortran.html), with the C interface by including the [`dftd3.h` header](https://dftd3.readthedocs.io/en/latest/api/c.html), or in Python by using the [`dftd3` module](https://dftd3.readthedocs.io/en/latest/api/python.html). +Additionally, the geometric counter-poise correction (see [*JCP* **136**, 154101 (2012)](https://dx.doi.org/10.1063/1.3700154)) +is available to correct for basis set superposition errors and construct composite electronic structure methods of the 3c family. ## Installation @@ -237,6 +239,11 @@ To evaluate a dispersion correction in C four objects are available: Damping parameter object determining the short-range behaviour of the dispersion correction. Standard damping parameters like the rational damping are independent of the molecular structure and can easily be reused for several structures or easily exchanged. +5. the counter-poise parameters: + + Counter-poise parameter object determining the basis set specific correction for basis set superposition error. + Recreating a structure object requires to recreate the counter-poise parameters as well as they are dependent on the basis definition for each element type. + The user is responsible for creating and deleting the objects to avoid memory leaks. diff --git a/app/cli.f90 b/app/cli.f90 index a3733368..cf6530aa 100644 --- a/app/cli.f90 +++ b/app/cli.f90 @@ -21,11 +21,11 @@ module dftd3_app_cli use dftd3, only : d3_param use dftd3_app_argument, only : argument_list, len use dftd3_app_help, only : prog_name, header, help_text, run_help_text, param_help_text, & - & version + & gcp_help_text, version implicit none private - public :: app_config, run_config, param_config, get_arguments + public :: app_config, run_config, param_config, gcp_config, get_arguments type, abstract :: app_config end type app_config @@ -37,6 +37,8 @@ module dftd3_app_cli integer, allocatable :: input_format !> Method name character(len=:), allocatable :: method + !> Basis name + character(len=:), allocatable :: basis !> Damping paramaters type(d3_param) :: inp logical :: json = .false. @@ -56,6 +58,7 @@ module dftd3_app_cli integer :: verbosity = 2 logical :: pair_resolved = .false. logical :: citation = .false. + logical :: gcp = .false. character(len=:), allocatable :: citation_output !> Parameter data base character(len=:), allocatable :: db @@ -70,6 +73,26 @@ module dftd3_app_cli character(len=:), allocatable :: damping end type param_config + type, extends(app_config) :: gcp_config + !> Geometry input file + character(len=:), allocatable :: input + !> Format of the geometry input + integer, allocatable :: input_format + !> Method name + character(len=:), allocatable :: method + !> Basis name + character(len=:), allocatable :: basis + logical :: json = .false. + character(len=:), allocatable :: json_output + logical :: wrap = .true. + logical :: tmer = .true. + logical :: grad = .false. + character(len=:), allocatable :: grad_output + integer :: verbosity = 2 + logical :: citation = .false. + character(len=:), allocatable :: citation_output + end type gcp_config + contains @@ -131,6 +154,9 @@ subroutine get_arguments(config, error) case("param") allocate(param_config :: config) exit + case("gcp") + allocate(gcp_config :: config) + exit end select end do if (allocated(error)) return @@ -146,6 +172,8 @@ subroutine get_arguments(config, error) call get_run_arguments(config, list, iarg, error) type is(param_config) call get_param_arguments(config, list, iarg, error) + type is(gcp_config) + call get_gcp_arguments(config, list, iarg, error) end select end subroutine get_arguments @@ -266,6 +294,17 @@ subroutine get_run_arguments(config, list, start, error) call get_argument_as_real(arg, config%inp%s9, error) if (allocated(error)) exit config%atm = .true. + case("--gcp") + config%gcp = .true. + iarg = iarg + 1 + call list%get(iarg, arg) + if (allocated(arg)) then + if (arg(1:1) == "-") then + iarg = iarg - 1 + cycle + end if + call move_alloc(arg, config%basis) + end if case("--zero") config%zero = .true. iarg = iarg + 1 @@ -489,6 +528,115 @@ subroutine get_param_arguments(config, list, start, error) end subroutine get_param_arguments +subroutine get_gcp_arguments(config, list, start, error) + + !> Configuation data + type(gcp_config), intent(out) :: config + + !> List of command line arguments + type(argument_list), intent(in) :: list + + !> First command line argument + integer, intent(in) :: start + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: iarg, narg + character(len=:), allocatable :: arg + + iarg = start + narg = len(list) + do while(iarg < narg) + iarg = iarg + 1 + call list%get(iarg, arg) + select case(arg) + case("--help") + call info_message(error, gcp_help_text) + exit + case("--version") + call version(output_unit) + stop + case("-v", "--verbose") + config%verbosity = config%verbosity + 1 + case("-s", "--silent") + config%verbosity = config%verbosity - 1 + case default + if (.not.allocated(config%input)) then + call move_alloc(arg, config%input) + cycle + end if + if (arg(1:1) == "-") then + call fatal_error(error, "Unknown argument encountered: '"//arg//"'") + else + call fatal_error(error, "Too many positional arguments present") + end if + exit + case("-i", "--input") + iarg = iarg + 1 + call list%get(iarg, arg) + if (.not.allocated(arg)) then + call fatal_error(error, "Missing argument for input format") + exit + end if + config%input_format = get_filetype("."//arg) + case("-l", "--level") + iarg = iarg + 1 + call list%get(iarg, arg) + if (.not.allocated(arg)) then + call fatal_error(error, "Missing argument for input format") + exit + end if + if (index(arg, "/") > 0) then + config%method = arg(1:index(arg, "/")-1) + config%basis = arg(index(arg, "/")+1:) + else + call move_alloc(arg, config%method) + end if + case("--json") + config%json = .true. + config%json_output = "dftd3.json" + iarg = iarg + 1 + call list%get(iarg, arg) + if (allocated(arg)) then + if (arg(1:1) == "-") then + iarg = iarg - 1 + cycle + end if + call move_alloc(arg, config%json_output) + end if + case("--nocpc") + config%tmer = .false. + case("--nowrap") + config%wrap = .false. + case("--grad") + config%grad = .true. + iarg = iarg + 1 + call list%get(iarg, arg) + if (allocated(arg)) then + if (arg(1:1) == "-") then + iarg = iarg - 1 + cycle + end if + call move_alloc(arg, config%grad_output) + end if + end select + end do + if (allocated(error)) return + + if (config%grad.and. .not.config%json .and. .not.allocated(config%grad_output)) then + config%grad_output = "dftd3.txt" + end if + + if (.not.allocated(config%input)) then + if (.not.allocated(error)) then + write(output_unit, '(a)') gcp_help_text + call fatal_error(error, "Insufficient arguments provided") + end if + end if + +end subroutine get_gcp_arguments + subroutine info_message(error, message) type(error_type), allocatable, intent(out) :: error character(len=*), intent(in) :: message diff --git a/app/driver.f90 b/app/driver.f90 index dc157da2..bf13f579 100644 --- a/app/driver.f90 +++ b/app/driver.f90 @@ -25,15 +25,16 @@ module dftd3_app_driver & get_optimizedpower_damping, optimizedpower_damping_param, & & new_optimizedpower_damping, new_d3_model, get_pairwise_dispersion, & & realspace_cutoff, get_lattice_points + use dftd3_gcp, only : gcp_param, get_gcp_param, get_geometric_counterpoise use dftd3_output, only : ascii_damping_param, ascii_atomic_radii, & & ascii_atomic_references, ascii_system_properties, ascii_energy_atom, & & ascii_results, ascii_pairwise, tagged_result, json_results, & - & turbomole_gradient, turbomole_gradlatt + & turbomole_gradient, turbomole_gradlatt, ascii_gcp_param use dftd3_utils, only : wrap_to_central_cell use dftd3_citation, only : format_bibtex, is_citation_present, citation_type, & & get_citation, doi_dftd3_0, doi_dftd3_bj, doi_dftd3_m, doi_dftd3_op, same_citation use dftd3_app_help, only : header - use dftd3_app_cli, only : app_config, run_config, param_config, get_arguments + use dftd3_app_cli, only : app_config, run_config, param_config, gcp_config, get_arguments use dftd3_app_toml, only : param_database implicit none private @@ -51,6 +52,8 @@ subroutine app_driver(config, error) call run_driver(config, error) type is(param_config) call param_driver(config, error) + type is(gcp_config) + call gcp_driver(config, error) end select end subroutine app_driver @@ -62,6 +65,7 @@ subroutine run_driver(config, error) class(damping_param), allocatable :: param type(d3_param) :: inp type(d3_model) :: d3 + type(gcp_param) :: gcp real(wp), allocatable :: energies(:), gradient(:, :), sigma(:, :) real(wp), allocatable :: pair_disp2(:, :), pair_disp3(:, :) real(wp), allocatable :: s9 @@ -179,9 +183,15 @@ subroutine run_driver(config, error) end block end if end if + if (config%gcp) then + call get_gcp_param(gcp, mol, config%method, config%basis) + end if - if (allocated(param) .and. config%verbosity > 0) then - call ascii_damping_param(output_unit, param, config%method) + if (config%verbosity > 0) then + if (allocated(param)) & + call ascii_damping_param(output_unit, param, config%method) + if (config%gcp) & + call ascii_gcp_param(output_unit, mol, gcp) end if if (allocated(param)) then @@ -200,6 +210,9 @@ subroutine run_driver(config, error) if (allocated(param)) then call get_dispersion(mol, d3, param, realspace_cutoff(), energies, & & gradient, sigma) + if (config%gcp) then + call get_geometric_counterpoise(mol, gcp, realspace_cutoff(), energies, gradient, sigma) + end if energy = sum(energies) if (config%pair_resolved) then @@ -217,50 +230,15 @@ subroutine run_driver(config, error) end if end if if (config%tmer) then - if (config%verbosity > 0) then - write(output_unit, '(a)') "[Info] Dispersion energy written to .EDISP" - end if - open(file=".EDISP", newunit=unit) - write(unit, '(f24.14)') energy - close(unit) + call tmer_writer(".EDISP", energy, "Dispersion", config%verbosity) end if if (config%grad) then if (allocated(config%grad_output)) then - open(file=config%grad_output, newunit=unit) - call tagged_result(unit, energy, gradient, sigma) - close(unit) - if (config%verbosity > 0) then - write(output_unit, '(a)') & - & "[Info] Dispersion results written to '"//config%grad_output//"'" - end if + call results_writer(config%grad_output, energy, gradient, sigma, & + & "Dispersion", config%verbosity) end if - inquire(file="gradient", exist=exist) - if (exist) then - call turbomole_gradient(mol, "gradient", energy, gradient, stat) - if (config%verbosity > 0) then - if (stat == 0) then - write(output_unit, '(a)') & - & "[Info] Dispersion gradient added to Turbomole gradient file" - else - write(output_unit, '(a)') & - & "[Warn] Could not add to Turbomole gradient file" - end if - end if - end if - inquire(file="gradlatt", exist=exist) - if (exist) then - call turbomole_gradlatt(mol, "gradlatt", energy, sigma, stat) - if (config%verbosity > 0) then - if (stat == 0) then - write(output_unit, '(a)') & - & "[Info] Dispersion virial added to Turbomole gradlatt file" - else - write(output_unit, '(a)') & - & "[Warn] Could not add to Turbomole gradlatt file" - end if - end if - end if + call turbomole_writer(mol, energy, gradient, sigma, config%verbosity, "Dispersion") end if if (config%json) then @@ -374,4 +352,185 @@ subroutine from_db(param, input, method, damping, error) end if end subroutine from_db +subroutine gcp_driver(config, error) + type(gcp_config), intent(in) :: config + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(gcp_param) :: param + integer :: unit + real(wp) :: energy + real(wp), allocatable :: energies(:), gradient(:, :), sigma(:, :) + + if (config%verbosity > 1) then + call header(output_unit) + end if + + if (config%input == "-") then + if (.not.allocated(config%input_format)) then + call read_structure(mol, input_unit, filetype%xyz, error) + else + call read_structure(mol, input_unit, config%input_format, error) + end if + else + call read_structure(mol, config%input, error, config%input_format) + end if + if (allocated(error)) return + if (config%wrap) then + call wrap_to_central_cell(mol%xyz, mol%lattice, mol%periodic) + end if + + call get_gcp_param(param, mol, config%method, config%basis) + if (config%verbosity > 0) then + call ascii_gcp_param(output_unit, mol, param) + end if + + allocate(energies(mol%nat), source=0.0_wp) + if (config%grad) then + allocate(gradient(3, mol%nat), source=0.0_wp) + allocate(sigma(3, 3), source=0.0_wp) + end if + + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energies, gradient, sigma) + energy = sum(energies) + + if (config%verbosity > 0) then + if (config%verbosity > 2) then + call ascii_energy_atom(output_unit, mol, energies, label="counter-poise") + end if + call ascii_results(output_unit, mol, energy, gradient, sigma, label="Counter-poise") + end if + + if (config%tmer) then + call tmer_writer(".CPC", energy, "Counter-poise", config%verbosity) + end if + + if (config%grad) then + if (allocated(config%grad_output)) then + call results_writer(config%grad_output, energy, gradient, sigma, & + & "Counter-poise", config%verbosity) + end if + + call turbomole_writer(mol, energy, gradient, sigma, config%verbosity, "Counter-poise") + end if + + if (config%json) then + open(file=config%json_output, newunit=unit) + call json_results(unit, " ", energy=energy, gradient=gradient, sigma=sigma) + close(unit) + if (config%verbosity > 0) then + write(output_unit, '(a)') & + & "[Info] JSON dump of results written to '"//config%json_output//"'" + end if + end if +end subroutine gcp_driver + +!> Write the energy to a file +subroutine tmer_writer(filename, energy, label, verbosity) + !> File name to write to + character(len=*), intent(in) :: filename + + !> Energy to write + real(wp), intent(in) :: energy + + !> Label for the output + character(len=*), intent(in) :: label + + !> Printout verbosity + integer, intent(in) :: verbosity + + integer :: unit + + if (verbosity > 0) then + if (verbosity > 1) then + write(output_unit, '(a)') "[Info] Writing "//label//" energy to '"//filename//"'" + end if + open(file=filename, newunit=unit) + write(unit, '(f24.14)') energy + close(unit) + end if +end subroutine tmer_writer + +!> Write the results to a file +subroutine results_writer(filename, energy, gradient, sigma, label, verbosity) + !> File name to write to + character(len=*), intent(in) :: filename + + !> Energy to write + real(wp), intent(in) :: energy + + !> Gradient to write + real(wp), intent(in) :: gradient(:, :) + + !> Virial tensor to write + real(wp), intent(in) :: sigma(:, :) + + !> Label for the output + character(len=*), intent(in) :: label + + !> Printout verbosity + integer, intent(in) :: verbosity + + integer :: unit + + open(file=filename, newunit=unit) + call tagged_result(unit, energy, gradient, sigma) + close(unit) + if (verbosity > 0) then + write(output_unit, '(a)') "[Info] "//label//" results written to '"//filename//"'" + end if +end subroutine results_writer + +!> Write the gradient and virial tensor to Turbomole files +subroutine turbomole_writer(mol, energy, gradient, sigma, verbosity, label) + + !> Molecular structure data + class(structure_type), intent(in) :: mol + + !> Energy of the system + real(wp), intent(in) :: energy + + !> Gradient of the system + real(wp), intent(in) :: gradient(:, :) + + !> Virial tensor of the system + real(wp), intent(in) :: sigma(:, :) + + !> Printout verbosity + integer, intent(in) :: verbosity + + !> Label for the output + character(len=*), intent(in) :: label + + logical :: exist + integer :: stat + + inquire(file="gradient", exist=exist) + if (exist) then + call turbomole_gradient(mol, "gradient", energy, gradient, stat) + if (verbosity > 0) then + if (stat == 0) then + write(output_unit, '(a)') & + & "[Info] "//label//" gradient added to Turbomole gradient file" + else + write(output_unit, '(a)') & + & "[Warn] Could not add to Turbomole gradient file" + end if + end if + end if + inquire(file="gradlatt", exist=exist) + if (exist) then + call turbomole_gradlatt(mol, "gradlatt", energy, sigma, stat) + if (verbosity > 0) then + if (stat == 0) then + write(output_unit, '(a)') & + & "[Info] "//label//" virial added to Turbomole gradlatt file" + else + write(output_unit, '(a)') & + & "[Warn] Could not add to Turbomole gradlatt file" + end if + end if + end if +end subroutine turbomole_writer + end module dftd3_app_driver diff --git a/app/help.f90 b/app/help.f90 index da619c24..5fb60c79 100644 --- a/app/help.f90 +++ b/app/help.f90 @@ -20,7 +20,7 @@ module dftd3_app_help private public :: prog_name, header, version - public :: help_text, run_help_text, param_help_text + public :: help_text, run_help_text, param_help_text, gcp_help_text character(len=*), parameter :: prog_name = "s-dftd3" @@ -46,6 +46,7 @@ module dftd3_app_help " expected order is s6, s8, a1, a2, bet (requires five arguments)"//nl//& " --atm Use ATM three-body dispersion"//nl//& " --atm-scale Use scaled ATM three-body dispersion"//nl//& + " --gcp Include geometric counter-poise correction for given basis set"//nl//& " --db Load parameters from external data file"//nl//& " --noedisp Disable writing of dispersion energy to .EDISP file"//nl//& " --json [file] Dump results to JSON output (default: dftd3.json)"//nl//& @@ -109,8 +110,29 @@ module dftd3_app_help " d3.op = {s6=1.0, s8=1.44765, a1=0.600, a2=2.50, bet=0.0}"//nl//& "" + character(len=*), parameter :: gcp_help_text = & + "Usage: "//prog_name//" gcp [options] "//nl//& + ""//nl//& + "Takes an geometry input to calculate the geometric counter-poise correction."//nl//& + "Periodic calculations are performed automatically for periodic input formats."//nl//& + "Specify the level of theory to select the correct parameters."//nl//& + ""//nl//& + "-i,--input Hint for the format of the input file"//nl//& + "-l,--level [/]"//nl//& + " Level of theory, basis set is inferred for 3c methods"//nl//& + " --nocpc Disable writing of dispersion energy to .CPC file"//nl//& + " --json [file] Dump results to JSON output (default: gcp.json)"//nl//& + " --grad [file] Request gradient evaluation,"//nl//& + " write results to file (default: gcp.txt),"//nl//& + " attempts to add to Turbomole gradient and gradlatt files"//nl//& + "-v,--verbose Show more, can be used multiple times"//nl//& + "-s,--silent Show less, use twice to supress all output"//nl//& + " --version Print program version and exit"//nl//& + " --help Show this help message"//nl//& + "" + character(len=*), parameter :: help_text = & - "Usage: "//prog_name//" [run|param] [options] ..."//nl//& + "Usage: "//prog_name//" [run|param|gcp] [options] ..."//nl//& ""//nl//& "Commands"//nl//& ""//nl//& diff --git a/doc/_static/references.bib b/doc/_static/references.bib index b903a85d..fcd348c4 100644 --- a/doc/_static/references.bib +++ b/doc/_static/references.bib @@ -214,6 +214,93 @@ @article{hourahine2020 url = {https://doi.org/10.1063/1.5143190}, } +@article{kruse2012, + title={A geometrical correction for the inter-and intra-molecular basis set superposition error in {Hartree-Fock} and density functional theory calculations for large systems}, + author={Kruse, Holger and Grimme, Stefan}, + journal={J. Chem. Phys.}, + volume={136}, + number={15}, + year={2012}, + publisher={AIP Publishing}, + doi = {10.1063/1.3700154}, + url = {https://doi.org/10.1063/1.3700154} +} + +@article{brandenburg2013, + title={Geometrical correction for the inter-and intramolecular basis set superposition error in periodic density functional theory calculations}, + author={Brandenburg, Jan Gerit and Alessio, Maristella and Civalleri, Bartolomeo and Peintinger, Michael F and Bredow, Thomas and Grimme, Stefan}, + journal={J. Phys. Chem. A}, + volume={117}, + number={38}, + pages={9282--9292}, + year={2013}, + publisher={ACS Publications}, + doi={10.1021/jp406658y}, + url={https://doi.org/10.1021/jp406658y} +} + +@article{sure2013, + title={Corrected small basis set {Hartree-Fock} method for large systems}, + author={Sure, Rebecca and Grimme, Stefan}, + journal={J. Computat. Chem.}, + volume={34}, + number={19}, + pages={1672--1685}, + year={2013}, + publisher={Wiley Online Library}, + doi={10.1002/jcc.23317}, + url={https://doi.org/10.1002/jcc.23317} +} + +@article{brandenburg2018, + title={{B97-3c}: A revised low-cost variant of the B97-D density functional method}, + author={Brandenburg, Jan Gerit and Bannwarth, Christoph and Hansen, Andreas and Grimme, Stefan}, + journal={J. Chem. Phys.}, + volume={148}, + number={6}, + year={2018}, + publisher={AIP Publishing}, + doi={10.1063/1.5012601}, + url={https://doi.org/10.1063/1.5012601} +} + +@article{grimme2021, + title={{r2SCAN-3c}: A ``Swiss army knife'' composite electronic-structure method}, + author={Grimme, Stefan and Hansen, Andreas and Ehlert, Sebastian and Mewes, Jan-Michael}, + journal={J. Chem. Phys.}, + volume={154}, + number={6}, + year={2021}, + publisher={AIP Publishing}, + doi={10.1063/5.0040021}, + url={https://doi.org/10.1063/5.0040021} +} + +@article{grimme2015, + title={Consistent structures and interactions by density functional theory with small atomic orbital basis sets}, + author={Grimme, Stefan and Brandenburg, Jan Gerit and Bannwarth, Christoph and Hansen, Andreas}, + journal={J. Chem. Phys.}, + volume={143}, + number={5}, + year={2015}, + publisher={AIP Publishing}, + doi={10.1063/1.4927476}, + url={https://doi.org/10.1063/1.4927476} +} + +@article{brandenburg2016, + title={Screened exchange hybrid density functional for accurate and efficient structures and interaction energies}, + author={Brandenburg, Jan Gerit and Caldeweyher, Eike and Grimme, Stefan}, + journal={Phys. Chem. Chem. Phys.}, + volume={18}, + number={23}, + pages={15519--15523}, + year={2016}, + publisher={Royal Society of Chemistry}, + doi={10.1039/c6cp01697a}, + url={https://doi.org/10.1039/c6cp01697a} +} + @article{golokesh2022, author = {Santra, Golokesh and Semidalas, Emmanouil and Mehta, Nisha and Karton, Amir and Martin, Jan M. L.}, title = {{S66x8} noncovalent interactions revisited: new benchmark and performance of composite localized coupled-cluster methods}, diff --git a/doc/api/c.rst b/doc/api/c.rst index b8e009e1..c9835207 100644 --- a/doc/api/c.rst +++ b/doc/api/c.rst @@ -5,7 +5,7 @@ The C API bindings are provided by using the ``iso_c_binding`` intrinsic module. Generally, objects are exported as opaque pointers and can only be manipulated within the library. The API user is required delete all objects created in the library by using the provided deconstructor functions to avoid mamory leaks. -Overall four classes of objects are provided by the library +Overall five classes of objects are provided by the library - error handlers (:c:type:`dftd3_error`), used to communicate exceptional conditions and errors from the library to the user @@ -16,6 +16,8 @@ Overall four classes of objects are provided by the library general model for calculating dispersion releated properties - damping function objects (:c:type:`dftd3_param`) polymorphic objects to represent the actual method parametrisation +- counter-poise parameter objects (:c:type:`dftd3_gcp`), + short range correction for compensating basis set superposition errors .. note:: @@ -263,6 +265,41 @@ Standard damping parameters like the rational damping are independent of the mol Delete damping parameters. The handle is set to NULL after deletion. +Geometrical counter-poise correction +------------------------------------ + +.. c:type:: struct _dftd3_gcp* dftd3_gcp; + + Counter-poise parameter class + +The counter-poise parameter object provides an additional short ranged correction to account for basis set superposition error in small basis sets. + +.. c:function:: dftd3_gcp dftd3_load_gcp_param(dftd3_error error, dftd3_structure mol, char* method, char* basis); + + :param error: Error handle + :param mol: Molecular structure data handle + :param method: Name of the method to load parameters for + :param basis: Name of the basis to load parameters for + :returns: New counter-poise parameter handle + + Load geometrical counter-poise parameters from internal storage + +.. c:function:: void dftd3_set_gcp_realspace_cutoff(dftd3_error error, dftd3_gcp gcp, double bas, double srb); + + :param error: Error handle + :param model: Dispersion model handle + :param bas: Cutoff for basis set superposition correction + :param srb: Cutoff for short-range bond correction + + Set realspace cutoffs for usage in the counter-poise calculation + +.. c:function:: void dftd3_delete_gcp(dftd3_gcp* gcp); + + :param param: Counter-poise parameter handle + + Delete counter-poise parameters. The handle is set to NULL after deletion. + + Calculation entrypoints ----------------------- @@ -291,6 +328,17 @@ To evaluate dispersion energies or related properties the :c:func:`dftd3_get_dis Evaluate the pairwise representation of the dispersion energy +.. c:function:: void dftd3_get_counterpoise(dftd3_error error, dftd3_structure mol, dftd3_gcp gcp, double* energy, double* gradient, double* sigma); + + :param error: Error handle + :param mol: Molecular structure data handle + :param gcp: Counter-poise parameter handle + :param energy: Dispersion energy + :param gradient: Dispersion gradient [natoms, 3] (optional) + :param sigma: Dispersion strain derivatives [3, 3] (optional) + + Evaluate the counter-poise energy and its derivatives. + Memory management ----------------- diff --git a/doc/index.rst b/doc/index.rst index 09ef0eec..b6b44d74 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -1,8 +1,9 @@ Reimplementation of D3 dispersion correction ============================================ -This pages describe the usage and functionality of `s-dftd3`_ reimplementation of the D3 dispersion correction. +This pages describe the usage and functionality of `s-dftd3`_ reimplementation of the D3 dispersion correction and geometric counter-poise correction. The *s-dftd3* project aims to provide a user-friendly and uniform interface to the D3 dispersion model and for the calculation of DFT-D3 dispersion corrections. +Additionally, it provides the geometric counter-poise correction to create composite electronic structure methods of the 3c family. .. _s-dftd3: https://github.com/dftd3/simple-dftd3 diff --git a/fpm.toml b/fpm.toml index 9cc003ab..a45930b8 100644 --- a/fpm.toml +++ b/fpm.toml @@ -24,3 +24,8 @@ dependencies.toml-f.git = "https://github.com/toml-f/toml-f" [[test]] name = "tester" source-dir = "test/unit" + +[[test]] +name = "api-tester" +source-dir = "test/api" +main = "api-test.c" diff --git a/include/s-dftd3.h b/include/s-dftd3.h index 78a57b26..45d11e38 100644 --- a/include/s-dftd3.h +++ b/include/s-dftd3.h @@ -29,6 +29,7 @@ #define SDFTD3_API_SUFFIX__V_0_3 #define SDFTD3_API_SUFFIX__V_0_4 #define SDFTD3_API_SUFFIX__V_0_5 +#define SDFTD3_API_SUFFIX__V_1_3 /// Error handle class typedef struct _dftd3_error* dftd3_error; @@ -39,6 +40,9 @@ typedef struct _dftd3_structure* dftd3_structure; /// Dispersion model class typedef struct _dftd3_model* dftd3_model; +/// Counter-poisecorrection parameters class +typedef struct _dftd3_gcp* dftd3_gcp; + /// Damping parameter class typedef struct _dftd3_param* dftd3_param; @@ -50,7 +54,8 @@ typedef struct _dftd3_param* dftd3_param; dftd3_error: dftd3_delete_error, \ dftd3_structure: dftd3_delete_structure, \ dftd3_model: dftd3_delete_model, \ - dftd3_param: dftd3_delete_param \ + dftd3_param: dftd3_delete_param, \ + dftd3_gcp: dftd3_delete_gcp \ )(&ptr) /* @@ -218,6 +223,28 @@ dftd3_load_optimizedpower_damping(dftd3_error /* error */, SDFTD3_API_ENTRY void SDFTD3_API_CALL dftd3_delete_param(dftd3_param* /* param */) SDFTD3_API_SUFFIX__V_0_2; +/* + * Counter-poise correction parameters +**/ + +/// Load geometric counter-poise parameters from internal storage +SDFTD3_API_ENTRY dftd3_gcp SDFTD3_API_CALL +dftd3_load_gcp_param(dftd3_error /* error */, + dftd3_structure /* mol */, + char* /* method */, + char* /* basis */) SDFTD3_API_SUFFIX__V_1_3; + +/// Set realspace cutoffs (quantities in Bohr) +SDFTD3_API_ENTRY void SDFTD3_API_CALL +dftd3_set_gcp_realspace_cutoff(dftd3_error /* error */, + dftd3_gcp /* gcp */, + double /* bas */, + double /* srb */) SDFTD3_API_SUFFIX__V_1_3; + +/// Delete counter-poise parameters +SDFTD3_API_ENTRY void SDFTD3_API_CALL +dftd3_delete_gcp(dftd3_gcp* /* gcp */) SDFTD3_API_SUFFIX__V_1_3; + /* * Perform dispersion calculations **/ @@ -240,3 +267,16 @@ dftd3_get_pairwise_dispersion(dftd3_error /* error */, dftd3_param /* param */, double* /* pair_energy2[n][n] */, double* /* pair_energy3[n][n] */) SDFTD3_API_SUFFIX__V_0_5; + +/* + * Perform geometric counterpoise calculations +**/ + +/// Evaluate the dispersion energy and its derivatives +SDFTD3_API_ENTRY void SDFTD3_API_CALL +dftd3_get_counterpoise(dftd3_error /* error */, + dftd3_structure /* mol */, + dftd3_gcp /* gcp */, + double* /* energy */, + double* /* gradient[n][3] */, + double* /* sigma[3][3] */) SDFTD3_API_SUFFIX__V_1_3; \ No newline at end of file diff --git a/src/dftd3/CMakeLists.txt b/src/dftd3/CMakeLists.txt index f15715f0..b47e687c 100644 --- a/src/dftd3/CMakeLists.txt +++ b/src/dftd3/CMakeLists.txt @@ -16,6 +16,7 @@ add_subdirectory("damping") add_subdirectory("data") +add_subdirectory("gcp") set(dir "${CMAKE_CURRENT_SOURCE_DIR}") @@ -26,6 +27,7 @@ list( "${dir}/damping.f90" "${dir}/data.f90" "${dir}/disp.f90" + "${dir}/gcp.f90" "${dir}/model.f90" "${dir}/ncoord.f90" "${dir}/output.f90" diff --git a/src/dftd3/api.f90 b/src/dftd3/api.f90 index 432e2db9..82f907e0 100644 --- a/src/dftd3/api.f90 +++ b/src/dftd3/api.f90 @@ -31,6 +31,7 @@ module dftd3_api use dftd3_damping_zero, only : zero_damping_param, new_zero_damping use dftd3_damping, only : damping_param use dftd3_disp, only : get_dispersion, get_pairwise_dispersion + use dftd3_gcp, only : gcp_param, get_gcp_param, get_geometric_counterpoise use dftd3_model, only : d3_model, new_d3_model use dftd3_param, only : d3_param, get_rational_damping, get_zero_damping, & & get_mrational_damping, get_mzero_damping, get_optimizedpower_damping @@ -58,6 +59,10 @@ module dftd3_api public :: new_optimizedpower_damping_api, load_optimizedpower_damping_api public :: delete_param_api + public :: vp_gcp + public :: load_gcp_param_api, delete_gcp_api, set_gcp_realspace_cutoff + public :: get_counterpoise_api + !> Void pointer to error handle type :: vp_error @@ -85,6 +90,14 @@ module dftd3_api class(damping_param), allocatable :: ptr end type vp_param + !> Void pointer to counter-poise parameters + type :: vp_gcp + !> Actual payload + type(gcp_param) :: ptr + !> Additional real space cutoff + type(realspace_cutoff), allocatable :: cutoff + end type vp_gcp + character(len=*), parameter :: namespace = "dftd3_" @@ -817,6 +830,138 @@ subroutine get_pairwise_dispersion_api(verror, vmol, vdisp, vparam, & end subroutine get_pairwise_dispersion_api +!> Create new error handle object +function load_gcp_param_api(verror, vmol, cmethod, cbasis) & + & result(vgcp) & + & bind(C, name=namespace//"load_gcp_param") + type(c_ptr), value :: verror + type(vp_error), pointer :: error + type(c_ptr), value :: vmol + type(vp_structure), pointer :: mol + character(kind=c_char), intent(in), optional :: cmethod(*) + character(len=:, kind=c_char), allocatable :: method + character(kind=c_char), intent(in), optional :: cbasis(*) + character(len=:, kind=c_char), allocatable :: basis + type(vp_gcp), pointer :: gcp + type(c_ptr) :: vgcp + + vgcp = c_null_ptr + + if (.not.c_associated(verror)) return + call c_f_pointer(verror, error) + + if (.not.c_associated(vmol)) then + call fatal_error(error%ptr, "Molecular structure data is missing") + return + end if + call c_f_pointer(vmol, mol) + + if (present(cmethod)) call c_f_character(cmethod, method) + if (present(cbasis)) call c_f_character(cbasis, basis) + + allocate(gcp) + call get_gcp_param(gcp%ptr, mol%ptr, method, basis) + vgcp = c_loc(gcp) + +end function load_gcp_param_api + + +subroutine set_gcp_realspace_cutoff(verror, vgcp, bas, srb) & + & bind(C, name=namespace//"set_gcp_realspace_cutoff") + type(c_ptr), value :: verror + type(vp_error), pointer :: error + type(c_ptr), value :: vgcp + type(vp_gcp), pointer :: gcp + real(c_double), value, intent(in) :: bas + real(c_double), value, intent(in) :: srb + + if (.not.c_associated(verror)) return + call c_f_pointer(verror, error) + + if (.not.c_associated(vgcp)) then + call fatal_error(error%ptr, "Counter-poise parameters are missing") + return + end if + call c_f_pointer(vgcp, gcp) + + gcp%cutoff = realspace_cutoff(gcp=bas, srb=srb) +end subroutine set_gcp_realspace_cutoff + + +!> Delete counter-poise parameter handle object +subroutine delete_gcp_api(vgcp) & + & bind(C, name=namespace//"delete_gcp") + type(c_ptr), intent(inout) :: vgcp + type(vp_error), pointer :: gcp + + if (c_associated(vgcp)) then + call c_f_pointer(vgcp, gcp) + + deallocate(gcp) + vgcp = c_null_ptr + end if + +end subroutine delete_gcp_api + + +!> Calculate dispersion +subroutine get_counterpoise_api(verror, vmol, vgcp, & + & energy, c_gradient, c_sigma) & + & bind(C, name=namespace//"get_counterpoise") + type(c_ptr), value :: verror + type(vp_error), pointer :: error + type(c_ptr), value :: vmol + type(vp_structure), pointer :: mol + type(c_ptr), value :: vgcp + type(vp_gcp), pointer :: gcp + real(c_double), intent(out) :: energy + real(c_double), intent(out), optional :: c_gradient(3, *) + real(wp), allocatable :: gradient(:, :) + real(c_double), intent(out), optional :: c_sigma(3, 3) + real(wp), allocatable :: sigma(:, :) + type(realspace_cutoff) :: cutoff + + if (.not.c_associated(verror)) return + call c_f_pointer(verror, error) + + if (.not.c_associated(vmol)) then + call fatal_error(error%ptr, "Molecular structure data is missing") + return + end if + call c_f_pointer(vmol, mol) + + if (.not.c_associated(vgcp)) then + call fatal_error(error%ptr, "Counter-poise parameters are missing") + return + end if + call c_f_pointer(vgcp, gcp) + + if (present(c_gradient)) then + gradient = c_gradient(:3, :mol%ptr%nat) + endif + + if (present(c_gradient)) then + sigma = c_sigma(:3, :3) + endif + + cutoff = realspace_cutoff() + if (allocated(gcp%cutoff)) then + cutoff = gcp%cutoff + end if + call get_geometric_counterpoise(mol%ptr, gcp%ptr, cutoff, & + & energy, gradient, sigma) + + if (present(c_gradient)) then + c_gradient(:3, :mol%ptr%nat) = gradient + endif + + if (present(c_gradient)) then + c_sigma(:3, :3) = sigma + endif + +end subroutine get_counterpoise_api + + subroutine f_c_character(rhs, lhs, len) character(kind=c_char), intent(out) :: lhs(*) character(len=*), intent(in) :: rhs diff --git a/src/dftd3/cutoff.f90 b/src/dftd3/cutoff.f90 index b5501b0b..8d99a18c 100644 --- a/src/dftd3/cutoff.f90 +++ b/src/dftd3/cutoff.f90 @@ -30,6 +30,12 @@ module dftd3_cutoff !> Three-body interaction cutoff real(wp), parameter :: disp3_default = 40.0_wp + !> Counter-poise correction cutoff + real(wp), parameter :: gcp_default = 60.0_wp + + !> Short-range bond correction cutoff + real(wp), parameter :: srb_default = 60.0_wp + !> Collection of real space cutoffs type :: realspace_cutoff @@ -44,6 +50,12 @@ module dftd3_cutoff !> Three-body interaction cutoff real(wp) :: disp3 = disp3_default + !> Counter-poise correction cutoff + real(wp) :: gcp = gcp_default + + !> Short-range bond correction cutoff + real(wp) :: srb = srb_default + end type realspace_cutoff diff --git a/src/dftd3/gcp.f90 b/src/dftd3/gcp.f90 new file mode 100644 index 00000000..293ba66a --- /dev/null +++ b/src/dftd3/gcp.f90 @@ -0,0 +1,1431 @@ +! This file is part of s-dftd3. +! SPDX-Identifier: LGPL-3.0-or-later +! +! s-dftd3 is free software: you can redistribute it and/or modify it under +! the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! s-dftd3 is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with s-dftd3. If not, see . + +module dftd3_gcp + use mctc_env, only : wp + use mctc_io, only : structure_type + use dftd3_gcp_param, only : gcp_param, get_gcp_param + use dftd3_cutoff, only : realspace_cutoff, get_lattice_points + implicit none + private + + public :: gcp_param, get_gcp_param, get_geometric_counterpoise + + !> Geometric counterpoise correction + interface get_geometric_counterpoise + module procedure get_geometric_counterpoise + module procedure get_geometric_counterpoise_atomic + end interface get_geometric_counterpoise + + + +contains + + +!> Geometric counterpoise correction +subroutine get_geometric_counterpoise(mol, param, cutoff, energy, gradient, sigma) + + !> Molecular structure data + class(structure_type), intent(in) :: mol + + !> Geometric counterpoise parameters + type(gcp_param), intent(in) :: param + + !> Realspace cutoffs + type(realspace_cutoff), intent(in) :: cutoff + + !> Counter-poise energy + real(wp), intent(inout) :: energy + + !> Counter-poise gradient + real(wp), intent(inout), optional :: gradient(:, :) + + !> Counter-poise virial + real(wp), intent(inout), optional :: sigma(:, :) + + real(wp), allocatable :: energies(:) + + allocate(energies(mol%nat), source=0.0_wp) + call get_geometric_counterpoise_atomic(mol, param, cutoff, energies, gradient, sigma) + energy = energy + sum(energies) +end subroutine get_geometric_counterpoise + + +!> Geometric counterpoise correction with atom-resolved energies +subroutine get_geometric_counterpoise_atomic(mol, param, cutoff, energies, gradient, sigma) + + !> Molecular structure data + class(structure_type), intent(in) :: mol + + !> Geometric counterpoise parameters + type(gcp_param), intent(in) :: param + + !> Realspace cutoffs + type(realspace_cutoff), intent(in) :: cutoff + + !> Dispersion energy + real(wp), intent(inout) :: energies(:) + + !> Dispersion gradient + real(wp), intent(inout), optional :: gradient(:, :) + + !> Dispersion virial + real(wp), intent(inout), optional :: sigma(:, :) + + real(wp), allocatable :: lattr(:, :) + real(wp), parameter :: rexp_base = 0.75_wp, zexp_base = 1.5_wp, rexp_srb = -1.0_wp, zexp_srb = 0.5_wp + logical :: grad + + grad = present(gradient) .and. present(sigma) + + if (allocated(param%emiss) .and. allocated(param%slater) .and. allocated(param%xv)) then + call get_lattice_points(mol%periodic, mol%lattice, cutoff%gcp, lattr) + if (grad) then + call gcp_deriv(mol, lattr, cutoff%gcp, param%zeff, param%emiss, param%slater, & + & param%xv, param%rvdw, param%sigma, param%alpha, param%beta, param%damp, & + & param%dmp_scal, param%dmp_exp, energies, gradient, sigma) + else + call gcp_energy(mol, lattr, cutoff%gcp, param%zeff, param%emiss, param%slater, & + & param%xv, param%rvdw, param%sigma, param%alpha, param%beta, param%damp, & + & param%dmp_scal, param%dmp_exp, energies) + end if + end if + + if (param%srb .or. param%base) then + call get_lattice_points(mol%periodic, mol%lattice, cutoff%srb, lattr) + end if + + if (param%srb) then + if (grad) then + call srb_deriv(mol, lattr, cutoff%srb, mol%num, param%rvdw_srb, param%rscal, param%qscal, & + & rexp_srb, zexp_srb, energies, gradient, sigma) + else + call srb_energy(mol, lattr, cutoff%srb, mol%num, param%rvdw_srb, param%rscal, param%qscal, & + & rexp_srb, zexp_srb, energies) + end if + end if + + if (param%base) then + if (grad) then + call srb_deriv(mol, lattr, cutoff%srb, param%zeff, param%rvdw, param%rscal, param%qscal, & + & rexp_base, zexp_base, energies, gradient, sigma) + else + call srb_energy(mol, lattr, cutoff%srb, param%zeff, param%rvdw, param%rscal, param%qscal, & + & rexp_base, zexp_base, energies) + end if + end if +end subroutine get_geometric_counterpoise_atomic + + +!> Geometric counterpoise correction +subroutine gcp_energy(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, alpha, beta, & + & damp, dmp_scal, dmp_exp, energies) + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Translation vectors + real(wp), intent(in) :: trans(:, :) + !> Distance cutoff + real(wp), intent(in) :: cutoff + !> Effective nuclear charges + integer, intent(in) :: iz(:) + !> Basis set superposition error per atom + real(wp), intent(in) :: emiss(:) + !> Slater exponents + real(wp), intent(in) :: slater(:) + !> Number of virtual orbitals + real(wp), intent(in) :: xv(:) + !> Van der Waals radii + real(wp), intent(in) :: rvdw(:, :) + !> Scaling factor + real(wp), intent(in) :: escal + !> Exponential factor + real(wp), intent(in) :: alpha + !> Power factor + real(wp), intent(in) :: beta + !> Damping flag + logical, intent(in) :: damp + !> Damping scaling factor + real(wp), intent(in) :: dmp_scal + !> Damping exponent + real(wp), intent(in) :: dmp_exp + !> Atom-resolved energy + real(wp), intent(inout) :: energies(:) + + integer :: iat, jat, jtr, izp, jzp + real(wp) :: argv + real(wp) :: xvi, xvj + real(wp) :: r1, rscal, rscalexp, r0 + real(wp) :: sij, expv, bsse, emi, emj + real(wp) :: vec(3) + real(wp) :: dampval, grd_dmp + real(wp) :: dE + + do iat = 1, mol%nat + izp = mol%id(iat) + xvi = merge(1.0_wp / sqrt(xv(izp)), 0.0_wp, xv(izp) >= 0.5_wp) + + ! the BSSE due to atom jat, Loop over all j atoms + do jat = 1, iat + jzp = mol%id(jat) + xvj = merge(1.0_wp / sqrt(xv(jzp)), 0.0_wp, xv(jzp) >= 0.5_wp) + + emi = emiss(izp)*xvj*escal + emj = emiss(jzp)*xvi*escal + r0 = rvdw(izp, jzp) + do jtr = 1, size(trans, 2) + vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) + r1 = norm2(vec) + + if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle + ! calulate slater overlap sij + call ssovl(r1, izp, jzp, iz, slater(izp), slater(jzp), sij) + ! evaluate gcp central expression + argv = -alpha*r1**beta + expv = exp(argv) + bsse = expv/sqrt(sij) + if (damp) then + rscal = r1/r0 + rscalexp = dmp_scal*rscal**dmp_exp + dampval = (1.0_wp-1.0_wp/(1.0_wp+rscalexp)) + else + dampval = 1.0_wp + endif + dE = bsse*dampval + energies(iat) = energies(iat) + emi*dE + if (iat /= jat) then + energies(jat) = energies(jat) + emj*dE + end if + end do + end do + end do + +end subroutine gcp_energy + + +!> Geometric counterpoise correction +subroutine gcp_deriv(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, alpha, beta, & + & damp, dmp_scal, dmp_exp, energies, gradient, sigma) + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Translation vectors + real(wp), intent(in) :: trans(:, :) + !> Distance cutoff + real(wp), intent(in) :: cutoff + !> Effective nuclear charges + integer, intent(in) :: iz(:) + !> Basis set superposition error per atom + real(wp), intent(in) :: emiss(:) + !> Slater exponents + real(wp), intent(in) :: slater(:) + !> Number of virtual orbitals + real(wp), intent(in) :: xv(:) + !> Van der Waals radii + real(wp), intent(in) :: rvdw(:, :) + !> Scaling factor + real(wp), intent(in) :: escal + !> Exponential factor + real(wp), intent(in) :: alpha + !> Power factor + real(wp), intent(in) :: beta + !> Damping flag + logical, intent(in) :: damp + !> Damping scaling factor + real(wp), intent(in) :: dmp_scal + !> Damping exponent + real(wp), intent(in) :: dmp_exp + !> Atom-resolved energy + real(wp), intent(inout) :: energies(:) + !> Molecular gradient + real(wp), intent(inout) :: gradient(:, :) + !> Virial + real(wp), intent(inout) :: sigma(:, :) + + integer :: iat, jat, jtr, izp, jzp + real(wp) :: argv, argd, ovlpd + real(wp) :: xvi, xvj + real(wp) :: r1, rscal, rscalexp, rscalexpm1, r0 + real(wp) :: sij, expv, expd, bsse, emi, emj, emij + real(wp) :: vec(3), gs(3), gij + real(wp) :: dampval, grd_dmp + real(wp) :: dE, dG(3), dS(3, 3) + + do iat = 1, mol%nat + izp = mol%id(iat) + xvi = merge(1.0_wp / sqrt(xv(izp)), 0.0_wp, xv(izp) >= 0.5_wp) + + ! the BSSE due to atom jat, Loop over all j atoms + do jat = 1, iat + jzp = mol%id(jat) + xvj = merge(1.0_wp / sqrt(xv(jzp)), 0.0_wp, xv(jzp) >= 0.5_wp) + + emi = emiss(izp)*xvj*escal + emj = emiss(jzp)*xvi*escal + emij = emi + emj + r0 = rvdw(izp, jzp) + do jtr = 1, size(trans, 2) + vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) + r1 = norm2(vec) + + if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle + ! calulate slater overlap sij + call ssovl(r1, izp, jzp, iz, slater(izp), slater(jzp), sij) + ! evaluate gcp central expression + argv = -alpha*r1**beta + expv = exp(argv) + bsse = expv/sqrt(sij) + if (damp) then + rscal = r1/r0 + rscalexp = dmp_scal*rscal**dmp_exp + dampval = (1.0_wp-1.0_wp/(1.0_wp+rscalexp)) + else + dampval = 1.0_wp + endif + + call gsovl(r1, izp, jzp, iz, slater(izp), slater(jzp), gij) + + gs(:) = gij*vec + expd = exp(-alpha*r1**beta)*(-0.5_wp) + argd = 2d0*alpha*beta*r1**beta*sij/r1 + ovlpd = r1*sij**1.5_wp + + if(damp) then + rscalexpm1 = rscal**(dmp_exp-1) + grd_dmp = dmp_scal*dmp_exp*rscalexpm1/r0 + grd_dmp = grd_dmp/(rscalexp+1.0_wp)**2 + endif + + dE = bsse*dampval + dG = expd*(argd*vec + gs)/ovlpd*emij + if(damp) then + dG = dG*dampval+bsse*grd_dmp*(vec/r1)*emij + endif + dS = spread(dG, 1, 3) * spread(vec, 2, 3) * 0.5_wp + + energies(iat) = energies(iat) + emi*dE + if (iat /= jat) then + energies(jat) = energies(jat) + emj*dE + end if + sigma(:, :) = sigma(:, :) + dS + if (iat /= jat) then + gradient(:, iat) = gradient(:, iat) + dG + gradient(:, jat) = gradient(:, jat) - dG + sigma(:, :) = sigma(:, :) + dS + end if + end do + end do + end do + +end subroutine gcp_deriv + + +!> Short-range bond length correction for HF-3c +subroutine srb_energy(mol, trans, cutoff, iz, r0ab, rscal, qscal, rexp, zexp, energies) + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Translation vectors + real(wp), intent(in) :: trans(:, :) + !> Distance cutoff + real(wp), intent(in) :: cutoff + !> Effective nuclear charges + integer, intent(in) :: iz(:) + !> Van der Waals radii + real(wp), intent(in) :: r0ab(:, :) + !> Radii scaling factor + real(wp), intent(in) :: rscal + !> Prefactor for the SRB potential + real(wp), intent(in) :: qscal + !> Exponent for radii + real(wp), intent(in) :: rexp + !> Exponent for charges + real(wp), intent(in) :: zexp + !> Atom-resolved energy + real(wp), intent(inout) :: energies(:) + + real(wp) :: fi, fj, ff, r1, expt + real(wp) :: r0, vec(3), dE + integer iat, jat, jtr, izp, jzp + + do iat = 1, mol%nat + izp = mol%id(iat) + fi = real(iz(izp), wp) + do jat = 1, iat + jzp = mol%id(jat) + r0 = rscal*r0ab(izp, jzp)**rexp + fj = real(iz(jzp), wp) + ff = -(fi*fj)**zexp + + do jtr = 1, size(trans, 2) + vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) + r1 = norm2(vec) + if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle + expt = exp(-r0*r1) + dE = qscal*ff*expt*0.5_wp + + energies(iat) = energies(iat) + dE + if (iat /= jat) then + energies(jat) = energies(jat) + dE + end if + end do + end do + end do + +end subroutine srb_energy + + +!> Short-range bond length correction for HF-3c +subroutine srb_deriv(mol, trans, cutoff, iz, r0ab, rscal, qscal, rexp, zexp, energies, gradient, sigma) + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Translation vectors + real(wp), intent(in) :: trans(:, :) + !> Distance cutoff + real(wp), intent(in) :: cutoff + !> Effective nuclear charges + integer, intent(in) :: iz(:) + !> Van der Waals radii + real(wp), intent(in) :: r0ab(:, :) + !> Radii scaling factor + real(wp), intent(in) :: rscal + !> Prefactor for the SRB potential + real(wp), intent(in) :: qscal + !> Exponent for radii + real(wp), intent(in) :: rexp + !> Exponent for charges + real(wp), intent(in) :: zexp + !> Atom-resolved energy + real(wp), intent(inout) :: energies(:) + !> Molecular gradient + real(wp), intent(inout) :: gradient(:, :) + !> Molecular virial + real(wp), intent(inout) :: sigma(:, :) + + real(wp) :: fi, fj, ff, rf, r1, expt + real(wp) :: r0, vec(3), dE, dG(3), dS(3, 3) + integer iat, jat, jtr, izp, jzp + + do iat = 1, mol%nat + izp = mol%id(iat) + fi = real(iz(izp), wp) + do jat = 1, iat + jzp = mol%id(jat) + r0 = rscal*r0ab(izp, jzp)**rexp + fj = real(iz(jzp), wp) + ff = -(fi*fj)**zexp + + do jtr = 1, size(trans, 2) + vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) + r1 = norm2(vec) + if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle + expt = exp(-r0*r1) + rf = qscal/r1 + + dE = qscal*ff*expt*0.5_wp + dG(:) = -ff*r0*vec*expt*rf + dS(:, :) = spread(dG, 1, 3) * spread(vec, 2, 3) * 0.5_wp + + energies(iat) = energies(iat) + dE + if (iat /= jat) then + energies(jat) = energies(jat) + dE + end if + sigma(:, :) = sigma + dS + if (iat /= jat) then + gradient(:, iat) = gradient(:, iat) + dG + gradient(:, jat) = gradient(:, jat) - dG + sigma(:, :) = sigma + dS + end if + + end do + end do + end do + +end subroutine srb_deriv + + +!****************************************************************************** +!* calculates the s-type overlap integral over 1s, 2s and 3s slater functions +!* added support for 3s functions +!* ovl = overlap integral +!* za = slater exponent atom A +!* zb = slater exponent atom B +!* R = distance between atom A and B +!* Inspired by mopac7.0 +!****************************************************************************** +subroutine ssovl(r, iat, jat, iz, xza, xzb, ovl) + integer ii, shell(72) + logical debug + real(wp) za, zb, R, ovl, ax, bx, norm, R05 + integer na, nb + real(wp) Bxx0, Bxx1, Bxx2, xx, Bxx4, Bxx6 + real(wp) Bxx3, Bxx5 + data shell/ & + ! h, he + 1, 1 & + ! li-ne + , 2, 2, 2, 2, 2, 2, 2, 2, & + ! na-ar + 3, 3, 3, 3, 3, 3, 3, 3, & + ! 4s, 5s will be treated as 3s + ! k-rn , no f-elements + 54*3/ + ! ... + real(wp) xza, xzb + integer iat, jat, iz(*) + + za = xza + zb = xzb + na = iz(iat) + nb = iz(jat) + debug = .false. +!debug = .true. + +! ii selects kind of ovl by multiplying the shell +! kind <1s|1s> <2s|1s> <2s|2s> <1s|3s> <2s|3s> <3s|3s> +! case: 1 2 4 3 6 9 +! + ii = shell(na)*shell(nb) + if(debug) write(*, *) 'shell', ii + + R05 = R*0.5 + ax = (za+zb)*R05 + bx = (zb-za)*R05 + + ! same elements + if(abs(za-zb) < 0.1) then + select case (ii) + case (1) + ovl = 0.25d0*sqrt((za*zb*R*R)**3)*(A2(ax)*Bint(bx, 0)-Bint(bx, 2)*A0(ax)) + case (2) + ovl = SQRT(1._wp/3._wp) + if(shell(na) < shell(nb)) then + ! <1s|2s> + norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp + ovl = ovl*norm*(A3(ax)*Bint(bx, 0)-Bint(bx, 3)*A0(ax)+A2(ax)*Bint(bx, 1)-Bint(bx, 2)*A1(ax)) + else + ! switch za/zb to get <2s|1s> + xx = za + za = zb + zb = xx + ax = (za+zb)*R05 + bx = (zb-za)*R05 + norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp + ovl = ovl*norm*(A3(ax)*Bint(bx, 0)-Bint(bx, 3)*A0(ax)+A2(ax)*Bint(bx, 1)-Bint(bx, 2)*A1(ax)) + endif + case (4) + norm = SQRT((ZA*ZB)**5)*(R**5)*0.0625d0 + ovl = norm* (A4(ax)*Bint(bx, 0)+Bint(bx, 4)*A0(ax)-2.0d0*A2(ax)*Bint(bx, 2))*(1d0/3d0) + case(3) + if(shell(na) < shell(nb)) then + norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp + ovl = norm*(A4(ax)*Bint(bx, 0)-Bint(bx, 4)*A0(ax)+2.d0*(A3(ax)*Bint(bx, 1)-Bint(bx, 3)*A1(ax)))/sqrt(3.d0) + else + xx = za + za = zb + zb = xx + ax = (za+zb)*R05 + bx = (zb-za)*R05 + norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp + ovl = norm*(A4(ax)*Bint(bx, 0)-Bint(bx, 4)*A0(ax)+2.d0*(A3(ax)*Bint(bx, 1)-Bint(bx, 3)*A1(ax)))/sqrt(3.d0) + endif + case(6) + if(shell(na) < shell(nb)) then + norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp + ovl = norm*(A5(ax)*Bint(bx, 0)+A4(ax)*Bint(bx, 1) & + & -2d0*(A3(ax)*Bint(bx, 2)+A2(ax)*Bint(bx, 3)) & + & +A1(ax)*Bint(bx, 4)+A0(ax)*Bint(bx, 5))/3.d0 + else + xx = za + za = zb + zb = xx + ax = (za+zb)*R05 + bx = (zb-za)*R05 + norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp + ovl = norm*(A5(ax)*Bint(bx, 0)+A4(ax)*Bint(bx, 1) & + & -2d0*(A3(ax)*Bint(bx, 2)+A2(ax)*Bint(bx, 3)) & + & +A1(ax)*Bint(bx, 4)+A0(ax)*Bint(bx, 5))/3.d0 + endif + case(9) + norm = sqrt((ZA*ZB*R*R)**7)/480.d0 + ovl = norm*(A6(ax)*Bint(bx, 0)-3.d0*(A4(ax)*Bint(bx, 2) & + & -A2(ax)*Bint(bx, 4))-A0(ax)*Bint(bx, 6))/3._wp + end select + else ! different elements + select case (ii) + case (1) + norm = 0.25d0*sqrt((za*zb*R*R)**3) + ovl = (A2(ax)*B0(bx)-B2(bx)*A0(ax))*norm + case (2) + ovl = SQRT(1._wp/3._wp) + if(shell(na) < shell(nb)) then + ! <1s|2s> + norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp + ovl = ovl*norm*(A3(ax)*B0(bx)-B3(bx)*A0(ax)+A2(ax)*B1(bx)-B2(bx)*A1(ax)) + else + ! switch za/zb to get <2s|1s> + xx = za + za = zb + zb = xx + ax = (za+zb)*R05 + bx = (zb-za)*R05 + norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp + ovl = ovl*norm*(A3(ax)*B0(bx)-B3(bx)*A0(ax)+A2(ax)*B1(bx)-B2(bx)*A1(ax)) + endif + case (4) ! <2s|2s> + norm = SQRT((ZA*ZB)**5)*(R**5)*0.0625_wp + ovl = norm* (A4(ax)*B0(bx)+B4(bx)*A0(ax)-2.0_wp*A2(ax)*B2(bx))*(1d0/3d0) + case(3) ! <1s|3s> + <3s|1s> + if(shell(na) < shell(nb)) then + norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp + ovl = norm*(A4(ax)*B0(bx)-B4(bx)*A0(ax)+2.d0*(A3(ax)*B1(bx)-B3(bx)*A1(ax)))/sqrt(3.d0) + else + xx = za + za = zb + zb = xx + ax = (za+zb)*R05 + bx = (zb-za)*R05 + norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp + ovl = norm*(A4(ax)*B0(bx)-B4(bx)*A0(ax)+2.d0*(A3(ax)*B1(bx)-B3(bx)*A1(ax)))/sqrt(3.d0) + endif + case(6) ! <2s|3s> + <3s|2s> + if(shell(na) < shell(nb)) then + norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp + ovl = norm*(A5(ax)*B0(bx)+A4(ax)*B1(bx)-2d0*(A3(ax)*B2(bx)+A2(ax)*B3(bx))+A1(ax)*B4(bx)+A0(ax)*B5(bx))/3.d0 + else + xx = za + za = zb + zb = xx + ax = (za+zb)*R05 + bx = (zb-za)*R05 + norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp + ovl = norm*(A5(ax)*B0(bx)+A4(ax)*B1(bx)-2.0_wp*(A3(ax)*B2(bx)+A2(ax)*B3(bx))+A1(ax)*B4(bx)+A0(ax)*B5(bx))/3.d0 + endif + case(9) ! <3s|3> + norm = sqrt((ZA*ZB*R*R)**7)/1440.d0 + ovl = norm*(A6(ax)*B0(bx)-3.d0*(A4(ax)*B2(bx)-A2(ax)*B4(bx))-A0(ax)*B6(bx)) + end select + endif +end subroutine ssovl + + +!**************************************** +!* A(x) auxiliary integrals * +!* Quantenchemie - Ein Lehrgang Vol 5 * +!* p. 570 eq. 11.4.14 * +!**************************************** + +real(wp) pure function A0(x) +! Hilfsintegral A_0 +implicit none +real(wp), intent(in) :: x +A0 = exp(-x)/x +return +end function + +real(wp) pure function A1(x) +! Hilfsintegral A_1 +implicit none +real(wp), intent(in) :: x +A1 = ((1+x)*exp(-x))/(x**2) +return +end function + + +real(wp) pure function A2(x) +! Hilfsintegral A_2 +implicit none +real(wp), intent(in) :: x +A2 = ((2d0+2d0*x+x**2)*exp(-x))/x**3 +return +end function + + +real(wp) pure function A3(x) +! Hilfsintegral A_3 +implicit none +real(wp), intent(in) :: x +real(wp) xx +real(wp) x2, x3, x4 +x2 = x*x +x3 = x2*x +x4 = x3*x +xx = (6d0+6d0*x+3d0*x2+x3) +A3 = (xx*exp(-x))/x4 +return +end function + + +real(wp) pure function A4(x) +! Hilfsintegral A_4 +implicit none +real(wp), intent(in) :: x +real(wp) xx +real(wp) x2, x3, x4, x5 +x2 = x*x +x3 = x2*x +x4 = x3*x +x5 = x4*x +xx = (24d0+24d0*x+12d0*x2+4d0*x3+x4) +A4 = (xx*exp(-x))/x5 +return +end function + +real(wp) pure function A5(x) +! Hilfsintegral A_5 +implicit none +real(wp), intent(in) :: x +real(wp) xx +real(wp) x2, x3, x4, x5, x6 +x2 = x*x +x3 = x2*x +x4 = x3*x +x5 = x4*x +x6 = x5*x +xx = (120d0+120d0*x+60d0*x2+20d0*x3+5d0*x4+x5) +A5 = (xx*exp(-x))/x6 +return +end function + +real(wp) pure function A6(x) +! Hilfsintegral A_6 +implicit none +real(wp), intent(in) :: x +real(wp) xx +real(wp) x2, x3, x4, x5, x6, x7 +x2 = x*x +x3 = x2*x +x4 = x3*x +x5 = x4*x +x6 = x5*x +x7 = x6*x +xx = (720d0+720d0*x+360d0*x2+120d0*x3+30d0*x4+6d0*x5+x6) +A6 = (xx*exp(-x))/x7 +return +end function + + + +!************************************** +!* B(x) auxiliary integrals * +!* Quantenchemie - Ein Lehrgang Vol 5 * +!* p. 570 eq. 11.4.14b * +!************************************** + + +real(wp) pure function B0(x) + real(wp), intent(in) :: x + B0 = (exp(x)-exp(-x))/x +end function + +real(wp) pure function B1(x) + real(wp), intent(in) :: x + real(wp) x2, x3 + x2 = x*x + x3 = x2*x + B1 = ((1.0_wp-x)*exp(x)-(1.0_wp+x)*exp(-x))/x2 +end function + +real(wp) pure function B2(x) + real(wp), intent(in) :: x + real(wp) x2, x3 + x2 = x*x + x3 = x2*x + B2 = (((2.0_wp-2*x+x2)*exp(x)) - ((2.0_wp+2.0_wp*x+x2)*exp(-x)))/x3 +end function + +real(wp) pure function B3(x) + real(wp), intent(in) :: x + real(wp) xx, yy + real(wp) x2, x3, x4 + x2 = x*x + x3 = x2*x + x4 = x3*x + xx = (6.0_wp-6.0_wp*x+3.0_wp*x2-x3)*exp(x)/x4 + yy = (6.0_wp+6.0_wp*x+3.0_wp*x2+x3)*exp(-x)/x4 + B3 = xx-yy +end function + + +real(wp) pure function B4(x) + real(wp), intent(in) :: x + real(wp) xx, yy + real(wp) x2, x3, x4, x5 + x2 = x*x + x3 = x2*x + x4 = x3*x + x5 = x4*x + xx = (24.0_wp-24.0_wp*x+12.0_wp*x2-4.0_wp*x3+x4)*exp(x)/x5 + yy = (24.0_wp+24.0_wp*x+12.0_wp*x2+4.0_wp*x3+x4)*exp(-x)/x5 + B4 = xx-yy +end function + +real(wp) pure function B5(x) + real(wp), intent(in) :: x + real(wp) xx, yy + real(wp) x2, x3, x4, x5, x6 + x2 = x*x + x3 = x2*x + x4 = x3*x + x5 = x4*x + x6 = x5*x + xx = (120.0_wp-120*x+60*x2-20*x3+5*x4-x5)*exp(x)/x6 + yy = (120.0_wp+120*x+60*x2+20*x3+5*x4+x5)*exp(-x)/x6 + B5 = xx-yy +end function + +real(wp) function B6(x) + real(wp), intent(in) :: x + real(wp) x2, x3, x4, x5, x6, x7, yy, xx + x2 = x*x + x3 = x2*x + x4 = x3*x + x5 = x4*x + x6 = x5*x + x7 = x6*x + xx = (720.0_wp - 720.0_wp*x+ 360.0_wp*x2 - 120.0_wp*x3 + 30.0_wp*x4 - 6.0_wp*x5 + x6)*exp(x)/x7 + yy = (720.0_wp + 720.0_wp*x + 360.0_wp*x2 + 120.0_wp*x3 + 30.0_wp*x4 + 6.0_wp*x5 + x6)*exp(-x)/x7 + B6 = xx-yy +end function + + +real(wp) function bint(x, k) +! calculates B_k(x) +! general summation formula +! 'infinite' sum is numerically unstable. 12 terms seem +! accurate enough +implicit none +real(wp), intent(in) :: x +real(wp) xx, yy +integer, intent(in) :: k +integer i +bint = 0 + +if(abs(x).lt.1e-6) then +do i = 0, k + bint = (1.d0+(-1d0)**i)/(dble(i)+1.d0) +end do +return +endif + +do i = 0, 12 +xx = 1d0-((-1d0)**(k+i+1)) +yy = dble(fact(i))*dble((k+i+1)) +bint = bint+xx/yy*(-x)**i +enddo + + +end function bint + + +! faculty function +integer(wp) function fact(N) +implicit none +integer j, n +fact = 1 +do j = 2, n + fact = fact*j +enddo +return + +end + + + + +subroutine gsovl(r, iat, jat, iz, xza, xzb, g) + ! GRADIENT + ! calculates the s-type overlap integral over 1s, 2s, 3s slater functions + ! ovl = overlap integral + ! za = slater exponent atom A + ! zb = slater exponent atom B + ! R = distance between atom A and B + integer ii, shell(72) + logical debug + real(wp) ax, bx, R05, za, zb, R + integer na, nb + data shell/ & + ! h, he + 1, 1 & + ! li-ne + , 2, 2, 2, 2, 2, 2, 2, 2, & + ! na-ar + 3, 3, 3, 3, 3, 3, 3, 3, & + ! 4s, 5s will be treated as 3s + ! k-rn , no f-elements + 54*3/ + ! ... + real(wp) g, Fa, Fb + !--------------------- set exponents --------------------------------------- + real(wp) xza, xzb + real(wp) xx + integer iat, jat, iz(*) + logical lsame + + za = xza + zb = xzb + na = iz(iat) + nb = iz(jat) + !---------------------------------------------------------------------------- + + debug = .false. + !debug = .true. + + ! ii selects kind of ovl by multiplying the shell + ! kind <1s|1s> <2s|1s> <2s|2s> <1s|3s> <2s|3s> <3s|3s> + ! case: 1 2 4 3 6 9 + ! + ii = shell(na)*shell(nb) + if(debug) write(*, *) 'gshell', ii + R05 = R*0.5_wp + ax = (za+zb)*R05 + Fa = (za+zb) + bx = (zb-za)*R05 + Fb = (zb-za) + lsame = .false. + ! + ! same elements + if(abs(za-zb) < 0.1) then + lsame = .true. + ! call arguments: gtype(exponents, argumentDeriv., distance, gradient, (Switch shell), sameElement) + select case (ii) + case (1) + call g1s1s(za, zb, Fa, Fb, R, g, lsame) + case (2) + if(shell(na) < shell(nb)) then + call g2s1s(za, zb, Fa, Fb, R, g, .false., lsame) + else + xx = za + za = zb + zb = xx + call g2s1s(za, zb, Fa, Fb, R, g, .true., lsame) + end if + case (4) + call g2s2s(za, zb, Fa, Fb, R, g, lsame) + case(3) + if(shell(na) < shell(nb)) then + call g1s3s(za, zb, Fa, Fb, R, g, .false., lsame) + else + xx = za + za = zb + zb = xx + call g1s3s(za, zb, Fa, Fb, R, g, .true., lsame) + end if + case(6) + if(shell(na) < shell(nb)) then + call g2s3s(za, zb, Fa, Fb, R, g, .false., lsame) + else + xx = za + za = zb + zb = xx + call g2s3s(za, zb, Fa, Fb, R, g, .true., lsame) + end if + case(9) + call g3s3s(za, zb, Fa, Fb, R, g, lsame) + end select + else ! different elements + lsame = .false. + select case (ii) + case (1) + call g1s1s(za, zb, Fa, Fb, R, g, lsame) + case (2) ! <1s|2s> + if(shell(na) < shell(nb)) then + call g2s1s(za, zb, Fa, Fb, R, g, .false., lsame) + else + xx = za + za = zb + zb = xx + call g2s1s(za, zb, Fa, Fb, R, g, .true., lsame) + end if + case (4) ! <2s|2s> + call g2s2s(za, zb, Fa, Fb, R, g, lsame) + case(3) ! <1s|3s> + <3s|1s> + if(shell(na) < shell(nb)) then + call g1s3s(za, zb, Fa, Fb, R, g, .false., lsame) + else + xx = za + za = zb + zb = xx + call g1s3s(za, zb, Fa, Fb, R, g, .true., lsame) + end if + case(6) ! <2s|3s> + <3s|2s> + if(shell(na) < shell(nb)) then + call g2s3s(za, zb, Fa, Fb, R, g, .false., lsame) + else + xx = za + za = zb + zb = xx + call g2s3s(za, zb, Fa, Fb, R, g, .true., lsame) + end if + case(9) ! <3s|3> + call g3s3s(za, zb, Fa, Fb, R, g, lsame) + end select + end if + +end subroutine gsovl + + +!------------------------------------------------------------- +! Maple was used to find the analy. derivatives of +! the slater integrals (expressions over A, B aux. integrals) +! Optimized fortran code by maple with some human-corrections +!------------------------------------------------------------- +subroutine g1s1s(za, zb, Fa, Fb, R, g, sameElement) + ! slater overlap derv. + ! derivative of explicit integral expression + ! using maple + implicit real(wp) (t) + real(wp) za, zb, Fa, Fb + real(wp) g, R + logical sameElement + + if(sameElement) then + t1 = za ** 2 + t3 = zb ** 2 + t5 = t1 * za * t3 * zb + t6 = R ** 2 + t7 = t6 ** 2 + t10 = Fa * R + t14 = exp(-0.5_wp * t10) + t17 = sqrt(t5 * t7 * t6) + g = -(1.0_wp/3.0_wp) * t5 * t7 / Fa * (0.2D1 + t10) * t14 / t17 + else + t1 = za ** 2 + t3 = zb ** 2 + t5 = t1 * za * t3 * zb + t6 = Fb ** 2 + t7 = Fb * R + t8 = 0.5_wp * t7 + t9 = exp(t8) + t12 = exp(-t8) + t15 = t6 * Fa + t22 = Fa ** 2 + t23 = t22 * t9 + t27 = t22 * t12 + t31 = t6 * Fb + t32 = R * t31 + t37 = t22 * Fa + t38 = R * t37 + t43 = R ** 2 + t44 = t43 * t31 + t51 = t43 * t37 + t56 = 0.4D1 * t6 * t9 - 0.4D1 * t6 * t12 + 0.2D1 * t15 * R * t9 - & + 0.2D1 * t15 * R * t12 - 0.4D1 * t23 + 0.2D1 * t23 * t7 + 0.4D1 * t27 & + + 0.2D1 * t27 * t7 - 0.2D1 * t32 * t9 - 0.2D1 * t32 * t12 - 0.2D1 * t38 * & + t9 + 0.2D1 * t38 * t12 - 0.1D1 * t44 * Fa * t9 - 0.1D1 & + * t44 * Fa * t12 + t51 * t9 * Fb + t51 * t12 * Fb + t61 = exp(-0.5_wp * Fa * R) + t62 = t43 ** 2 + t65 = sqrt(t5 * t62 * t43) + g = -0.2D1 * t5 * R * t56 * t61 / t65 / t31 / t37 + end if + + +end subroutine g1s1s + + +subroutine g2s1s(za, zb, Fa, Fb, R, g, switch, lsame) + ! slater overlap derv. + ! derivative of explicit integral expression + ! using maple + implicit real(wp) (t) + real(wp) za, zb, Fa, Fb + real(wp) g, R, norm + logical switch + logical lsame + norm = (1d0/24d0)*sqrt(za**3*zb**5*3d0) + + if(switch) then + Fb = -Fb + endif + + if(lsame) then + + t1 = Fa * R + t3 = exp(-0.5000000000_wp * t1) + t6 = Fa ** 2 + t7 = R ** 2 + g = -0.1000000000D-8 * R * t3 * (0.5333333380D10 + 0.2666666670D10 & + * t1 + 0.1333333333D10 * t6 * t7) / t6 + g = g*norm + + else + + t3 = exp(-0.5000000000_wp * Fa * R) + t4 = Fa ** 2 + t5 = t4 * Fa + t6 = Fb * R + t7 = 0.5000000000_wp * t6 + t8 = exp(t7) + t9 = t5 * t8 + t11 = Fb ** 2 + t12 = t11 * Fa + t15 = exp(-t7) + t18 = t4 ** 2 + t19 = R * t18 + t22 = t11 ** 2 + t29 = Fb * t4 + t36 = R ** 2 + t37 = t36 * t18 + t44 = t36 * R + t48 = -0.12D2 * t9 + 0.4D1 * t12 * t8 - 0.4D1 * t12 * t15 & + - 0.6D1 * t19 * t8 - 0.6D1 * t22 * t8 * R - 0.6D1 * t22 * t15 & + * R + 0.4D1 * t29 * t15 - 0.4D1 * t29 * t8 + 0.6D1 * t19 * t15 & + + 0.2D1 * t37 * t8 * Fb + 0.4D1 * t37 * t15 * Fb + t44 * t18 * t15 * t11 + t49 = t5 * t15 + t51 = t11 * Fb + t58 = t51 * Fa + t59 = R * t8 + t76 = t36 * t15 + t79 = t22 * Fa + t87 = 0.12D2 * t49 - 0.12D2 * t51 * t15 - 0.1D1 * t22 * t4 * t15 * t44 & + + 0.4D1 * t58 * t59 - 0.8D1 * t58 * R * t15 + 0.4D1 * t9 * t6 + 0.8D1 * & + t49 * t6 + 0.2D1 * t49 * t11 * t36 + 0.4D1 * t11 * t4 * t59 - 0.2D1 * t51 & + * t4 * t76 - 0.2D1 * t79 * t36 * t8 - 0.4D1 * t79 * t76 + 0.12D2 * t51 * t8 + g = -0.16D2 * t3 * (t48 + t87) / t36 / t22 / t18 + g = g*norm + end if + +end subroutine g2s1s + +subroutine g2s2s(za, zb, Fa, Fb, R, g, SameElement) + ! slater overlap derv. + ! derivative of explicit integral expression + ! using maple + implicit real(wp) (t) + real(wp) za, zb, Fa, Fb + real(wp) g, R, norm + logical SameElement + + norm = 1d0/(16d0*3d0)*SQRT((ZA*ZB)**5) + + if(SameElement) then + + t2 = R ** 2 + t5 = Fa ** 2 + t9 = t5 * Fa + t10 = t2 ** 2 + t16 = exp(-Fa * R / 0.2D1) + g = (-0.4266666666D2 * R - 0.2133333333D2 * Fa * t2 - 0.2133333333D1 & + * t5 * t2 * R - 0.1066666666D1 * t9 * t10) * t16 / t9 + g = g*norm + + + else + t1 = R ** 2 + t3 = 0.3840000000D3 * t1 * Fb + t4 = t1 * R + t5 = Fb ** 2 + t7 = 0.6400000000D2 * t4 * t5 + t8 = 0.7680000000D3 * R + t10 = Fa ** 2 + t11 = t10 ** 2 + t12 = t11 * Fa + t14 = Fb * R + t15 = 0.768000000D3 * t14 + t17 = 0.1280000000D3 * t5 * t1 + t21 = 0.256000000D3 * t5 * R + t22 = t5 * Fb + t24 = 0.1280000000D3 * t22 * t1 + t26 = t10 * Fa + t28 = t5 ** 2 + t30 = 0.1280000000D3 * t1 * t28 + t32 = 0.256000000D3 * t22 * R + t33 = 0.512000000D3 * t5 + t34 = t28 * Fb + t36 = 0.6400000000D2 * t4 * t34 + t40 = 0.768000000D3 * t28 * R + t42 = 0.3840000000D3 * t1 * t34 + t45 = 0.1536000000D4 * t28 + t47 = 0.7680000000D3 * t34 * R + t51 = exp(-0.5_wp * Fa * R) + t53 = 0.5_wp * t14 + t54 = exp(-t53) + t68 = exp(t53) + g = (((t3 + t7 + t8) * t12 + (0.1536000000D4 + t15 + t17) * t11 + & + (-t21 - t24) * t26 + (t30 - t32 - t33 + t36) * t10 + (t40 + t42) *& + Fa + t45 + t47) * t51 * t54 + ((t3 - t8 - t7) * t12 + (-0.1536000000D4 & + + t15 - t17) * t11 + (-t24 + t21) * t26 + (-t30 + t33 - t32 & + + t36) * t10 + (-t40 + t42) * Fa + t47 - t45) * t51 * t68) / t1 / & + t12 / t34 + + + g = g*norm + endif + +end subroutine g2s2s + + +subroutine g1s3s(za, zb, Fa, Fb, R, g, switch, lsame) + ! slater overlap derv. + ! derivative of explicit integral expression + ! using maple + implicit real(wp) (t) + real(wp) za, zb, Fa, Fb + real(wp) g, R, norm + logical switch + logical lsame + + if(switch) Fb = -Fb + + norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)/(16d0*sqrt(3d0)) + + if(lsame) then + + t1 = Fa * R + t3 = exp(-0.5000000000_wp * t1) + t4 = R ** 2 + g = -0.1600000000D1 * t3 * t4 * R * (0.2D1 + t1) / Fa + g = g*norm + else + + t3 = exp(-0.5000_wp * Fa * R) + t4 = Fb ** 2 + t5 = t4 ** 2 + t6 = t5 * Fb + t7 = t6 * Fa + t8 = R ** 2 + t9 = Fb * R + t10 = 0.50_wp * t9 + t11 = exp(t10) + t15 = exp(-t10) + t16 = t8 * t15 + t19 = Fa ** 2 + t21 = t8 * R + t22 = t21 * t15 + t25 = t19 * Fa + t27 = t8 ** 2 + t31 = t19 ** 2 + t32 = t31 * Fa + t33 = t8 * t32 + t45 = t4 * Fb + t48 = t31 * t15 + t55 = t4 * t25 + t56 = t11 * R + t59 = t15 * R + t62 = t5 * Fa + t73 = -0.6D1 * t7 * t8 * t11 - 0.18D2 * t7 * t16 - 0.6D1 * t6 * t19 & + * t22 - 0.1D1 * t6 * t25 * t27 * t15 + 0.6D1 * t33 * t11 * Fb + 0.18D2 & + * t33 * t15 * Fb + 0.6D1 * t21 * t32 * t15 * t4 + t27 * t32* t15 * t45 & + + 0.2D1 * t48 * t45 * t21 + 0.12D2 * t48 * t4 * t8 + 0.12D2 * t55 * t56 & + + 0.12D2 * t55 * t59 + 0.12D2 * t62 * t56 - 0.36D2 * t62 * t59 - 0.12D2 & + * t5 * t19 * t16 - 0.2D1 * t5 * t25 * t22 + t74 = t31 * t11 + t79 = t45 * t19 + t92 = R * t32 + t95 = t45 * Fa + t100 = Fb * t25 + t111 = 0.12D2 * t74 * t9 + 0.36D2 * t48 * t9 + 0.12D2 * t79 * t56 - 0.12D2 & + * t79 * t59 + 0.48D2 * t5 * t11 - 0.24D2 * t6 * t11 * R - 0.24D2 * t6 * t15 & + * R - 0.24D2 * t92 * t11 + 0.24D2 * t95 * t11 - 0.24D2 * t95 * t15 + 0.24D2 & + * t100 * t15 + 0.24D2 * t92 * t15 - 0.24D2 * t100 * t11 - 0.48D2 * t5 * t15 & + - 0.48D2 * t74 + 0.48D2 * t48 + g = -0.32D2 * t3 * (t73 + t111) / t8 / t6 / t32 + + g = g*norm + endif + +end subroutine g1s3s + + +subroutine g2s3s(za, zb, Fa, Fb, R, g, switch, lsame) + ! slater overlap derv. + ! derivative of explicit integral expression + ! using maple + implicit real(wp) (t) + real(wp) za, zb, Fa, Fb + real(wp) g, R, norm + logical switch + logical lsame + norm = sqrt((za**5)*(zb**7)/7.5_wp)/96.d0 + + if(switch) Fb = -Fb + + if(lsame) then + t1 = Fa * R + t3 = exp(-0.5000000000_wp * t1) + t6 = Fa ** 2 + t7 = R ** 2 + t14 = t6 ** 2 + t15 = t7 ** 2 + g = -0.2000000000D-8 * R * t3 * (0.1280000000D12 + 0.6400000000D11 & + * t1 + 0.1280000000D11 * t6 * t7 + 0.1066666670D10 * t6 * Fa * t7 & + * R + 0.533333333D9 * t14 * t15) / t14 + g = g*norm + else + + t3 = exp(-0.5_wp * Fa * R) + t4 = Fb ** 2 + t5 = t4 ** 2 + t6 = Fa ** 2 + t7 = t6 * Fa + t8 = t5 * t7 + t9 = R ** 2 + t11 = 0.50_wp * Fb * R + t12 = exp(t11) + t13 = t9 * t12 + t16 = t6 ** 2 + t17 = t16 * Fa + t18 = exp(-t11) + t21 = t5 * Fb + t28 = t9 * t18 + t32 = t9 * R + t33 = t32 * t18 + t36 = t5 * t4 + t38 = t9 ** 2 + t39 = t38 * t18 + t41 = t21 * Fa + t42 = R * t12 + t45 = t16 * t6 + t46 = t4 * Fb + t49 = t46 * t16 + t52 = -0.6D1 * t8 * t13 + 0.120D3 * t17 * t18 + 0.120D3 * t21 * t18 & + - 0.120D3 * t17 * t12 - 0.120D3 * t21 * t12 - 0.6D1 * t8 * t28 - 0.2D1 & + * t5 * t16 * t33 + t36 * t7 * t39 - 0.48D2 * t41 * t42 + t45 * t46 * t39 - 0.6D1 * t49 * t13 + t54 = R * t18 + t60 = t46 * t6 + t63 = Fb * t16 + t66 = t5 * Fa + t69 = t4 * t7 + t72 = t36 * t6 + t75 = t32 * t12 + t78 = Fb * t9 + t84 = Fb * t17 + t87 = -0.24D2 * t46 * t7 * t54 - 0.24D2 * t5 * t6 * t42 + 0.24D2 *& + t60 * t12 + 0.24D2 * t63 * t18 - 0.24D2 * t66 * t12 - 0.24D2 * t69 & + * t18 + 0.9D1 * t72 * t33 + 0.3D1 * t72 * t75 + 0.24D2 * t78 * t45 & + * t12 - 0.6D1 * t49 * t28 + 0.48D2 * t84 * t42 + t102 = t21 * t6 + t105 = t4 * t17 + t113 = t45 * t4 + t118 = 0.72D2 * t84 * t54 + 0.72D2 * t41 * t54 + 0.36D2 * t78 * t45 & + * t18 + 0.2D1 * t46 * t17 * t33 + 0.24D2 * t4 * t16 * t42 - 0.6D1 & + * t102 * t13 - 0.6D1 * t105 * t13 + 0.18D2 * t105 * t28 + 0.2D1 & + * t21 * t7 * t33 - 0.3D1 * t113 * t75 + 0.9D1 * t113 * t33 + t121 = t36 * Fa + t130 = R * t45 + t145 = 0.18D2 * t102 * t28 + 0.24D2 * t121 * t13 + 0.36D2 * t121 * & + t28 - 0.24D2 * t60 * t18 - 0.24D2 * t63 * t12 + 0.60D2 * t130 * t18 & + + 0.60D2 * t36 * t18 * R + 0.24D2 * t69 * t12 + 0.60D2 * t36 * t12 * & + R - 0.60D2 * t130 * t12 + 0.24D2 * t66 * t18 + g = 0.128D3 * t3 * (t52 + t87 + t118 + t145) / t9 / t36 / t45 + + g = g*norm + endif + +end subroutine g2s3s + + +subroutine g3s3s(za, zb, Fa, Fb, R, g, SameElement) + ! slater overlap derv. + ! derivative of explicit integral expression + ! using maple + implicit real(wp) (t) + real(wp) za, zb, Fa, Fb + real(wp) g, R, norm + logical SameElement + + norm = sqrt((ZA*ZB)**7)/1440.d0 + + if(SameElement) then + + t1 = Fa * R + t3 = exp(-0.5000000000_wp * t1) + t5 = Fa ** 2 + t6 = t5 ** 2 + t7 = t6 * Fa + t8 = R ** 2 + t9 = t8 ** 2 + g = -0.2000000000D-8 * t3 * R * (0.457142857D9 * t7 * t9 & + * R + 0.7680000000D12 * t1 + 0.1536000000D12 * t5 * t8 & + + 0.1280000000D11 * t5 * Fa * t8 * R + 0.914285715D9 * t6 * t9 + 0.1536000000D13) / t7 + + g = g*norm + else + + t3 = exp(-0.5000000000_wp * Fa * R) + t4 = Fa ** 2 + t5 = t4 ** 2 + t6 = t5 * t4 + t7 = Fb * R + t8 = 0.5000000000_wp * t7 + t9 = exp(-t8) + t10 = t6 * t9 + t13 = Fb ** 2 + t14 = t13 * Fb + t15 = t13 ** 2 + t16 = t15 * t14 + t17 = R ** 2 + t18 = t17 * R + t19 = t16 * t18 + t23 = exp(t8) + t24 = t6 * t23 + t27 = t5 * Fa + t28 = t27 * t13 + t29 = R * t23 + t32 = t6 * t13 + t33 = t17 * t9 + t36 = t15 * Fb + t37 = t4 * t36 + t38 = t9 * R + t43 = t17 * t23 + t46 = t4 * Fa + t47 = t5 * t46 + t48 = t47 * t18 + t52 = t47 * t17 + t65 = 0.120D3 * t10 * t7 - 0.12D2 * t19 * t4 * t9 + 0.120D3 & + * t24 * t7 + 0.24D2 * t28 * t29 + 0.24D2 * t32 * t33 + 0.24D2 * t37 & + * t38 - 0.24D2 * t28 * t38 - 0.24D2 * t32 * t43 - 0.12D2 * t48 * t13 & + * t23 + 0.60D2 * t52 * t23 * Fb + 0.12D2 * t48 * t13 * t9 + 0.60D2 & + * t52 * t9 * Fb - 0.12D2 * t19 * t4 * t23 + t66 = t17 ** 2 + t67 = t16 * t66 + t74 = t27 * t14 + t77 = t6 * t14 + t78 = t18 * t23 + t81 = t46 * t15 + t86 = t27 * t15 + t89 = t5 * t36 + t90 = t18 * t9 + t97 = t46 * t36 + t104 = -0.1D1 * t67 * t46 * t9 - 0.1D1 * t67 * t46 * t23 - 0.12D2 & + * t74 * t43 + 0.2D1 * t77 * t78 - 0.24D2 * t81 * t29 + 0.24D2 * t81 & + * t38 + 0.2D1 * t86 * t78 + 0.2D1 * t89 * t90 - 0.2D1 * t86 * t90 & + + 0.24D2 * t37 * t29 + 0.12D2 * t97 * t33 + 0.2D1 * t89 * t78 - 0.12D2 * t74 * t33 + t108 = t5 * t14 + t111 = t15 * t13 + t112 = t111 * t4 + t117 = t111 * t46 + t122 = t111 * Fa + t129 = t4 * t15 + t132 = t47 * R + t139 = 0.2D1 * t77 * t90 - 0.24D2 * t108 * t38 + 0.24D2 * t112 * t43 & + - 0.24D2 * t112 * t33 + 0.2D1 * t117 * t78 - 0.2D1 * t117 * t90 + 0.120D3 & + * t122 * t29 - 0.120D3 * t122 * t38 + 0.12D2 * t97 * t43 - 0.48D2 * t129 & + * t23 + 0.120D3 * t132 * t9 - 0.120D3 * t132 * t23 + 0.240D3 * t111 * t23 + t140 = t47 * t66 + t145 = t16 * R + t150 = t16 * t17 + t160 = t5 * t13 + t170 = t140 * t14 * t23 + t140 * t14 * t9 - 0.120D3 * t145 * t9 - 0.24D2 & + * t108 * t29 - 0.60D2 * t150 * Fa * t23 - 0.240D3 * t111 * t9 - 0.240D3 & + * t24 + 0.240D3 * t10 + 0.48D2 * t129 * t9 - 0.48D2 * t160 * t9 + 0.48D2 & + * t160 * t23 - 0.120D3 * t145 * t23 - 0.60D2 * t150 * Fa * t9 + g = -0.768D3 * t3 * (t65 + t104 + t139 + t170) / t17 / t47 / t16 + + g = g*norm + endif + +end subroutine g3s3s + +end module dftd3_gcp \ No newline at end of file diff --git a/src/dftd3/gcp/CMakeLists.txt b/src/dftd3/gcp/CMakeLists.txt new file mode 100644 index 00000000..dbe22936 --- /dev/null +++ b/src/dftd3/gcp/CMakeLists.txt @@ -0,0 +1,24 @@ +# This file is part of s-dftd3. +# SPDX-Identifier: LGPL-3.0-or-later +# +# s-dftd3 is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# s-dftd3 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with s-dftd3. If not, see . + +set(dir "${CMAKE_CURRENT_SOURCE_DIR}") + +list( + APPEND srcs + "${dir}/param.f90" +) + +set(srcs "${srcs}" PARENT_SCOPE) diff --git a/src/dftd3/gcp/meson.build b/src/dftd3/gcp/meson.build new file mode 100644 index 00000000..28e21f78 --- /dev/null +++ b/src/dftd3/gcp/meson.build @@ -0,0 +1,19 @@ +# This file is part of s-dftd3. +# SPDX-Identifier: LGPL-3.0-or-later +# +# s-dftd3 is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# s-dftd3 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with s-dftd3. If not, see . + +srcs += files( + 'param.f90', +) diff --git a/src/dftd3/gcp/param.f90 b/src/dftd3/gcp/param.f90 new file mode 100644 index 00000000..5463afe1 --- /dev/null +++ b/src/dftd3/gcp/param.f90 @@ -0,0 +1,1084 @@ +! This file is part of s-dftd3. +! SPDX-Identifier: LGPL-3.0-or-later +! +! s-dftd3 is free software: you can redistribute it and/or modify it under +! the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! s-dftd3 is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with s-dftd3. If not, see . + +module dftd3_gcp_param + use mctc_env, only : wp + use mctc_io, only : structure_type + use dftd3_data, only : get_vdw_rad + implicit none + private + + public :: gcp_param, get_gcp_param + + !> Parameters for the geometric counter-poise correction + type :: gcp_param + !> Basis set superposition error correction + real(wp), allocatable :: emiss(:) + !> Number of virtual orbitals + real(wp), allocatable :: xv(:) + !> Slater exponents + real(wp), allocatable :: slater(:) + !> Van der Waals radii for effective nuclear charges + real(wp), allocatable :: rvdw(:, :) + !> Van der Waals radii for true nuclear charges + real(wp), allocatable :: rvdw_srb(:, :) + !> Effective nuclear charges + integer, allocatable :: zeff(:) + !> Scaling factor for the counter-poise correction + real(wp) :: sigma = 0.0_wp + !> Exponential parameter + real(wp) :: alpha = 0.0_wp + !> Power parameter + real(wp) :: beta = 0.0_wp + !> Damping enabled + logical :: damp = .false. + !> Damping scaling factor + real(wp) :: dmp_scal = 4.0_wp + !> Damping exponent + real(wp) :: dmp_exp = 6.0_wp + !> Short-range bond correction + logical :: srb = .false. + !> Short-range bond correction radii scaling factor + real(wp) :: rscal = 0.0_wp + !> Short-range bond correction scaling factor + real(wp) :: qscal = 0.0_wp + !> Short-range bond correction for HF-3c + logical :: base = .false. + end type + + enum, bind(c) + enumerator :: & + p_unknown_bas, & + p_sv_bas, & + p_sv_p_bas, & + p_svx_bas, & + p_svp_bas, & + p_svp_old_bas, & + p_minis_bas, & + p_631gd_bas, & + p_tz_bas, & + p_deftzvp_bas, & + p_ccdz_bas, & + p_accdz_bas, & + p_pobtz_bas, & + ! p_pobdzvp_bas, & + p_minix_bas, & + p_gcore_bas, & + p_2g_bas, & + p_dzp_bas, & + ! p_hsv_bas, & + p_dz_bas, & + p_msvp_bas, & + p_lanl2_bas, & + p_pbeh3c_bas, & + p_def2mtzvpp_bas, & + p_def2mtzvp_bas, & + p_vmb_bas + end enum + + enum, bind(c) + enumerator :: & + p_unknown_method, & + p_hf_method, & + p_dft_method, & + p_hyb_method, & + p_gga_method, & + p_b3lyp_method, & + p_blyp_method, & + p_pbe_method, & + p_tpss_method, & + p_pw6b95_method, & + p_hf3c_method, & + p_pbeh3c_method, & + p_hse3c_method, & + p_b973c_method, & + p_b3pbe3c_method, & + p_r2scan3c_method + end enum + + real(wp), parameter :: emiss_hf_sv(*) = [ & + & 0.009037_wp, 0.008843_wp, & ! He,He + & 0.204189_wp, 0.107747_wp, 0.049530_wp, 0.055482_wp, 0.072823_wp, 0.100847_wp, 0.134029_wp, 0.174222_wp, & ! Li-Ne + & 0.315616_wp, 0.261123_wp, 0.168568_wp, 0.152287_wp, 0.146909_wp, 0.168248_wp, 0.187882_wp, 0.211160_wp, & !Na -Ar + & 0.374252_wp, 0.460972_wp, & ! K-Ca + & 0.444886_wp, 0.404993_wp, 0.378406_wp, 0.373439_wp, 0.361245_wp, & + & 0.360014_wp, 0.362928_wp, 0.243801_wp, 0.405299_wp, 0.396510_wp, & ! 3d-TM + & 0.362671_wp, 0.360457_wp, 0.363355_wp, 0.384170_wp, 0.399698_wp, 0.417307_wp] !Ga-Kr + + real(wp), parameter :: emiss_hf_minis(*) = [ & + & 0.042400_wp, 0.028324_wp, & + & 0.252661_wp, 0.197201_wp, 0.224237_wp, 0.279950_wp, 0.357906_wp, 0.479012_wp, 0.638518_wp, 0.832349_wp, & + & 1.232920_wp, 1.343390_wp, 1.448280_wp, 1.613360_wp, 1.768140_wp, 1.992010_wp, 2.233110_wp, 2.493230_wp, & + & 3.029640_wp, 3.389980_wp, & ! H-Ca + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp] + + real(wp), parameter :: emiss_hf_631gd(*) = [ &! H-Ca + Br (no 3d) + & 0.010083_wp, 0.008147_wp, & + & 0.069260_wp, 0.030540_wp, 0.032736_wp, 0.021407_wp, 0.024248_wp, 0.036746_wp, 0.052733_wp, 0.075120_wp, & + & 0.104255_wp, 0.070691_wp, 0.100260_wp, 0.072534_wp, 0.054099_wp, 0.056408_wp, 0.056025_wp, 0.057578_wp, & + & 0.079198_wp, 0.161462_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.000000_wp, 0.000000_wp, 0.000000_wp, 0.000000_wp, 0.381049_wp, 0.000000_wp] + + real(wp), parameter :: emiss_hf_svp_old(*) = [real(wp):: & ! Li,Na,Mg,K had wrong parameters + & 0.008107,0.008045,& + & 0.142957,0.028371,0.049369,0.055376,0.072785,0.100310,0.133273,0.173600,& + & 0.191109,0.222839,0.167188,0.149843,0.145396,0.164308,0.182990,0.205668,& + & 0.221189,0.299661,& + & 0.325995,0.305488,0.291723,0.293801,0.29179,0.296729,0.304603,0.242041,0.354186,0.350715,& + & 0.350021,0.345779,0.349532,0.367305,0.382008,0.399709] + + real(wp), parameter :: emiss_hf_svp(*) = [ & ! H-Kr + & 0.008107,0.008045,& + & 0.113583,0.028371,0.049369,0.055376,0.072785,0.100310,0.133273,0.173600,& + & 0.181140,0.125558,0.167188,0.149843,0.145396,0.164308,0.182990,0.205668,& + & 0.200956,0.299661, & + & 0.325995,0.305488,0.291723,0.293801,0.29179,0.296729,0.304603,0.242041,0.354186,0.350715,& + & 0.350021,0.345779,0.349532,0.367305,0.382008,0.399709] + + real(wp), parameter :: emiss_hf_sv_p(*) = [ & ! H-Kr + & emiss_hf_sv(1), emiss_hf_svp(2:)] + + real(wp), parameter :: emiss_hf_svx(*) = [ & ! H-Kr + & emiss_hf_sv(1), emiss_hf_svp(2:5), emiss_hf_sv(6), emiss_hf_svp(7:)] + + real(wp), parameter :: emiss_hf_tz(*) = [ & ! H-Kr !def2-TZVP + & 0.007577,0.003312,& + & 0.086763,0.009962,0.013964,0.005997,0.004731,0.005687,0.006367,0.007511,& + & 0.077721,0.050003,0.068317,0.041830,0.025796,0.025512,0.023345,0.022734,& + & 0.097241,0.099167,& + & 0.219194,0.189098,0.164378,0.147238,0.137298,0.12751,0.118589,0.0318653,0.120985,0.0568313, & + & 0.090996,0.071820,0.063562,0.064241,0.061848,0.061021] + + real(wp), parameter :: emiss_hf_def2mtzvp(*) = [& ! m def2-TZVP, no f for B-Ne + & 0.007930,0.003310,& + & 0.086760,0.009960,0.013960,0.006000,0.003760,0.004430,0.005380,0.006750,& + & 0.077720,0.050000,0.068320,0.041830,0.025800,0.025510,0.023340,0.022730,& + & 0.097240,0.099170,& + & 0.219190,0.189100,0.164380,0.147240,0.137300,0.127510,0.118590,0.031870,0.120990,0.056830,& + & 0.091000,0.071820,0.063560,0.064240,0.061850,0.061020] + + real(wp), parameter :: emiss_hf_vmb(*) = [ & + & 0.042400_wp, 0.028324_wp, & + & 0.252661_wp, 0.197201_wp, 0.156009_wp, 0.164586_wp, 0.169273_wp, 0.214704_wp, 0.729138_wp, 0.336072_wp, & + & 0.262329_wp, 0.231722_wp, 0.251169_wp, 0.287795_wp, 0.213590_wp, 0.250524_wp, 0.728579_wp, 0.260658_wp, & + & 0.0_wp, 0.0_wp,& + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp] + + real(wp), parameter :: emiss_hf_minisd(*) = [& !Al-Ar MINIS + Ahlrichs "P" funktions + & 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 1.446950_wp, 1.610980_wp, 1.766610_wp, 1.988230_wp, 2.228450_wp, 2.487960_wp, & + & 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp] + + real(wp), parameter :: emiss_hf_minix(*) = [ & + & emiss_hf_minis(1:2), 0.177871_wp, 0.171596_wp, emiss_hf_minis(5:10), 1.114110_wp, 1.271150_wp, & + & emiss_hf_minisd(13:18), emiss_hf_sv(19:30), emiss_hf_svp(31:36)] + + real(wp), parameter :: emiss_hf_lanl2(*) = [ & ! LANL2TZ+ vs LANL2DZ (ORCA), only Sc-Zn + & emiss_hf_631gd(1:20), & + & 0.102545_wp, 0.0719529_wp, 0.0491798_wp, 0.0362976_wp, 0.0266369_wp, & + & 0.0235484_wp, 0.0171578_wp, 0.0438906_wp, 0.0100259_wp, 0.016208_wp, & + & emiss_hf_631gd(31:)] + + real(wp), parameter :: emiss_hf_pobtz(*) = [ & ! H-Kr, no RG + & 0.010077_wp, 0.000000_wp, & + & 0.173239_wp, 0.101973_wp, 0.131181_wp, 0.032234_wp, 0.011630_wp, 0.008475_wp, 0.011673_wp, 0.000000_wp, & + & 0.240653_wp, 0.661819_wp, 0.522306_wp, 0.141630_wp, 0.052456_wp, 0.184547_wp, 0.040837_wp, 0.000000_wp, & + & 0.318136_wp, 0.564721_wp, & + & 0.523194_wp, 0.767449_wp, 0.620122_wp, 0.390227_wp, 0.237571_wp, & + & 0.247940_wp, 0.249589_wp, 0.117864_wp, 0.325725_wp, 0.197183_wp, & + & 0.264489_wp, 0.180375_wp, 0.111262_wp, 0.147239_wp, 0.081747_wp, 0.000000_wp] + + real(wp), parameter :: emiss_hf_pobdzvp(*) = [ & ! FIXME: https://github.com/grimme-lab/gcp/issues/22 + & 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & + & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp] + + real(wp), parameter :: emiss_hf_2gcore(*) = [real(wp):: & ! only HCNOF yet + & 0.000539,0.000000,& + & 0.000000,0.000000,0.000000,0.173663,0.269952,0.364341,0.384923,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,& + & 0.000000,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000] + + real(wp), parameter :: emiss_hf_def1tzvp(*) = [real(wp):: & ! org + & 0.007577,0.003312,& + & 0.136371,0.011163,0.017129,0.008140,0.005826,0.006777,0.007108,0.008132,& + & 0.134992,0.147417,0.085253,0.054238,0.033790,0.032862,0.029038,0.026555,& + & 0.141595,0.207980,& + & 0.223252,0.193038,0.167892,0.148726,0.140473,0.130220,0.121166,0.113839,0.121855,0.107138,& + & 0.105637,0.086639,0.075084,0.075089,0.070868,0.068706] + + real(wp), parameter :: emiss_hf_def2mtzvpp(*) = [real(wp):: & !SG + & 0.027000,0.000000,& + & 0.000000,0.000000,0.200000,0.020000,0.180000,0.080000,0.070000,0.065000,& + & 0.000000,0.000000,0.000000,0.200000,0.600000,0.600000,0.600000,0.300000,& + & 0.000000,0.000000,& + & 0.300000,0.300000,0.300000,0.300000,0.300000,0.300000,0.300000,0.300000,0.300000,0.300000,& + & 0.300000,0.300000,0.300000,0.300000,0.300000,0.000000] + + real(wp), parameter :: emiss_hf_2g(*) = [real(wp):: & !no ne, ar ecp + & 0.0539181,0.161846,& + & 0.1581960,0.214318,0.808955,0.470398,0.724457,1.260960,2.014430,0.000000,& + & 0.3072290,0.265975,0.354139,0.307818,0.356962,0.461661,0.588346,0.000000,& + & 0.000000,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000] + + real(wp), parameter :: emiss_hf_ccdz(*) = [real(wp):: & + & 0.007907,0.008287,& + & 0.047380,0.014240,0.022133,0.014999,0.018148,0.028240,0.042261,0.061485,& + & 0.073185,0.056218,0.082660,0.052975,0.033874,0.034056,0.031433,0.030034,& + & 0.000000,0.078016,& !no k cc-pVDZ Basis + & 0.036885,0.038540,0.036474,0.036061,0.030289,0.027959,0.025177,0.022709,0.027386,0.015816,& + & 0.135176,0.115515,0.102761,0.102967,0.097891,0.097331] + + real(wp), parameter :: emiss_hf_accdz(*) = [real(wp):: & !for li,be,na,mg,k-zn energy below def2-QZVPD reference + & 0.001183,0.005948,& + & 0.000000,0.000000,0.005269,0.006380,0.011700,0.021199,0.034160,0.051481,& + & 0.000000,0.000000,0.016018,0.009268,0.010076,0.015153,0.016889,0.018563,& + & 0.000000,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,& + & 0.069963,0.065687,0.072944,0.077585,0.078777,0.080746] + + real(wp), parameter :: emiss_hf_dzp(*) = [real(wp):: & + & 0.008107,0.008045,& + & 0.136751,0.016929,0.026729,0.021682,0.027391,0.040841,0.058747,0.082680,& + & 0.153286,0.162296,0.102704,0.073144,0.056217,0.061333,0.065045,0.071398,& + & 0.145642,0.212865,& + & 0.232821,0.204796,0.182933,0.169554,0.164701,0.160112,0.157723,0.158037,0.179104,0.169782,& + & 0.159396,0.140611,0.129645,0.132664,0.132121,0.134081] + + real(wp), parameter :: emiss_hf_hsv(*) = [real(wp):: & + & 0.030224,0.028324,& + & 0.125379,0.064094,0.059751,0.079387,0.108929,0.167264,0.245786,0.347818,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,& + & 0.000000,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,& + & 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000] + + real(wp), parameter :: emiss_hf_dz(*) = [real(wp):: & + & 0.009037,0.008843,& + & 0.198254,0.000000,0.026921,0.021817,0.027458,0.041391,0.059495,0.083286,& + & 0.268608,0.202374,0.104146,0.075686,0.057826,0.065300,0.069912,0.076845,& + & 0.296046,0.370399,& + & 0.349482,0.302284,0.267639,0.244306,0.232237,0.221488,0.214153,0.032694,0.226865,0.213902,& + & 0.172296,0.155496,0.143646,0.149642,0.149871,0.151705] + + real(wp), parameter :: emiss_hf_msvp(*) = [real(wp):: & !H-Kr modified Ahlrichs DZ, supplemented by def2-SV(P) + & 0.000000_wp,0.000000_wp,& !RG,H set to zero, F adjusted empirically, Be corrected due to ROHF problems + & 0.107750_wp,0.020000_wp,0.026850_wp,0.021740_wp,0.027250_wp,0.039930_wp,0.030000_wp,0.000000_wp,& + & 0.153290_wp,0.162300_wp,0.102700_wp,0.073140_wp,0.056220_wp,0.061330_wp,0.065040_wp,0.000000_wp,& + & 0.200960_wp,0.299660_wp,& + & 0.325990_wp,0.305490_wp,0.291720_wp,0.293800_wp,0.291790_wp,0.296730_wp,0.304600_wp,0.242040_wp,0.354190_wp,0.350720_wp,& + & 0.350020_wp,0.345780_wp,0.349530_wp,0.367310_wp,0.382010_wp,0.000000_wp] + + real(wp), parameter :: emiss_hf_pbeh3c(*) = [ & + & emiss_hf_msvp(1:18), emiss_hf_dzp(19:35), 0.0_wp] + + ! ********************* + ! * nr. of basis fkt * + ! ********************* + integer, parameter :: nbas_sv(*) = [ & + & 2, 2, & + & 3, 3, 9, 9, 9, 9, 9, 9, & + & 7, 7, 13, 13, 13, 13, 13, 13, & + & 11, 11, & + & 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, & + & 27, 27, 27, 27, 27, 27] + integer, parameter :: nbas_minis(*) = [ & + & 1, 1, & + & 2, 2, 5, 5, 5, 5, 5, 5, & + & 6, 6, 9, 9, 9, 9, 9, 9, & + & 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, & + & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 0, 0] + integer, parameter :: nbas_631gd(*) = [ & + & 2, 5, & + & 14, 14, 14, 14, 14, 14, 14, 14, & + & 18, 18, 18, 18, 18, 18, 18, 18, & + & 22, 22, & + & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 32, 0] + integer, parameter :: nbas_svp(*) = [ & + & 5, 5, & + & 9, 9, 14, 14, 14, 14, 14, 14, & + & 15, 18, 18, 18, 18, 18, 18, 18, & + & 24, 24, & + & 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, & + & 32, 32, 32, 32, 32, 32] + integer, parameter :: nbas_sv_p(*) = [ & + & nbas_sv(1), nbas_svp(2:)] + integer, parameter :: nbas_svx(*) = [ & + & nbas_sv(1), nbas_svp(2:5), nbas_sv(6), nbas_svp(7:)] + integer, parameter :: nbas_svp_old(*) = [ & + & 5, 5, & + & 6, 9, 14, 14, 14, 14, 14, 14, & + & 10, 10, 18, 18, 18, 18, 18, 18, & + & 14, 24, & + & 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, & + & 32, 32, 32, 32, 32, 32] + integer, parameter :: nbas_tz(*) = [ & + & 6, 6, & + & 14, 19, 31, 31, 31, 31, 31, 31, & + & 32, 32, 37, 37, 37, 37, 37, 37, & + & 33, 36, & + & 45, 45, 45, 45, 45, 45, 45, 45, 45, 48, & + & 48, 48, 48, 48, 48, 48] + integer, parameter :: nbas_def2mtzvp(*) = [ & + & 3, 3, & + & 8, 11, 19, 19, 19, 24, 19, 19, & + & 14, 14, 22, 22, 22, 22, 22, 22, & + & 18, 28, & + & 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, & + & 36, 36, 36, 36, 36, 36] !def2-mTZVP + integer, parameter :: nbas_def2mtzvpp(*) = [ & + & 5, 5, & + & 9, 11, 19, 19, 24, 24, 24, 24, & + & 14, 14, 27, 27, 27, 27, 27, 27, & + & 18, 28, & + & 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, & + & 36, 36, 36, 36, 36, 36] !def2-mTZVPP + integer, parameter :: nbas_vmb(*) = [ & + & 1, 1, & + & 2, 2, 4, 4, 4, 4, 4, 4, & + & 1, 1, 4, 4, 4, 4, 4, 4, & + & 0, 0, & + & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 0, 0] ! minimal basis set with ECPs + integer, parameter :: nbas_minisd(*) = [ & + & 0, 0, & + & 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 14, 14, 14, 14, 14, 14, & + & 0, 0, & + & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 0, 0] + integer, parameter :: nbas_hsv(*) = [ & ! FIXME: https://github.com/grimme-lab/gcp/issues/23 + & 0, 0, & + & 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, & + & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 0, 0] + integer, parameter :: nbas_minix(*) = [ & + & nbas_minis(1:2), 5, 5, nbas_minis(5:10), 9, 9, & + & nbas_minisd(13:18), nbas_sv(19:30), nbas_svp(31:36)] + integer, parameter :: nbas_lanl2(*) = [ & + & nbas_631gd(1:20), & + & 22, 22, 22, 22, 22, 22, 22, 22, 22, 18, & + & nbas_631gd(31:)] ! Sc-Zn LANL2DZ + integer, parameter :: nbas_pobtz(*) = [ & ! H-Kr no RG + & 6, 0, & + & 7, 7, 18, 18, 18, 18, 18, 0, & + & 19, 19, 22, 22, 22, 22, 22, 0, & + & 23, 23, & + & 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, & + & 41, 41, 41, 41, 41, 0 ] + + integer, parameter :: nbas_pobdzvp(*) = [ & ! H-KR, no RG + & 5, 0, & + & 6, 6, 14, 14, 14, 14, 14, 0, & + & 15, 18, 18, 18, 18, 18, 18, 0, & + & 19, 19, & + & 31, 31, 31, 31, 31, 31, 31, 31, 31, 31,& + & 32, 32, 32, 32, 32, 0] + + integer, parameter :: nbas_2gcore(*) = [ & ! Only HCNOF yet + & 1, 0, & + & 5, 5, 5, 5, 5, 5, 5, 5, & + & 9, 9, 14, 14, 14, 14, 14, 14, & + & 0, 0,& + & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 0, 0] + + integer, parameter :: nbas_2g(*) = [ & + & 1, 1, & + & 5, 5, 5, 5, 5, 5, 5, 5, & + & 9, 9, 14, 14, 14, 14, 14, 14, & + & 0, 0, & + & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + & 0, 0, 0, 0, 0, 0] + + integer, parameter :: nbas_def1tzvp(*) = [ & ! FIXME: https://github.com/grimme-lab/gcp/issues/24 + & 6, 6, & + & 8, 11, 19, 19, 19, 19, 19, 19, & + & 14, 14, 22, 22, 22, 22, 22, 22, & + & 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, & + & 18, 28,& + & 36, 36, 36, 36, 36, 36] + + integer, parameter :: nbas_ccdz(*) = [ & + & 5, 5, & + & 14, 14, 14, 14, 14, 14, 14, 14, & + & 18, 18, 18, 18, 18, 18, 18, 18, & + & 0,27,& + & 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, & + & 27, 27, 27, 27, 27, 27] + + integer, parameter :: nbas_accdz(*) = [ & + & 9, 9,& + & 23, 23, 23, 23, 23, 23, 23, 23,& + & 27, 27, 27, 27, 27, 27, 27, 27,& + & 0, 0,& + & 59, 59, 59, 59, 59, 59, 59, 59, 59, 59,& + & 36, 36, 36, 36, 36, 36] + + integer, parameter :: nbas_dzp(*) = [ & + & 5, 5,& + & 7, 10, 15, 15, 15, 15, 15, 15,& + & 15, 15, 23, 23, 23, 23, 23, 23,& + & 26, 36,& + & 41, 41, 41, 41, 41, 41, 41, 41, 41, 41,& + & 41, 41, 41, 41, 41, 41] + + integer, parameter :: nbas_dz(*) = [ & + & 2, 2, & + & 4, 0, 10, 10, 10, 10, 10, 10, & + & 12, 12, 18, 18, 18, 18, 18, 18, & + & 23, 23, & + & 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, & + & 36, 36, 36, 36, 36, 36] + + integer, parameter :: nbas_msvp(*) = [ & ! modified Ahlrichs DZ, supplemented by def2-SV(P) + & 2, 2,& + & 10, 10, 15, 15, 15, 15, 15, 15,& + & 15, 18, 18, 18, 18, 18, 18, 18,& + & 24, 24, & + & 31, 31, 31, 31, 31, 31, 31, 31, 31, 31,& + & 32, 32, 32, 32, 32, 32] + + real(wp), parameter :: slater_s(*) = [real(wp) :: & + & 1.2000,1.6469,0.6534,1.0365,1.3990,1.7210,2.0348,2.2399,2.5644,2.8812,& + & 0.8675,1.1935,1.5143,1.7580,1.9860,2.1362,2.3617,2.5796,0.9362,1.2112,& + & 1.2870,1.3416,1.3570,1.3804,1.4761,1.5465,1.5650,1.5532,1.5781,1.7778,& + & 2.0675,2.2702,2.4546,2.5680,2.7523,2.9299] + real(wp), parameter :: slater_p(*) = [real(wp) :: & + & 0.0000,0.0000,0.5305,0.8994,1.2685,1.6105,1.9398,2.0477,2.4022,2.7421,& + & 0.6148,0.8809,1.1660,1.4337,1.6755,1.7721,2.0176,2.2501,0.6914,0.9329,& + & 0.9828,1.0104,0.9947,0.9784,1.0641,1.1114,1.1001,1.0594,1.0527,1.2448,& + & 1.5073,1.7680,1.9819,2.0548,2.2652,2.4617] + real(wp), parameter :: slater_d(*) = [real(wp) :: & + & 0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,& + & 0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,& + & 2.4341,2.6439,2.7809,2.9775,3.2208,3.4537,3.6023,3.7017,3.8962,2.0477,& + & 2.4022,2.7421,0.6148,0.8809,1.1660,1.4337] + real(wp), parameter :: slater_exp(*) = [ & + slater_s(1:2), (slater_s(3:20) + slater_p(3:20))/2, & + (slater_s(21:30) + slater_d(21:30) + slater_d(21:30))/3, slater_s(31:36)] + +contains + +subroutine get_gcp_param(param, mol, method, basis, eta) + type(gcp_param), intent(out) :: param + type(structure_type), intent(in) :: mol + character(len=*), intent(in), optional :: method + character(len=*), intent(in), optional :: basis + real(wp), intent(in), optional :: eta + + real(wp) :: eta_, eta_spec + real(wp), allocatable :: nbas(:) + integer :: isp, izp, jsp, basis_id, method_id + integer, allocatable :: nel(:) + logical :: valence_minimal_basis + + method_id = p_unknown_method + if (present(method)) then + method_id = get_method_id(method) + basis_id = get_basis_id(method) + else + basis_id = p_unknown_bas + end if + if (present(basis)) basis_id = get_basis_id(basis) + if (basis_id == p_unknown_bas .and. method_id == p_unknown_method) return + + eta_ = 0.0_wp + if (present(eta)) eta_ = eta + eta_spec = 0.0_wp + + allocate(param%zeff(mol%nid)) + do isp = 1, mol%nid + izp = mol%num(isp) + param%zeff(isp) = effective_atomic_number(izp) + end do + + allocate(param%rvdw(mol%nid, mol%nid)) + do isp = 1, mol%nid + do jsp = 1, mol%nid + param%rvdw(isp, jsp) = get_vdw_rad(param%zeff(isp), param%zeff(jsp)) + end do + end do + + allocate(param%rvdw_srb(mol%nid, mol%nid)) + do isp = 1, mol%nid + do jsp = 1, mol%nid + param%rvdw_srb(isp, jsp) = get_vdw_rad(mol%num(isp), mol%num(jsp)) + end do + end do + + valence_minimal_basis = basis_id == p_vmb_bas + + allocate(nel(mol%nid)) + do isp = 1, mol%nid + izp = param%zeff(isp) + nel(isp) = number_of_electrons(izp, valence_minimal_basis) + end do + + select case(basis_id) + case(p_vmb_bas) + param%emiss = emiss_hf_vmb(param%zeff) + nbas = nbas_vmb(param%zeff) + + case(p_sv_bas) + param%emiss = emiss_hf_sv(param%zeff) + nbas = nbas_sv(param%zeff) + + select case(method_id) + case (p_hf_method) ! RMS=0.3218975 + param%sigma = 0.1724_wp + eta_ = 1.2804_wp + param%alpha = 0.8568_wp + param%beta = 1.2342_wp + case (p_hyb_method, p_b3lyp_method, p_pw6b95_method) ! RMS= 0.557 + param%sigma = 0.4048_wp + eta_ = 1.1626_wp + param%alpha = 0.8652_wp + param%beta = 1.2375_wp + case (p_gga_method, p_tpss_method, p_blyp_method) ! RMS = 0.6652 + param%sigma = 0.2727_wp + eta_ = 1.4022_wp + param%alpha = 0.8055_wp + param%beta = 1.3000_wp + end select + + case(p_sv_p_bas) + param%emiss = emiss_hf_sv_p(param%zeff) + nbas = nbas_sv_p(param%zeff) + + if (method_id == p_hf_method) then !RMS=0.3502 + param%sigma = 0.1373_wp + eta_ = 1.4271_wp + param%alpha = 0.8141_wp + param%beta = 1.2760_wp + elseif (is_dft_method(method_id)) then ! RMS= 0.57 ! def2-SV(P) + param%sigma = 0.2424_wp + eta_ = 1.2371_wp + param%alpha = 0.6076_wp + param%beta = 1.4078_wp + end if + + case(p_svx_bas) ! RMS= ! def2-SV(P/h,c) = SV at h,c + param%emiss = emiss_hf_svx(param%zeff) + nbas = nbas_svx(param%zeff) + + if (is_dft_method(method_id)) then ! RMS= 0.56 ! def2-SV(P/h,c) = SV at h,c + param%sigma = 0.1861_wp + eta_ = 1.3200_wp + param%alpha = 0.6171_wp + param%beta = 1.4019_wp + end if + + case(p_svp_bas) + param%emiss = emiss_hf_svp(param%zeff) + nbas = nbas_svp(param%zeff) + + if (method_id == p_hf_method) then ! RMS=0.4065 + param%sigma = 0.2054_wp + eta_ = 1.3157_wp + param%alpha = 0.8136_wp + param%beta = 1.2572_wp + elseif (method_id == p_tpss_method) then ! RMS= 0.618 + param%sigma = 0.6647_wp + eta_ = 1.3306_wp + param%alpha = 1.0792_wp + param%beta = 1.1651_wp + elseif (method_id == p_pw6b95_method) then ! RMS = 0.58312 + param%sigma = 0.3098_wp + eta_ = 1.2373_wp + param%alpha = 0.6896_wp + param%beta = 1.3347_wp + elseif (is_hyb_method(method_id)) then ! RMS=0.6498 + param%sigma = 0.2990_wp + eta_ = 1.2605_wp + param%alpha = 0.6438_wp + param%beta = 1.3694_wp + elseif (is_gga_method(method_id)) then ! RMS= + param%sigma = 0.6823_wp + eta_ = 1.2491_wp + param%alpha = 0.8225_wp + param%beta = 1.2811_wp + end if + + case(p_svp_old_bas) + param%emiss = emiss_hf_svp_old(param%zeff) + nbas = nbas_svp_old(param%zeff) + + if (method_id == p_hf_method) then ! RMS=0.4065 + param%sigma = 0.2054_wp + eta_ = 1.3157_wp + param%alpha = 0.8136_wp + param%beta = 1.2572_wp + elseif (is_dft_method(method_id)) then ! RMS=0.6498 + param%sigma = 0.2990_wp + eta_ = 1.2605_wp + param%alpha = 0.6438_wp + param%beta = 1.3694_wp + end if + + case(p_minis_bas) + param%emiss = emiss_hf_minis(param%zeff) + nbas = nbas_minis(param%zeff) + + if (method_id == p_hf_method) then ! RMS= 0.3040 + param%sigma = 0.1290_wp + eta_ = 1.1526_wp + param%alpha = 1.1549_wp + param%beta = 1.1763_wp + elseif (method_id == p_tpss_method) then ! RMS= + param%sigma = 0.22982_wp + eta_ = 1.35401_wp + param%alpha = 1.47633_wp + param%beta = 1.11300_wp + elseif (method_id == p_pw6b95_method) then ! RMS = 0.3279929 + param%sigma = 0.21054_wp + eta_ = 1.25458_wp + param%alpha = 1.35003_wp + param%beta = 1.14061_wp + elseif (is_gga_method(method_id)) then ! RMS= 0.3462 + param%sigma = 0.1566_wp + eta_ = 1.0271_wp + param%alpha = 1.0732_wp + param%beta = 1.1968_wp + elseif (is_hyb_method(method_id)) then ! RMS= 0.3400 + param%sigma = 0.2059_wp + eta_ = 0.9722_wp + param%alpha = 1.1961_wp + param%beta = 1.1456_wp + end if + + case(p_631gd_bas) + param%emiss = emiss_hf_631gd(param%zeff) + nbas = nbas_631gd(param%zeff) + + if (method_id == p_hf_method) then ! RMS= 0.40476 + param%sigma = 0.2048_wp + eta_ = 1.5652_wp + param%alpha = 0.9447_wp + param%beta = 1.2100_wp + elseif (is_dft_method(method_id)) then ! RMS= 0.47856 + param%sigma = 0.3405_wp + eta_ = 1.6127_wp + param%alpha = 0.8589_wp + param%beta = 1.2830_wp + end if + + case(p_tz_bas) + param%emiss = emiss_hf_tz(param%zeff) + nbas = nbas_tz(param%zeff) + + if (method_id == p_hf_method) then ! RMS= 0.1150 + param%sigma = 0.3127_wp + eta_ = 1.9914_wp + param%alpha = 1.0216_wp + param%beta = 1.2833_wp + elseif (is_hyb_method(method_id)) then ! RMS=0.19648 + param%sigma = 0.2905_wp + eta_ = 2.2495_wp + param%alpha = 0.8120_wp + param%beta = 1.4412_wp + elseif (is_gga_method(method_id)) then !RMS = 0.21408 + param%sigma = 0.1182_wp + eta_ = 1.0631_wp + param%alpha = 1.0510_wp + param%beta = 1.1287_wp + end if + + case(p_deftzvp_bas) + param%emiss = emiss_hf_def1tzvp(param%zeff) + nbas = nbas_def1tzvp(param%zeff) + + if (method_id == p_hf_method) then ! RMS=0.209 + param%sigma = 0.2600_wp + eta_ = 2.2448_wp + param%alpha = 0.7998_wp + param%beta = 1.4381_wp + elseif (is_dft_method(method_id)) then ! RMS=0.1817 + param%sigma = 0.2393_wp + eta_ = 2.2247_wp + param%alpha = 0.8185_wp + param%beta = 1.4298_wp + end if + + case(p_ccdz_bas) + param%emiss = emiss_hf_ccdz(param%zeff) + nbas = nbas_ccdz(param%zeff) + + if (method_id == p_hf_method) then ! RMS=0.4968 + param%sigma = 0.4416_wp + eta_ = 1.5185_wp + param%alpha = 0.6902_wp + param%beta = 1.3713_wp + elseif (is_dft_method(method_id)) then ! RMS=0.7610 + param%sigma = 0.5383_wp + eta_ = 1.6482_wp + param%alpha = 0.6230_wp + param%beta = 1.4523_wp + end if + + case(p_accdz_bas) + param%emiss = emiss_hf_accdz(param%zeff) + nbas = nbas_accdz(param%zeff) + + if (method_id == p_hf_method) then !RMS=0.2222 + param%sigma = 0.0748_wp + eta_ = 0.0663_wp + param%alpha = 0.3811_wp + param%beta = 1.0155_wp + elseif (is_dft_method(method_id)) then ! RMS=0.1840 + param%sigma = 0.1465_wp + eta_ = 0.0500_wp + param%alpha = 0.6003_wp + param%beta = 0.8761_wp + end if + + case(p_pobtz_bas) + param%emiss = emiss_hf_pobtz(param%zeff) + nbas = nbas_pobtz(param%zeff) + + if (is_dft_method(method_id)) then + param%sigma = 0.1300_wp + eta_ = 1.3743_wp + param%alpha = 0.4792_wp + param%beta = 1.3962_wp + end if + + ! case(p_pobdzvp_bas) + ! param%emiss = emiss_hf_pobdzvp(param%zeff) + ! nbas = nbas_pobdzvp(param%zeff) + + case(p_minix_bas) + param%emiss = emiss_hf_minix(param%zeff) + nbas = nbas_minix(param%zeff) + + if (method_id == p_hf_method .or. method_id == p_hf3c_method) then + param%sigma = 0.1290_wp + eta_ = 1.1526_wp + param%alpha = 1.1549_wp + param%beta = 1.1763_wp + param%base = method_id == p_hf3c_method + if (param%base) then + param%rscal = 0.7_wp + param%qscal = 0.03_wp + end if + elseif (is_dft_method(method_id)) then + param%sigma = 0.2059_wp + eta_ = 0.9722_wp + param%alpha = 1.1961_wp + param%beta = 1.1456_wp + end if + + case(p_gcore_bas) + param%emiss = emiss_hf_2gcore(param%zeff) + nbas = nbas_2gcore(param%zeff) + + case(p_2g_bas) + param%emiss = emiss_hf_2g(param%zeff) + nbas = nbas_2g(param%zeff) + + if (method_id == p_hf_method) then + param%sigma = 0.2461_wp + eta_ = 1.1616_wp + param%alpha = 0.7335_wp + param%beta = 1.4709_wp + end if + + case(p_dzp_bas) + param%emiss = emiss_hf_dzp(param%zeff) + nbas = nbas_dzp(param%zeff) + + if (method_id == p_hf_method) then !RMS=0.4571 + param%sigma = 0.1443_wp + eta_ = 1.4547_wp + param%alpha = 0.3711_wp + param%beta = 1.6300_wp + elseif (is_dft_method(method_id)) then !RMS=0.7184 + param%sigma = 0.2687_wp + eta_ = 1.4634_wp + param%alpha = 0.3513_wp + param%beta = 1.6880_wp + end if + + ! case(p_hsv_bas) + ! param%emiss = emiss_hf_hsv(param%zeff) + ! nbas = nbas_hsv(param%zeff) + + case(p_dz_bas) + param%emiss = emiss_hf_dz(param%zeff) + nbas = nbas_dz(param%zeff) + + if (method_id == p_hf_method) then !RMS=0.3754 + param%sigma = 0.1059_wp + eta_ = 1.4554_wp + param%alpha = 0.3711_wp + param%beta = 1.6342_wp + elseif (is_dft_method(method_id)) then + param%sigma = 0.2687_wp + eta_ = 1.4634_wp + param%alpha = 0.3513_wp + param%beta = 1.6880_wp + end if + + case(p_msvp_bas) + param%emiss = emiss_hf_msvp(param%zeff) + nbas = nbas_msvp(param%zeff) + + case(p_def2mtzvp_bas) + param%emiss = emiss_hf_def2mtzvp(param%zeff) + nbas = nbas_def2mtzvp(param%zeff) + + if (method_id == p_b3pbe3c_method) then + param%sigma = 1.0000_wp + eta_ = 2.98561_wp + param%alpha = 0.3011_wp + param%beta = 2.4405_wp + end if + + case(p_pbeh3c_bas) + param%emiss = emiss_hf_pbeh3c(param%zeff) + nbas = nbas_msvp(param%zeff) + + if (method_id == p_pbeh3c_method) then + param%sigma = 1.00000_wp + eta_ = 1.32492_wp + param%alpha = 0.27649_wp + param%beta = 1.95600_wp + param%damp = .true. + elseif (method_id == p_hse3c_method) then + param%sigma = 1.00000_wp + eta_ = 1.32378_wp + param%alpha = 0.28314_wp + param%beta = 1.94527_wp + param%damp = .true. + end if + + case(p_lanl2_bas) + param%emiss = emiss_hf_lanl2(param%zeff) + nbas = nbas_lanl2(param%zeff) + + if (is_dft_method(method_id)) then + param%sigma = 0.3405_wp + eta_ = 1.6127_wp + param%alpha = 0.8589_wp + param%beta = 1.2830_wp + end if + + case(p_def2mtzvpp_bas) + param%emiss = emiss_hf_def2mtzvpp(param%zeff) + nbas = nbas_def2mtzvpp(param%zeff) + + if (method_id == p_r2scan3c_method) then + param%sigma = 1.0000_wp + eta_ = 1.3150_wp + eta_spec = 1.15_wp + param%alpha = 0.9410_wp + param%beta = 1.4636_wp + param%damp = .true. + end if + + end select + if (allocated(nbas)) then + param%xv = nbas - 0.5_wp * nel + if (basis_id == p_def2mtzvpp_bas) then + param%xv(:) = 1.0_wp + where(param%zeff == 6) param%xv = 3.0_wp + where(param%zeff == 7) param%xv = 0.5_wp + where(param%zeff == 8) param%xv = 0.5_wp + end if + end if + + if (eta_ > 0.0_wp) then + param%slater = eta_ * slater_exp(param%zeff) + if (eta_spec > 0.0_wp) then + where(param%zeff > 10) param%slater = eta_spec * param%slater + end if + end if + + if (method_id == p_b973c_method) then + param%srb=.true. + param%rscal=10.00_wp + param%qscal=0.08_wp + end if +end subroutine get_gcp_param + +pure function get_method_id(method) result(id) + character(len=*), intent(in) :: method + integer :: id + select case(method) + case default + id = p_unknown_method + case('hf') + id = p_hf_method + case('dft') + id = p_dft_method + case('b3lyp') + id = p_b3lyp_method + case('blyp') + id = p_blyp_method + case('gga') + id = p_gga_method + case('tpss') + id = p_tpss_method + case('pw6b95') + id = p_pw6b95_method + case('b973c') + id = p_b973c_method + case('r2scan3c') + id = p_r2scan3c_method + case('pbe') + id = p_pbe_method + case('pbeh3c') + id = p_pbeh3c_method + case('hse3c') + id = p_hse3c_method + case('hf3c') + id = p_hf3c_method + case('b3pbe3c') + id = p_b3pbe3c_method + end select +end function get_method_id + +pure function is_dft_method(method_id) result(res) + integer, intent(in) :: method_id + logical :: res + res = method_id == p_dft_method .or. method_id == p_b3lyp_method .or. & + method_id == p_blyp_method .or. method_id == p_gga_method .or. & + method_id == p_tpss_method .or. method_id == p_pw6b95_method +end function is_dft_method + +pure function is_hyb_method(method_id) result(res) + integer, intent(in) :: method_id + logical :: res + res = method_id == p_b3lyp_method .or. method_id == p_pw6b95_method +end function is_hyb_method + +pure function is_gga_method(method_id) result(res) + integer, intent(in) :: method_id + logical :: res + res = method_id == p_gga_method .or. method_id == p_tpss_method .or. & + method_id == p_blyp_method +end function is_gga_method + +pure function get_basis_id(basis) result(id) + character(len=*), intent(in) :: basis + integer :: id + select case(basis) + case default + id = p_unknown_bas + case('sv') + id = p_sv_bas + case('sv(p)','def2sv(p)', 'sv_p', 'def2sv_p') + id = p_sv_p_bas + case ('svx') ! RMS= ! def2-SV(P/h,c) = SV at h,c + id = p_svx_bas + case('svp') + id = p_svp_bas + case('minis') + id = p_minis_bas + case('631gd', '631gs') + id = p_631gd_bas + case('tz', 'def2tzvp') + id = p_tz_bas + case('deftzvp', 'def1tzvp') + id = p_deftzvp_bas + case('ccdz', 'ccpvdz') + id = p_ccdz_bas + case('accdz', 'augccpvdz', 'accpvdz') + id = p_accdz_bas + case('pobtz', 'pobtzvp') + id = p_pobtz_bas + ! case('pobdzvp') + ! id = p_pobdzvp_bas + case ('minix', 'hf3c') + id = p_minix_bas + case('gcore') + id = p_gcore_bas + case('2g', 'twog', 'fitg') + id = p_2g_bas + case('dzp') + id = p_dzp_bas + ! case('hsv') + ! id = p_hsv_bas + case('lanl') + id = p_lanl2_bas + case('dz') + id = p_dz_bas + case('msvp', 'def2msvp') + id = p_msvp_bas + case('pbeh3c', 'hse3c') + id = p_pbeh3c_bas + case('mtzvp', 'def2mtzvp') + id = p_def2mtzvp_bas + case('mtzvpp', 'def2mtzvpp', 'r2scan3c') + id = p_def2mtzvpp_bas + end select +end function get_basis_id + +pure function effective_atomic_number(number) result(zeff) + integer, intent(in) :: number + integer :: zeff + + select case (number) + case default + zeff = number + case(37:54) + zeff = number - 18 + case(55:57) + zeff = number - 18*2 + case(58:71, 90:94) + zeff = 21 + case(72:89) + zeff = number - 2*18 - 14 + end select +end function effective_atomic_number + +pure function number_of_electrons(number, valence_minimal_basis) result(nel) + integer, intent(in) :: number + logical, intent(in) :: valence_minimal_basis + integer :: nel + + if (valence_minimal_basis) then + select case(number) + case default + nel = number + case(5:10) + nel = number - 2 + case(11:18) + nel = number - 10 + end select + else + nel = number + end if +end function number_of_electrons + +endmodule dftd3_gcp_param \ No newline at end of file diff --git a/src/dftd3/meson.build b/src/dftd3/meson.build index 97b81163..1cd9cae6 100644 --- a/src/dftd3/meson.build +++ b/src/dftd3/meson.build @@ -16,6 +16,7 @@ subdir('damping') subdir('data') +subdir('gcp') srcs += files( 'citation.f90', @@ -23,6 +24,7 @@ srcs += files( 'damping.f90', 'data.f90', 'disp.f90', + 'gcp.f90', 'model.f90', 'ncoord.f90', 'output.f90', diff --git a/src/dftd3/output.f90 b/src/dftd3/output.f90 index 1ddbde06..0dc3104c 100644 --- a/src/dftd3/output.f90 +++ b/src/dftd3/output.f90 @@ -24,6 +24,7 @@ module dftd3_output use dftd3_damping_optimizedpower, only : optimizedpower_damping_param use dftd3_damping_rational, only : rational_damping_param use dftd3_damping_zero, only : zero_damping_param + use dftd3_gcp, only : gcp_param use dftd3_model, only : d3_model use dftd3_version, only : get_dftd3_version implicit none @@ -31,7 +32,7 @@ module dftd3_output public :: ascii_atomic_radii, ascii_atomic_references, ascii_system_properties public :: ascii_energy_atom - public :: ascii_results, ascii_damping_param, ascii_pairwise + public :: ascii_results, ascii_damping_param, ascii_pairwise, ascii_gcp_param public :: turbomole_gradient, turbomole_gradlatt public :: json_results, tagged_result @@ -142,7 +143,7 @@ end subroutine ascii_system_properties !> Print atom-resolved dispersion energies -subroutine ascii_energy_atom(unit, mol, energies) +subroutine ascii_energy_atom(unit, mol, energies, label) !> Unit for output integer, intent(in) :: unit @@ -153,9 +154,16 @@ subroutine ascii_energy_atom(unit, mol, energies) !> Atom-resolved dispersion energies real(wp), allocatable, intent(in) :: energies(:) + !> Label for the output + character(len=*), intent(in), optional :: label + integer :: iat, isp + character(len=:), allocatable :: label_ + + label_ = "dispersion" + if (present(label)) label_ = label - write(unit, '(a,":")') "Atom-resolved dispersion energies" + write(unit, '(a,":")') "Atom-resolved "//label_//" energies" write(unit, '(50("-"))') write(unit, '(a6,1x,a4,1x,4x,a15,1x,a15)') "#", "Z", "[Hartree]", "[kcal/mol]" write(unit, '(50("-"))') @@ -170,7 +178,7 @@ subroutine ascii_energy_atom(unit, mol, energies) end subroutine ascii_energy_atom -subroutine ascii_results(unit, mol, energy, gradient, sigma) +subroutine ascii_results(unit, mol, energy, gradient, sigma, label) !> Unit for output integer, intent(in) :: unit @@ -182,14 +190,21 @@ subroutine ascii_results(unit, mol, energy, gradient, sigma) real(wp), intent(in), optional :: gradient(:, :) real(wp), intent(in), optional :: sigma(:, :) + !> Label for the output + character(len=*), intent(in), optional :: label + integer :: iat, isp logical :: grad character(len=1), parameter :: comp(3) = ["x", "y", "z"] + character(len=:), allocatable :: label_ + + label_ = "Dispersion" + if (present(label)) label_ = label grad = present(gradient) .and. present(sigma) write(unit, '(a,":", t25, es20.13, 1x, a)') & - & "Dispersion energy", energy, "Eh" + & label_//" energy", energy, "Eh" write(unit, '(a)') if (grad) then write(unit, '(a,":", t25, es20.13, 1x, a)') & @@ -371,6 +386,61 @@ subroutine ascii_damping_param(unit, param, method) end subroutine ascii_damping_param +subroutine ascii_gcp_param(unit, mol, param, method) + + !> Unit for output + integer, intent(in) :: unit + + !> Molecular structure data + type(structure_type), intent(in) :: mol + + !> Counter-poise parameters + type(gcp_param), intent(in) :: param + + !> Method name + character(len=*), intent(in), optional :: method + + integer :: isp + + write(unit, '(a,":")') "Global counter-poise parameters" + write(unit, '(20("-"))') + if (param%sigma > 0.0_wp .and. param%alpha > 0.0_wp .and. param%beta > 0.0_wp) then + write(unit, '(a6, t10, f10.4)') & + & "sigma", param%sigma, & + & "alpha", param%alpha, & + & "beta", param%beta + end if + if (param%damp) then + write(unit, '(a6, t10, f10.4)') & + & "dscal", param%dmp_scal, & + & "dexpo", param%dmp_exp + end if + if (param%srb) then + write(unit, '(a6, t10, f10.4)') & + & "rscal", param%rscal, & + & "qscal", param%qscal + end if + write(unit, '(20("-"))') + write(unit, '(a)') + + if (allocated(param%emiss) .and. allocated(param%xv) & + & .and. allocated(param%slater)) then + write(unit, '(a,":")') "Atomic counter-poise parameters" + write(unit, '(47("-"))') + write(unit, '(a4,5x,a4,*(1x,a10))') "Z", "Zeff", "Emiss[Eh]", "Virtual", "Slater" + write(unit, '(47("-"))') + do isp = 1, mol%nid + write(unit, '(i4, 1x, a4, i4, *(1x,f10.4))') & + & mol%num(isp), mol%sym(isp), param%zeff(isp), & + & param%emiss(isp), param%xv(isp), param%slater(isp) + end do + write(unit, '(47("-"))') + write(unit, '(a)') + end if + +end subroutine ascii_gcp_param + + subroutine turbomole_gradlatt(mol, fname, energy, sigma, stat) type(structure_type),intent(in) :: mol character(len=*),intent(in) :: fname diff --git a/test/api/api-test.c b/test/api/api-test.c index de977025..5d31bbfb 100644 --- a/test/api/api-test.c +++ b/test/api/api-test.c @@ -28,6 +28,22 @@ show_error(dftd3_error error) printf("[Message] %s\n", message); } +static inline dftd3_structure +get_test_structure(dftd3_error error) +{ + int const natoms = 7; + int const attyp[7] = {6,6,6,1,1,1,1}; + double const coord[21] = + {0.00000000000000, 0.00000000000000,-1.79755622305860, + 0.00000000000000, 0.00000000000000, 0.95338756106749, + 0.00000000000000, 0.00000000000000, 3.22281255790261, + -0.96412815539807,-1.66991895015711,-2.53624948351102, + -0.96412815539807, 1.66991895015711,-2.53624948351102, + 1.92825631079613, 0.00000000000000,-2.53624948351102, + 0.00000000000000, 0.00000000000000, 5.23010455462158}; + return dftd3_new_structure(error, natoms, attyp, coord, NULL, NULL); +} + int test_version (void) { @@ -56,12 +72,136 @@ test_uninitialized_structure (void) dftd3_update_structure(error, mol, xyz, NULL); if (!dftd3_check_error(error)) goto unexpected; + show_error(error); + + dftd3_delete(error); + return 0; + +unexpected: + printf("[Fatal] Unexpected pass for uninitialized-structure test\n"); + dftd3_delete(error); + return 1; +} + +int +test_uninitialized_param (void) +{ + printf("Start test: uninitialized parameters\n"); + dftd3_error error = NULL; + dftd3_structure mol = NULL; + dftd3_model disp = NULL; + dftd3_param param = NULL; + double energy; + + error = dftd3_new_error(); + + dftd3_get_dispersion(error, mol, disp, param, &energy, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + dftd3_get_pairwise_dispersion(error, mol, disp, param, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + mol = get_test_structure(error); + dftd3_get_dispersion(error, mol, disp, param, &energy, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + dftd3_get_pairwise_dispersion(error, mol, disp, param, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + disp = dftd3_new_d3_model(error, mol); + + dftd3_get_dispersion(error, mol, disp, param, &energy, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + dftd3_get_pairwise_dispersion(error, mol, disp, param, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + dftd3_delete(error); + dftd3_delete(disp); + dftd3_delete(mol); + return 0; + +unexpected: + printf("[Fatal] Unexpected pass for uninitialized-parameters test\n"); + dftd3_delete(error); + dftd3_delete(disp); + dftd3_delete(mol); + return 1; +} + +int +test_uninitialized_model (void) +{ + printf("Start test: uninitialized model\n"); + dftd3_error error = NULL; + dftd3_model disp = NULL; + + error = dftd3_new_error(); + + dftd3_set_model_realspace_cutoff(error, disp, 0.0, 0.0, 0.0); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + dftd3_delete(error); + dftd3_delete(disp); + return 0; + +unexpected: + printf("[Fatal] Unexpected pass for uninitialized-model test\n"); + dftd3_delete(error); + dftd3_delete(disp); + return 1; +} + +int +test_uninitialized_gcp (void) +{ + printf("Start test: uninitialized counter-poise\n"); + dftd3_error error = NULL; + dftd3_structure mol = NULL; + dftd3_gcp gcp = NULL; + double energy; + + error = dftd3_new_error(); + + dftd3_set_gcp_realspace_cutoff(error, gcp, 0.0, 0.0); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + dftd3_get_counterpoise(error, mol, gcp, &energy, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + + mol = get_test_structure(error); + + dftd3_get_counterpoise(error, mol, gcp, &energy, NULL, NULL); + if (!dftd3_check_error(error)) goto unexpected; + + show_error(error); + dftd3_delete(error); + dftd3_delete(mol); return 0; unexpected: - printf("[Fatal] Unexpected pass for unititalized-structure test\n"); + printf("[Fatal] Unexpected pass for uninitialized-counter-poise test\n"); dftd3_delete(error); + dftd3_delete(mol); return 1; } @@ -95,17 +235,7 @@ test_invalid_structure (void) } int -test (void) { - int const natoms = 7; - int const attyp[7] = {6,6,6,1,1,1,1}; - double const coord[21] = - {0.00000000000000, 0.00000000000000,-1.79755622305860, - 0.00000000000000, 0.00000000000000, 0.95338756106749, - 0.00000000000000, 0.00000000000000, 3.22281255790261, - -0.96412815539807,-1.66991895015711,-2.53624948351102, - -0.96412815539807, 1.66991895015711,-2.53624948351102, - 1.92825631079613, 0.00000000000000,-2.53624948351102, - 0.00000000000000, 0.00000000000000, 5.23010455462158}; +test_d3 (void) { double energy; double pair_disp2[49]; double pair_disp3[49]; @@ -120,7 +250,7 @@ test (void) { error = dftd3_new_error(); assert(!!error); - mol = dftd3_new_structure(error, natoms, attyp, coord, NULL, NULL); + mol = get_test_structure(error); if (dftd3_check_error(error)) {return 1;}; assert(!!mol); @@ -185,6 +315,58 @@ test (void) { return 0; } +int +test_gcp (void) { + double energy; + double pair_disp2[49]; + double pair_disp3[49]; + double gradient[21]; + double sigma[9]; + + dftd3_error error; + dftd3_structure mol; + dftd3_gcp gcp; + + error = dftd3_new_error(); + assert(!!error); + + mol = get_test_structure(error); + if (dftd3_check_error(error)) {return 1;}; + assert(!!mol); + + gcp = dftd3_load_gcp_param(error, mol, "pbeh3c", NULL); + if (dftd3_check_error(error)) {return 1;} + assert(!!gcp); + + dftd3_get_counterpoise(error, mol, gcp, &energy, NULL, NULL); + if (dftd3_check_error(error)) {return 1;} + dftd3_get_counterpoise(error, mol, gcp, &energy, gradient, sigma); + if (dftd3_check_error(error)) {return 1;} + dftd3_delete(gcp); + + gcp = dftd3_load_gcp_param(error, mol, "b973c", NULL); + if (dftd3_check_error(error)) {return 1;} + assert(!!gcp); + + dftd3_set_gcp_realspace_cutoff(error, gcp, 50.0, 30.0); + if (dftd3_check_error(error)) {return 1;} + + dftd3_get_counterpoise(error, mol, gcp, &energy, NULL, NULL); + if (dftd3_check_error(error)) {return 1;} + dftd3_get_counterpoise(error, mol, gcp, &energy, gradient, sigma); + if (dftd3_check_error(error)) {return 1;} + dftd3_delete(gcp); + + dftd3_delete(mol); + dftd3_delete(error); + + assert(!gcp); + assert(!mol); + assert(!error); + + return 0; +} + int main (void) { @@ -192,7 +374,11 @@ main (void) stat += test_version(); stat += test_uninitialized_error(); stat += test_uninitialized_structure(); + stat += test_uninitialized_model(); + stat += test_uninitialized_param(); + stat += test_uninitialized_gcp(); stat += test_invalid_structure(); - stat += test(); + stat += test_d3(); + stat += test_gcp(); return stat; } diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 85f63e25..c8c02410 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -23,6 +23,7 @@ set( "param" "pairwise" "periodic" + "gcp" "regression" ) set( diff --git a/test/unit/main.f90 b/test/unit/main.f90 index 284a67f7..b6feb51b 100644 --- a/test/unit/main.f90 +++ b/test/unit/main.f90 @@ -20,6 +20,7 @@ program tester use mctc_env_testing, only : run_testsuite, new_testsuite, testsuite_type, & & select_suite, run_selected use test_dftd3, only : collect_dftd3 + use test_gcp, only : collect_gcp use test_model, only : collect_model use test_ncoord, only : collect_ncoord use test_param, only : collect_param @@ -41,6 +42,7 @@ program tester & new_testsuite("param", collect_param), & & new_testsuite("pairwise", collect_pairwise), & & new_testsuite("periodic", collect_periodic), & + & new_testsuite("gcp", collect_gcp), & & new_testsuite("regression", collect_regression) & & ] diff --git a/test/unit/meson.build b/test/unit/meson.build index ca527b68..753b9e4f 100644 --- a/test/unit/meson.build +++ b/test/unit/meson.build @@ -33,6 +33,7 @@ tests = [ 'param', 'pairwise', 'periodic', + 'gcp', 'regression', ] diff --git a/test/unit/test_gcp.f90 b/test/unit/test_gcp.f90 new file mode 100644 index 00000000..384a4324 --- /dev/null +++ b/test/unit/test_gcp.f90 @@ -0,0 +1,1082 @@ +! This file is part of s-dftd3. +! SPDX-Identifier: LGPL-3.0-or-later +! +! s-dftd3 is free software: you can redistribute it and/or modify it under +! the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! s-dftd3 is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with s-dftd3. If not, see . + +module test_gcp + use mctc_env, only : wp + use mctc_env_testing, only : new_unittest, unittest_type, error_type, check, & + & test_failed + use mctc_io, only : structure_type + use mstore, only : get_structure + use dftd3_gcp + use dftd3_cutoff, only : realspace_cutoff + use dftd3_output, only : ascii_gcp_param + implicit none + private + + public :: collect_gcp + + real(wp), parameter :: thr = 100*epsilon(1.0_wp) + real(wp), parameter :: thr2 = sqrt(epsilon(1.0_wp)) * 10 + + +contains + + +!> Collect all exported unit tests +subroutine collect_gcp(testsuite) + + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + ! new_unittest("vmb", test_vmb), & + & new_unittest("sv", test_sv), & + & new_unittest("sv(p)", test_sv_p), & + & new_unittest("svx", test_svx), & + & new_unittest("svp", test_svp), & + & new_unittest("minis", test_minis), & + & new_unittest("minix", test_minix), & + & new_unittest("631gd", test_631gd), & + & new_unittest("def2tzvp", test_tz), & + & new_unittest("def1tzvp", test_deftzvp), & + & new_unittest("ccdz", test_ccdz), & + & new_unittest("accdz", test_accdz), & + & new_unittest("pobtzvp", test_pobtz), & + & new_unittest("2g", test_2g), & + & new_unittest("dz", test_dz), & + & new_unittest("dzp", test_dzp), & + & new_unittest("lanl", test_lanl), & + & new_unittest("msvp", test_msvp), & + & new_unittest("def2mtzvp", test_def2mtzvp), & + ! new_unittest("def2mtzvpp", test_def2mtzvpp), & + & new_unittest("hf3c", test_hf3c), & + & new_unittest("pbeh3c", test_pbeh3c), & + & new_unittest("hse3c", test_hse3c), & + & new_unittest("b973c", test_b973c), & + & new_unittest("r2scan3c", test_r2scan3c), & + & new_unittest("hf/sv", test_hf_sv), & + & new_unittest("hybrid/sv", test_hyb_sv), & + & new_unittest("gga/sv", test_gga_sv), & + & new_unittest("hf/sv(p)", test_hf_sv_p), & + & new_unittest("dft/sv(p)", test_dft_sv_p), & + & new_unittest("dft/svx", test_dft_svx), & + & new_unittest("hf/svp", test_hf_svp), & + & new_unittest("tpss/svp", test_tpss_svp), & + & new_unittest("pw6b95/svp", test_pw6b95_svp), & + & new_unittest("hybrid/svp", test_hyb_svp), & + & new_unittest("gga/svp", test_gga_svp), & + & new_unittest("hf/minis", test_hf_minis), & + & new_unittest("tpss/minis", test_tpss_minis), & + & new_unittest("pw6b95/minis", test_pw6b95_minis), & + & new_unittest("hybrid/minis", test_hyb_minis), & + & new_unittest("gga/minis", test_gga_minis), & + & new_unittest("hf/631gd", test_hf_631gd), & + & new_unittest("dft/631gd", test_dft_631gd), & + & new_unittest("hf/tz", test_hf_tz), & + & new_unittest("hybrid/tz", test_hyb_tz), & + & new_unittest("gga/tz", test_gga_tz), & + & new_unittest("hf/deftzvp", test_hf_deftzvp), & + & new_unittest("dft/deftzvp", test_dft_deftzvp), & + & new_unittest("hf/ccdz", test_hf_ccdz), & + & new_unittest("dft/ccdz", test_dft_ccdz), & + & new_unittest("hf/accdz", test_hf_accdz), & + & new_unittest("dft/accdz", test_dft_accdz), & + & new_unittest("dft/pobtz", test_dft_pobtz), & + & new_unittest("hf/minix", test_hf_minix), & + & new_unittest("dft/minix", test_dft_minix), & + & new_unittest("hf/2g", test_hf_2g), & + & new_unittest("hf/dzp", test_hf_dzp), & + & new_unittest("dft/dzp", test_dft_dzp), & + & new_unittest("hf/dz", test_hf_dz), & + & new_unittest("dft/dz", test_dft_dz), & + & new_unittest("b3pbe3c/def2mtzvp", test_b3pbe3c_def2mtzvp), & + & new_unittest("dft/lanl2", test_dft_lanl2), & + & new_unittest("grad:hf/dz", test_hf_dz_grad), & + ! new_unittest("grad:dft/sv", test_dft_sv_grad), & + & new_unittest("grad:hse3c", test_hse3c_grad), & + ! new_unittest("grad:hf3c", test_hf3c_grad), & + & new_unittest("grad:b973c", test_b973c_grad), & + & new_unittest("without-args", test_without_args) & + & ] + +end subroutine collect_gcp + + +subroutine test_without_args(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(gcp_param) :: param + + call get_structure(mol, "MB16-43", "01") + call get_gcp_param(param, mol) + + call check(error, .not.allocated(param%xv)) + if (allocated(error)) return + call check(error, .not.allocated(param%emiss)) + if (allocated(error)) return + call check(error, .not.allocated(param%slater)) + if (allocated(error)) return + call check(error, .not.allocated(param%zeff)) + if (allocated(error)) return +end subroutine test_without_args + + +subroutine test_hf3c(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(gcp_param) :: param + real(wp) :: energy + + call get_structure(mol, "MB16-43", "02") + call get_gcp_param(param, mol, method="hf3c") + + call check(error, allocated(param%xv)) + if (allocated(error)) return + call check(error, allocated(param%emiss)) + if (allocated(error)) return + call check(error, allocated(param%slater)) + if (allocated(error)) return + call check(error, allocated(param%zeff)) + if (allocated(error)) return + call check(error, param%base) + if (allocated(error)) return + + param%base = .false. + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, 0.0848896136_wp, thr=thr2) + if (allocated(error)) return + + param%base = .true. + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, 0.0292452232_wp, thr=thr2) + if (allocated(error)) return +end subroutine test_hf3c + + +subroutine test_pbeh3c(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(gcp_param) :: param + real(wp) :: energy + + call get_structure(mol, "MB16-43", "03") + call get_gcp_param(param, mol, method="pbeh3c") + + call check(error, allocated(param%xv)) + if (allocated(error)) return + call check(error, allocated(param%emiss)) + if (allocated(error)) return + call check(error, allocated(param%slater)) + if (allocated(error)) return + call check(error, allocated(param%zeff)) + if (allocated(error)) return + call check(error, .not.param%base) + if (allocated(error)) return + call check(error, .not.param%srb) + if (allocated(error)) return + + ! call ascii_gcp_param(6, mol, param) + + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, 0.0206352378_wp, thr=thr2) +end subroutine test_pbeh3c + + +subroutine test_hse3c(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(gcp_param) :: param + real(wp) :: energy + + call get_structure(mol, "MB16-43", "04") + call get_gcp_param(param, mol, method="hse3c") + + call check(error, allocated(param%xv), message="xv") + if (allocated(error)) return + call check(error, allocated(param%emiss), message="emiss") + if (allocated(error)) return + call check(error, allocated(param%slater), message="slater") + if (allocated(error)) return + call check(error, allocated(param%zeff), message="zeff") + if (allocated(error)) return + call check(error, .not.param%base) + if (allocated(error)) return + call check(error, .not.param%srb) + if (allocated(error)) return + + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, 0.0219213037_wp, thr=thr2) +end subroutine test_hse3c + + +subroutine test_b973c(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(gcp_param) :: param + real(wp) :: energy + + call get_structure(mol, "MB16-43", "05") + call get_gcp_param(param, mol, method="b973c") + + call check(error, .not.allocated(param%xv), message="xv") + if (allocated(error)) return + call check(error, .not.allocated(param%emiss), message="emiss") + if (allocated(error)) return + call check(error, .not.allocated(param%slater), message="slater") + if (allocated(error)) return + call check(error, allocated(param%zeff), message="zeff") + if (allocated(error)) return + call check(error, .not.param%base) + if (allocated(error)) return + call check(error, param%srb) + if (allocated(error)) return + + ! call ascii_gcp_param(6, mol, param) + + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, -0.0369737816_wp, thr=thr2) +end subroutine test_b973c + + +subroutine test_r2scan3c(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(gcp_param) :: param + real(wp) :: energy + + call get_structure(mol, "MB16-43", "06") + call get_gcp_param(param, mol, method="r2scan3c") + + call check(error, allocated(param%xv)) + if (allocated(error)) return + call check(error, allocated(param%emiss)) + if (allocated(error)) return + call check(error, allocated(param%slater)) + if (allocated(error)) return + call check(error, allocated(param%zeff)) + if (allocated(error)) return + + ! call ascii_gcp_param(6, mol, param) + + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, 0.0113040952_wp, thr=thr2) + if (allocated(error)) return +end subroutine test_r2scan3c + + +subroutine test_hf_sv(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "07") + call test_generic_energy(error, mol, "hf", "sv", 0.0459359335_wp) +end subroutine test_hf_sv + + +subroutine test_hyb_sv(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "08") + call test_generic_energy(error, mol, "b3lyp", "sv", 0.0817595304_wp) +end subroutine test_hyb_sv + + +subroutine test_gga_sv(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "09") + call test_generic_energy(error, mol, "blyp", "sv", 0.0527050316_wp) +end subroutine test_gga_sv + + +subroutine test_hf_sv_p(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "10") + call test_generic_energy(error, mol, "hf", "sv(p)", 0.0320180286_wp) +end subroutine test_hf_sv_p + + +subroutine test_dft_sv_p(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "11") + call test_generic_energy(error, mol, "dft", "sv(p)", 0.0481858368_wp) +end subroutine test_dft_sv_p + + +subroutine test_dft_svx(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "12") + call test_generic_energy(error, mol, "dft", "svx", 0.0442103186_wp) +end subroutine test_dft_svx + + +subroutine test_hf_svp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "13") + call test_generic_energy(error, mol, "hf", "svp", 0.0389762270_wp) +end subroutine test_hf_svp + + +subroutine test_tpss_svp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "14") + call test_generic_energy(error, mol, "tpss", "svp", 0.0361827169_wp) +end subroutine test_tpss_svp + + +subroutine test_pw6b95_svp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "15") + call test_generic_energy(error, mol, "pw6b95", "svp", 0.0308635906_wp) +end subroutine test_pw6b95_svp + + +subroutine test_hyb_svp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "16") + call test_generic_energy(error, mol, "b3lyp", "svp", 0.0310927236_wp) +end subroutine test_hyb_svp + + +subroutine test_gga_svp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "17") + call test_generic_energy(error, mol, "blyp", "svp", 0.0956240694_wp) +end subroutine test_gga_svp + + +subroutine test_hf_minis(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "18") + call test_generic_energy(error, mol, "hf", "minis", 0.0997731613_wp) +end subroutine test_hf_minis + + +subroutine test_tpss_minis(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "19") + call test_generic_energy(error, mol, "tpss", "minis", 0.0888263651_wp) +end subroutine test_tpss_minis + + +subroutine test_pw6b95_minis(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "20") + call test_generic_energy(error, mol, "pw6b95", "minis", 0.1107806553_wp) +end subroutine test_pw6b95_minis + + +subroutine test_hyb_minis(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "21") + call test_generic_energy(error, mol, "b3lyp", "minis", 0.1019823129_wp) +end subroutine test_hyb_minis + + +subroutine test_gga_minis(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "22") + call test_generic_energy(error, mol, "blyp", "minis", 0.1185302176_wp) +end subroutine test_gga_minis + + +subroutine test_hf_631gd(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "23") + call test_generic_energy(error, mol, "hf", "631gd", 0.0164851915_wp) +end subroutine test_hf_631gd + + +subroutine test_dft_631gd(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "24") + call test_generic_energy(error, mol, "dft", "631gd", 0.0691171349_wp) +end subroutine test_dft_631gd + + +subroutine test_hf_tz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "25") + call test_generic_energy(error, mol, "hf", "tz", 0.0079158812_wp) +end subroutine test_hf_tz + + +subroutine test_hyb_tz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "26") + call test_generic_energy(error, mol, "b3lyp", "tz", 0.0094799906_wp) +end subroutine test_hyb_tz + + +subroutine test_gga_tz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "27") + call test_generic_energy(error, mol, "blyp", "tz", 0.0010732743_wp) +end subroutine test_gga_tz + + +subroutine test_hf_deftzvp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "28") + call test_generic_energy(error, mol, "hf", "deftzvp", 0.0173709952_wp) +end subroutine test_hf_deftzvp + + +subroutine test_dft_deftzvp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "29") + call test_generic_energy(error, mol, "b3lyp", "deftzvp", 0.0162351039_wp) +end subroutine test_dft_deftzvp + + +subroutine test_hf_ccdz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "30") + call test_generic_energy(error, mol, "hf", "ccdz", 0.0531963593_wp) +end subroutine test_hf_ccdz + + +subroutine test_dft_ccdz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "31") + call test_generic_energy(error, mol, "b3lyp", "ccdz", 0.0369893530_wp) +end subroutine test_dft_ccdz + + +subroutine test_hf_accdz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "32") + call test_generic_energy(error, mol, "hf", "accdz", 0.0058924811_wp) +end subroutine test_hf_accdz + + +subroutine test_dft_accdz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "33") + call test_generic_energy(error, mol, "b3lyp", "accdz", 0.0065641768_wp) +end subroutine test_dft_accdz + + +subroutine test_dft_pobtz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "34") + call test_generic_energy(error, mol, "dft", "pobtz", 0.0630483570_wp) +end subroutine test_dft_pobtz + + +subroutine test_hf_minix(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "35") + call test_generic_energy(error, mol, "hf", "minix", 0.0629778826_wp) +end subroutine test_hf_minix + + +subroutine test_dft_minix(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "36") + call test_generic_energy(error, mol, "b3lyp", "minix", 0.1777559799_wp) +end subroutine test_dft_minix + + +subroutine test_hf_2g(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "37") + call test_generic_energy(error, mol, "hf", "2g", 0.3699081462_wp, threshold=thr2*10) +end subroutine test_hf_2g + + +subroutine test_hf_dzp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "38") + call test_generic_energy(error, mol, "hf", "dzp", 0.0237741643_wp) +end subroutine test_hf_dzp + + +subroutine test_dft_dzp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "39") + call test_generic_energy(error, mol, "b3lyp", "dzp", 0.0621712827_wp) +end subroutine test_dft_dzp + + +subroutine test_hf_dz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "40") + call test_generic_energy(error, mol, "hf", "dz", 0.0246769955_wp) +end subroutine test_hf_dz + + +subroutine test_dft_dz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "41") + call test_generic_energy(error, mol, "b3lyp", "dz", 0.0721913136_wp) +end subroutine test_dft_dz + + +subroutine test_b3pbe3c_def2mtzvp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "42") + call test_generic_energy(error, mol, "b3pbe3c", "def2mtzvp", 0.0712990259_wp) +end subroutine test_b3pbe3c_def2mtzvp + + +subroutine test_dft_lanl2(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "43") + call test_generic_energy(error, mol, "b3lyp", "lanl", 0.0471216257_wp) +end subroutine test_dft_lanl2 + + +subroutine test_generic_energy(error, mol, method, basis, reference_energy, threshold) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + !> Molecular structure data + type(structure_type), intent(in) :: mol + + !> Method name + character(len=*), intent(in) :: method + + !> Basis name + character(len=*), intent(in) :: basis + + !> Reference energy + real(wp), intent(in) :: reference_energy + + !> Threshold + real(wp), intent(in), optional :: threshold + + type(gcp_param) :: param + real(wp) :: energy, thr_ + + thr_ = thr2 + if (present(threshold)) thr_ = threshold + + call get_gcp_param(param, mol, method=method, basis=basis) + + call check(error, allocated(param%xv), "Missing number of virtual orbitals") + if (allocated(error)) return + call check(error, allocated(param%emiss), "Missing BSSE parameters") + if (allocated(error)) return + call check(error, allocated(param%slater), "Missing Slater parameters") + if (allocated(error)) return + call check(error, allocated(param%zeff), "Missing effective nuclear charges") + if (allocated(error)) return + + ! call ascii_gcp_param(6, mol, param) + + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, reference_energy, thr=thr_) + if (allocated(error)) then + print *, energy, reference_energy, energy - reference_energy + end if +end subroutine test_generic_energy + +subroutine test_sv(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "01") + call test_generic_basis(error, mol, "sv", 0.3076865293_wp, threshold=thr2*10) +end subroutine test_sv + +subroutine test_sv_p(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "02") + call test_generic_basis(error, mol, "sv(p)", 0.2177767255_wp) +end subroutine test_sv_p + +subroutine test_svx(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "03") + call test_generic_basis(error, mol, "svx", 0.1983990319_wp) +end subroutine test_svx + +subroutine test_svp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "04") + call test_generic_basis(error, mol, "svp", 0.1319325189_wp) +end subroutine test_svp + +subroutine test_minis(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "05") + call test_generic_basis(error, mol, "minis", 4.0137040360_wp, threshold=thr2*100) +end subroutine test_minis + +subroutine test_631gd(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "12") + call test_generic_basis(error, mol, "631gd", 0.0927430940_wp) +end subroutine test_631gd + +subroutine test_tz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "06") + call test_generic_basis(error, mol, "tz", 0.0263158967_wp) +end subroutine test_tz + +subroutine test_deftzvp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "07") + call test_generic_basis(error, mol, "deftzvp", 0.0328360246_wp) +end subroutine test_deftzvp + +subroutine test_ccdz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "08") + call test_generic_basis(error, mol, "ccdz", 0.0688930836_wp) +end subroutine test_ccdz + +subroutine test_accdz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "09") + call test_generic_basis(error, mol, "accdz", 0.0119586199_wp) +end subroutine test_accdz + +subroutine test_pobtz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "10") + call test_generic_basis(error, mol, "pobtz", 0.1094965965_wp) +end subroutine test_pobtz + +subroutine test_minix(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "11") + call test_generic_basis(error, mol, "minix", 2.0696864301_wp, threshold=thr2*10) +end subroutine test_minix + +subroutine test_2g(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "13") + call test_generic_basis(error, mol, "2g", 4.6190048085_wp, threshold=thr2*100) +end subroutine test_2g + +subroutine test_dzp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "14") + call test_generic_basis(error, mol, "dzp", 0.0675215544_wp) +end subroutine test_dzp + +subroutine test_dz(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "15") + call test_generic_basis(error, mol, "dz", 0.0892315292_wp) +end subroutine test_dz + +subroutine test_msvp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "16") + call test_generic_basis(error, mol, "msvp", 0.0927839764_wp) +end subroutine test_msvp + +subroutine test_def2mtzvp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "17") + call test_generic_basis(error, mol, "def2mtzvp", 0.0337359941_wp) +end subroutine test_def2mtzvp + +subroutine test_lanl(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "18") + call test_generic_basis(error, mol, "lanl", 0.1276953745_wp) +end subroutine test_lanl + +subroutine test_def2mtzvpp(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "19") + call test_generic_basis(error, mol, "def2mtzvpp", 0.0161854167_wp) +end subroutine test_def2mtzvpp + +subroutine test_vmb(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + call get_structure(mol, "MB16-43", "20") + call test_generic_basis(error, mol, "vmb", 1.6197120080_wp) +end subroutine test_vmb + +subroutine test_generic_basis(error, mol, basis, reference_energy, threshold) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + !> Molecular structure data + type(structure_type), intent(in) :: mol + + !> Basis name + character(len=*), intent(in) :: basis + + !> Reference energy + real(wp), intent(in) :: reference_energy + + !> Threshold + real(wp), intent(in), optional :: threshold + + type(gcp_param) :: param + real(wp) :: energy, thr_ + + thr_ = thr2 + if (present(threshold)) thr_ = threshold + + call get_gcp_param(param, mol, basis=basis, eta=1.0_wp) + param%sigma = 1.0_wp + param%alpha = 1.0_wp + param%beta = 1.0_wp + + call check(error, allocated(param%xv), "Missing number of virtual orbitals") + if (allocated(error)) return + call check(error, allocated(param%emiss), "Missing BSSE parameters") + if (allocated(error)) return + call check(error, allocated(param%slater), "Missing Slater parameters") + if (allocated(error)) return + call check(error, allocated(param%zeff), "Missing effective nuclear charges") + if (allocated(error)) return + + ! call ascii_gcp_param(6, mol, param) + + energy = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy) + + call check(error, energy, reference_energy, thr=thr_) + if (allocated(error)) then + print *, energy, reference_energy, energy - reference_energy + end if +end subroutine test_generic_basis + + +subroutine test_hf_dz_grad(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "03") + call test_numgrad(error, mol, "hf", "dz") + +end subroutine test_hf_dz_grad + + +subroutine test_dft_sv_grad(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "X23", "CO2") + call test_numgrad(error, mol, "dft", "sv") + +end subroutine test_dft_sv_grad + + +subroutine test_hse3c_grad(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "02") + call test_numgrad(error, mol, "hse3c") + +end subroutine test_hse3c_grad + + +subroutine test_hf3c_grad(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "X23", "ammonia") + call test_numgrad(error, mol, "hf3c") + +end subroutine test_hf3c_grad + + +subroutine test_b973c_grad(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "03") + call test_numgrad(error, mol, "b973c") + +end subroutine test_b973c_grad + + +subroutine test_numgrad(error, mol, method, basis) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Method name + character(len=*), intent(in) :: method + + !> Method name + character(len=*), intent(in), optional :: basis + + type(gcp_param) :: param + integer :: iat, ic, mat + real(wp) :: energy, er, el + real(wp), allocatable :: gradient(:, :), numgrad(:, :) + real(wp) :: sigma(3, 3) + logical :: pbc + real(wp), parameter :: step = 1.0e-6_wp + real(wp), parameter :: thrg = sqrt(epsilon(1.0_wp)) + + pbc = any(mol%periodic) + + allocate(gradient(3, mol%nat), numgrad(3, mol%nat)) + + call get_gcp_param(param, mol, method=method, basis=basis) + + if (pbc) then + mat = min(mol%nat, 3) + else + mat = mol%nat + end if + do iat = 1, mat + do ic = 1, 3 + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + er = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), er) + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + el = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), el) + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + numgrad(ic, iat) = 0.5_wp*(er - el)/step + end do + end do + + gradient(:, :) = 0.0_wp + call get_geometric_counterpoise(mol, param, realspace_cutoff(), energy, gradient, sigma) + + if (any(abs(gradient(:, :mat)-numgrad(:, :mat)) > thrg)) then + call test_failed(error, "Numerical and analytical gradient do not match") + do iat = 1, mat + print'(3es14.5,3x,a)', gradient(:, iat)-numgrad(:, iat), mol%sym(mol%id(iat)) + end do + end if + +end subroutine test_numgrad + + +end module test_gcp \ No newline at end of file diff --git a/test/validation/21-ad3.coord b/test/validation/21-ad3.coord new file mode 100644 index 00000000..7c0ff6df --- /dev/null +++ b/test/validation/21-ad3.coord @@ -0,0 +1,24 @@ +$coord + 0.79457714965951 3.55602905386715 0.00000000000000 c + -0.79457714965951 3.77432460885587 2.39660993619036 c + -0.79457714965951 3.77432460885587 -2.39660993619036 c + 2.25063800006750 5.02905587026014 0.00000000000000 h + 1.79951300232100 1.74494525518548 0.00000000000000 h + 0.37035808014807 3.61430066526315 4.09670630847074 h + -2.22853764673775 2.28531070488959 2.46581361633477 h + -1.77392043411708 5.59484192304179 2.46581361633477 h + 0.37035808014807 3.61430066526315 -4.09670630847074 h + -2.22853764673775 2.28531070488959 -2.46581361633477 h + -1.77392043411708 5.59484192304179 -2.46581361633477 h + -0.79457714965951 -3.55602905386715 0.00000000000000 c + 0.79457714965951 -3.77432460885587 -2.39660993619036 c + 0.79457714965951 -3.77432460885587 2.39660993619036 c + -1.79951300232100 -1.74494525518548 0.00000000000000 h + -2.25063800006750 -5.02905587026014 0.00000000000000 h + -0.37035808014807 -3.61430066526315 -4.09670630847074 h + 1.77392043411708 -5.59484192304179 -2.46581361633477 h + 2.22853764673775 -2.28531070488959 -2.46581361633477 h + -0.37035808014807 -3.61430066526315 4.09670630847074 h + 1.77392043411708 -5.59484192304179 2.46581361633477 h + 2.22853764673775 -2.28531070488959 2.46581361633477 h +$end \ No newline at end of file diff --git a/test/validation/21-gcp-pbeh3c.json b/test/validation/21-gcp-pbeh3c.json new file mode 100644 index 00000000..79ddd7f9 --- /dev/null +++ b/test/validation/21-gcp-pbeh3c.json @@ -0,0 +1,83 @@ +{ + "version": "1.2.1", + "energy": 2.4408581444610614E-02, + "virial": [ + -1.6574359800216241E-02, + 2.0119627561238601E-03, + 1.6534083130403943E-18, + 2.0119627561238618E-03, + -2.2463098851840998E-03, + -1.6940658945086007E-19, + 1.1113072267976420E-18, + -4.4045713257223618E-19, + -3.7753665063192894E-02 + ], + "gradient": [ + -2.8026033066000148E-03, + 2.4390779524848948E-04, + -1.0842021724855044E-19, + 1.1950451267053234E-03, + -2.3762562720102155E-04, + -2.0663035801345331E-03, + 1.1950451267053234E-03, + -2.3762562720102155E-04, + 2.0663035801345335E-03, + -6.8360546166405974E-04, + 1.0124060154628643E-06, + 0.0000000000000000E+00, + -7.2054440646042786E-04, + -1.4820156713158696E-05, + 0.0000000000000000E+00, + 4.5251185867576900E-04, + -6.3822465818580702E-05, + -2.4669490466976086E-04, + 1.9648267644481823E-04, + -1.9837755776525693E-04, + -7.0986894457992632E-04, + 2.2138861006725468E-04, + 1.0653376468129232E-04, + -6.9966824801623285E-04, + 4.5251185867576895E-04, + -6.3822465818580702E-05, + 2.4669490466976092E-04, + 1.9648267644481828E-04, + -1.9837755776525693E-04, + 7.0986894457992643E-04, + 2.2138861006725473E-04, + 1.0653376468129233E-04, + 6.9966824801623285E-04, + 2.8026033066000143E-03, + -2.4390779524848953E-04, + 1.0842021724855044E-19, + -1.1950451267053236E-03, + 2.3762562720102152E-04, + 2.0663035801345331E-03, + -1.1950451267053236E-03, + 2.3762562720102163E-04, + -2.0663035801345331E-03, + 7.2054440646042786E-04, + 1.4820156713158750E-05, + 0.0000000000000000E+00, + 6.8360546166405984E-04, + -1.0124060154628376E-06, + 0.0000000000000000E+00, + -4.5251185867576900E-04, + 6.3822465818580702E-05, + 2.4669490466976081E-04, + -2.2138861006725468E-04, + -1.0653376468129230E-04, + 6.9966824801623274E-04, + -1.9648267644481828E-04, + 1.9837755776525696E-04, + 7.0986894457992643E-04, + -4.5251185867576900E-04, + 6.3822465818580702E-05, + -2.4669490466976086E-04, + -2.2138861006725468E-04, + -1.0653376468129232E-04, + -6.9966824801623263E-04, + -1.9648267644481834E-04, + 1.9837755776525696E-04, + -7.0986894457992632E-04 + ] +} diff --git a/test/validation/21-gcp-pbeh3c.resp b/test/validation/21-gcp-pbeh3c.resp new file mode 100644 index 00000000..3cfae49e --- /dev/null +++ b/test/validation/21-gcp-pbeh3c.resp @@ -0,0 +1,6 @@ +gcp +$ORIGIN/21-ad3.coord +--level +pbeh3c +--grad +--nocpc \ No newline at end of file diff --git a/test/validation/22-benI-acetone.xyz b/test/validation/22-benI-acetone.xyz new file mode 100644 index 00000000..1227ac8f --- /dev/null +++ b/test/validation/22-benI-acetone.xyz @@ -0,0 +1,24 @@ +22 + +C -0.43601987190704 -0.00000150300001 1.03540729252862 +C 0.58976220857948 -0.00001321400006 1.97850007065348 +C 1.92224732040744 -0.00000620100003 1.57169559487422 +C 2.22826604374589 -0.00002197200010 0.21256912192973 +C 1.21108714329700 -0.00000640500003 -0.74015358323725 +H -1.46924203642611 -0.00000195000001 1.35222461891430 +H 0.34310109250064 -0.00000842500004 3.03162969825961 +H 2.71469423987341 0.00001150200005 2.30703292809040 +H 3.26026141025959 -0.00002655400011 -0.11155152048790 +H 1.45070781434504 -0.00001010800004 -1.79398552984646 +C -0.12040282452661 0.00002212200010 -0.32357594141524 +I -1.65822411125267 0.00022120300097 -1.73801997760168 +O -3.94993994027608 -0.00014345000063 -3.81589692468980 +C -3.80061832262298 0.00377150001650 -5.02581161398168 +C -2.41374687855714 0.01392866206092 -5.62884087361918 +C -5.00409199488668 -0.00151565300663 -5.94221514498980 +H -1.88103810222720 0.89133045189846 -5.26547224702990 +H -1.87007825717927 -0.85848332775480 -5.26971232504844 +H -2.42067773158745 0.01652053407226 -6.71412316936594 +H -5.61626987956420 0.87012827680573 -5.71741149300657 +H -5.60537142251653 -0.88176795885664 -5.72158319002481 +H -4.74004724773182 0.00263918901154 -6.99488255559392 diff --git a/test/validation/22-gcp-b973c.json b/test/validation/22-gcp-b973c.json new file mode 100644 index 00000000..876d4802 --- /dev/null +++ b/test/validation/22-gcp-b973c.json @@ -0,0 +1,13 @@ +{ + "version": "1.2.1", + "energy": -1.1487620744730206E-01, + "damping parameters": { + "damping": "rational", + "s6": 1.0000000000000000E+00, + "s8": 1.5000000000000000E+00, + "s9": 0.0000000000000000E+00, + "a1": 3.7000000000000000E-01, + "a2": 4.0999999999999996E+00, + "alp": 1.4000000000000000E+01 + } +} diff --git a/test/validation/22-gcp-b973c.resp b/test/validation/22-gcp-b973c.resp new file mode 100644 index 00000000..7ccf510a --- /dev/null +++ b/test/validation/22-gcp-b973c.resp @@ -0,0 +1,6 @@ +run +--bj +b973c +--gcp +--noedisp +$ORIGIN/22-benI-acetone.xyz \ No newline at end of file diff --git a/test/validation/23-444.sdf b/test/validation/23-444.sdf new file mode 100644 index 00000000..8f17520c --- /dev/null +++ b/test/validation/23-444.sdf @@ -0,0 +1,43 @@ + + 11302419103D + + 37 0 0 0 0 999 V2000 + 2.5903 -2.7264 1.1396 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.9003 -1.7974 0.4296 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.8263 -1.8644 -0.9034 O 0 0 0 0 0 0 0 0 0 0 0 0 + 3.1053 -0.9954 -1.2824 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.4013 -0.4684 0.9836 C 0 0 0 0 0 0 0 0 0 0 0 0 + 4.4123 -0.2904 0.6116 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.4233 -0.5444 2.0656 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.5623 0.6566 0.6136 N 0 0 0 0 0 0 0 0 0 0 0 0 + 1.8913 1.0566 1.2686 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.5763 1.1686 -0.6264 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.2463 0.6926 -1.5524 O 0 0 0 0 0 0 0 0 0 0 0 0 + 1.7003 2.3976 -0.8174 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.9023 2.7876 -1.8094 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.9593 3.1466 -0.0674 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2873 2.0946 -0.7124 N 0 0 0 0 0 0 0 0 0 0 0 0 + -0.3337 2.1456 -1.5084 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.3257 1.9366 0.4716 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2703 1.9106 1.5546 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8457 1.8206 0.3966 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.2987 1.8046 -0.9934 N 0 0 0 0 0 0 0 0 0 0 0 0 + -3.2037 2.2476 -1.0744 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.4147 0.8416 -1.2984 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.2107 2.7366 0.8656 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.3197 0.6246 1.2456 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8237 0.6816 2.2136 H 0 0 0 0 0 0 0 0 0 0 0 0 + -3.3947 0.7276 1.4056 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.0327 -0.6824 0.5636 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.0087 -1.2834 -0.2404 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.9927 -0.8344 -0.3024 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.7397 -2.4624 -0.9334 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.5087 -2.9154 -1.5444 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.4857 -3.0574 -0.8284 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.2737 -3.9754 -1.3594 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.5057 -2.4784 -0.0224 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4623 -2.9484 0.0826 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7787 -1.2964 0.6626 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.0207 -0.8594 1.3016 H 0 0 0 0 0 0 0 0 0 0 0 0 +M END +$$$$ diff --git a/test/validation/23-gcp-hf3c.json b/test/validation/23-gcp-hf3c.json new file mode 100644 index 00000000..b50fa6f1 --- /dev/null +++ b/test/validation/23-gcp-hf3c.json @@ -0,0 +1,4 @@ +{ + "version": "1.2.1", + "energy": -2.0199949882394205E-01 +} diff --git a/test/validation/23-gcp-hf3c.resp b/test/validation/23-gcp-hf3c.resp new file mode 100644 index 00000000..74988592 --- /dev/null +++ b/test/validation/23-gcp-hf3c.resp @@ -0,0 +1,4 @@ +gcp +--level +hf3c +$ORIGIN/23-444.sdf \ No newline at end of file diff --git a/test/validation/24-gcp-hf-deftzvp.json b/test/validation/24-gcp-hf-deftzvp.json new file mode 100644 index 00000000..62eedc18 --- /dev/null +++ b/test/validation/24-gcp-hf-deftzvp.json @@ -0,0 +1,13 @@ +{ + "version": "1.2.1", + "energy": 7.1769896849684205E-03, + "damping parameters": { + "damping": "zero", + "s6": 1.0000000000000000E+00, + "s8": 1.7460000000000000E+00, + "s9": 1.0000000000000000E+00, + "rs6": 1.1579999999999999E+00, + "rs8": 1.0000000000000000E+00, + "alp": 1.4000000000000000E+01 + } +} diff --git a/test/validation/24-gcp-hf-deftzvp.resp b/test/validation/24-gcp-hf-deftzvp.resp new file mode 100644 index 00000000..086f4358 --- /dev/null +++ b/test/validation/24-gcp-hf-deftzvp.resp @@ -0,0 +1,7 @@ +run +--zero +hf +--gcp +deftzvp +--atm +$ORIGIN/24-h2o_6.qchem \ No newline at end of file diff --git a/test/validation/24-h2o_6.qchem b/test/validation/24-h2o_6.qchem new file mode 100644 index 00000000..c26dc1e2 --- /dev/null +++ b/test/validation/24-h2o_6.qchem @@ -0,0 +1,21 @@ +$molecule + 0 1 +o 0.11448392281773 1.51123327008092 0.84351297623233 +h 0.95880573623723 1.47155799179987 0.32950022850208 +h 0.26129103645727 2.12410912352789 1.56850115482200 +o -0.02817040100173 -1.43708228218242 0.98511825009724 +h -0.03831273056231 -0.48623197337830 1.17189380200513 +h -0.86424990477834 -1.59839377519356 0.50558843087301 +o -2.35379807024322 1.32801810801878 -0.46182075721370 +h -1.48714230779134 1.50194038048355 -0.04531091525911 +h -2.34725343052677 1.79416500782283 -1.30125669084489 +o -2.46767106530964 -1.45760274129337 -0.41854609908859 +h -2.53975419218662 -0.48503716343007 -0.49641687571481 +h -3.26029002096919 -1.74564200881398 0.04080615100984 +o 2.42644955265113 1.14816004581118 -0.49586368573878 +h 2.48383495016489 0.16184556854355 -0.50914940516318 +h 2.52619219832975 1.43384884343363 -1.40718630625545 +o 2.27629136915678 -1.54900494733334 -0.42626722875407 +h 1.42256374614478 -1.63495392360958 0.06497690996263 +h 2.91672961140960 -2.08092952428755 0.05192006052832 +$end diff --git a/test/validation/25-g3.gen b/test/validation/25-g3.gen new file mode 100644 index 00000000..ebeb5a93 --- /dev/null +++ b/test/validation/25-g3.gen @@ -0,0 +1,26 @@ +24 C + C O H + 1 1 -1.12885064915344E-01 -1.41863844617349E+00 -7.80362498298278E-01 + 2 2 4.09592407421189E-01 -2.36548698407175E-01 -1.36364838541136E+00 + 3 1 6.19799977621229E-01 8.89905751402348E-01 -5.02538213519736E-01 + 4 1 -6.10116041805569E-01 1.20757778290933E+00 3.67780827198579E-01 + 5 1 -1.43900587170863E+00 -1.16243455996366E+00 -4.62672873964163E-02 + 6 1 -1.24760318997206E+00 -5.08890838053665E-02 1.00951572052529E+00 + 7 2 -4.86797482060935E-01 -4.99194756761102E-01 2.12454981438015E+00 + 8 2 -1.66043778109922E+00 1.79383387185583E+00 -4.22680367436187E-01 + 9 1 1.93115068908887E+00 7.74794098890431E-01 3.28828322366809E-01 + 10 2 2.88225415460348E+00 -1.16345429245288E-01 -2.56160344296218E-01 + 11 3 1.29113665504839E-01 -1.18582686926397E+00 1.80296413321587E+00 + 12 3 -1.33590286621933E+00 2.59430229650854E+00 -8.64241078868696E-01 + 13 3 2.97633759330056E+00 8.92087438907900E-02 -1.20071451738229E+00 + 14 2 7.93330183785820E-01 -1.98847595858884E+00 1.60863527434552E-01 + 15 2 -2.47074754876863E+00 -8.63512585694265E-01 -9.70219081919390E-01 + 16 3 -2.23180271908296E+00 2.45770664766718E-01 1.38586283862572E+00 + 17 3 -3.18410548179756E-01 1.89503777569582E+00 1.17421385796704E+00 + 18 3 7.51511780065613E-01 1.71194211464097E+00 -1.21902350171035E+00 + 19 3 2.36346506806486E+00 1.78005871805928E+00 4.37353614115734E-01 + 20 3 1.70348772650670E+00 -1.75009521752220E+00 -1.11972385813754E-01 + 21 3 -2.39388948446623E+00 7.67028777273814E-02 -1.21635880287680E+00 + 22 3 -2.65720769643726E-01 -2.08828234086970E+00 -1.63230513539516E+00 + 23 3 -1.73289577975631E+00 -2.08526093479100E+00 4.64247803351864E-01 + 24 3 1.74617190171553E+00 3.86370184738622E-01 1.33031114114306E+00 diff --git a/test/validation/25-gcp-b3lyp-631gd.json b/test/validation/25-gcp-b3lyp-631gd.json new file mode 100644 index 00000000..53a95b81 --- /dev/null +++ b/test/validation/25-gcp-b3lyp-631gd.json @@ -0,0 +1,89 @@ +{ + "version": "1.2.1", + "energy": 9.3584755369634995E-02, + "virial": [ + -8.1810660539116148E-02, + -9.6093525430042008E-05, + -1.3672586053845414E-04, + -9.6093525430033822E-05, + -8.1086316048136181E-02, + 2.6884889728128747E-03, + -1.3672586053845327E-04, + 2.6884889728128769E-03, + -6.9896275544401601E-02 + ], + "gradient": [ + 7.2728850731336030E-05, + 1.4180915671166015E-04, + 1.6261536689883918E-04, + -3.2812414702953552E-04, + -9.6291499805133672E-05, + 2.8585524050143050E-03, + 1.1530950793400656E-06, + 4.4394105707793303E-05, + 2.2859757014506683E-05, + 1.5704086214981611E-04, + -6.9760254879827824E-05, + 4.8776502098259354E-05, + 2.7803749513078053E-04, + 8.0445644544793331E-05, + 2.1536273375148256E-04, + -9.1001875781775095E-06, + -3.1382705708857148E-05, + -2.8908539135137200E-04, + 9.0521100022711862E-04, + -7.1275780796299670E-04, + -3.7621158749063439E-03, + 2.4600714427851058E-03, + -6.7297758581742154E-04, + 3.3247966137161005E-06, + -2.9861834386274064E-04, + -1.5141183384514436E-04, + 1.7158713453133080E-04, + -2.8018803569370364E-03, + 1.0769449589210006E-03, + -1.1034076359568937E-03, + -1.4942950694016729E-03, + 1.8988535063400643E-03, + -1.1230645409583207E-03, + -2.0482168026468152E-04, + -2.7953944201808637E-03, + 1.4740456014026243E-03, + -1.3710125284322350E-03, + -3.7898366040712731E-04, + 2.8764723118768682E-03, + 1.8885118020500004E-04, + 2.9908946774985300E-03, + -1.0400408033423422E-03, + 2.6727382615169670E-03, + 2.1580838440340890E-03, + 1.2568082812288058E-03, + 2.5547141679660656E-03, + -6.4673074021500733E-04, + -1.4790393947981288E-03, + -5.8467010454819131E-04, + -2.0336642369187391E-03, + -2.0365272866933945E-03, + -6.5013599295336030E-04, + -2.2455940125889404E-03, + 1.8185293287786258E-03, + -1.1638022762662889E-03, + -2.5594896062603911E-03, + -5.6359284689985419E-04, + -2.6910174215525641E-03, + 1.1131745039147612E-03, + 2.9384988490817972E-04, + 1.1343537456426651E-03, + -1.8440875951904231E-03, + 1.6511634393546679E-03, + 3.7969649933155322E-04, + 2.0849072919157260E-03, + 2.3577195538768436E-03, + 1.2255060371075487E-03, + 2.5563382240707719E-03, + -1.1503593321425839E-03, + -4.3262452904681314E-04, + 9.2680046121683160E-05, + -2.6644339902998199E-03 + ] +} diff --git a/test/validation/25-gcp-b3lyp-631gd.resp b/test/validation/25-gcp-b3lyp-631gd.resp new file mode 100644 index 00000000..8226a1c8 --- /dev/null +++ b/test/validation/25-gcp-b3lyp-631gd.resp @@ -0,0 +1,8 @@ +gcp +$ORIGIN/25-g3.gen +-i +gen +--nocpc +--level +b3lyp/631gd +--grad \ No newline at end of file diff --git a/test/validation/meson.build b/test/validation/meson.build index 0bc3775f..b3825d25 100644 --- a/test/validation/meson.build +++ b/test/validation/meson.build @@ -36,6 +36,11 @@ tests = [ '18-db-d3-bjm', '19-db-d3-zerom', '20-db-d3-op', + '21-gcp-pbeh3c', + '22-gcp-b973c', + '23-gcp-hf3c', + '24-gcp-hf-deftzvp', + '25-gcp-b3lyp-631gd', ] foreach t : tests diff --git a/test/validation/tester.py b/test/validation/tester.py index 73f2b17e..762b8118 100644 --- a/test/validation/tester.py +++ b/test/validation/tester.py @@ -23,7 +23,7 @@ args = [arg.replace('$ORIGIN', wdir) for arg in fd.read().strip().split("\n")] stat = subprocess.call( - [prog, "--json", os.path.basename(outp)] + args, + [prog, *args, "--json", os.path.basename(outp)], shell=False, stdin=None, stderr=subprocess.STDOUT,