From d12df6f3321fa3ff5020cc93811bc6705590aa22 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Mon, 18 May 2020 23:00:57 +0200 Subject: [PATCH 1/4] C-API version 6.3 - direct control over the Fortran objects - (almost) thread-safe implementation - separated APIs for IO, structure and calulations --- .github/workflows/fortran-build.yml | 2 +- .travis.yml | 2 +- CMakeLists.txt | 3 - TESTSUITE/CMakeLists.txt | 4 +- TESTSUITE/c_api_example.c | 60 +- TESTSUITE/gfn0.f90 | 36 +- TESTSUITE/gfn1.f90 | 130 ++-- TESTSUITE/gfn2.f90 | 119 +-- TESTSUITE/meson.build | 2 +- TESTSUITE/peeq.f90 | 40 +- TESTSUITE/symmetry.f90 | 2 +- cmake/CMakeLists.txt | 3 - include/xtb.h | 496 ++++++------- meson.build | 60 +- meson_options.txt | 2 +- python/xtb/__init__.py | 5 +- python/xtb/calculators.py | 265 ------- python/xtb/interface.py | 841 +++++++++++---------- python/xtb/solvation.py | 133 ---- python/xtb/test/test_calculators.py | 117 --- python/xtb/test/test_interface.py | 95 +-- python/xtb/test/test_solvation.py | 125 ---- src/CMakeLists.txt | 4 - src/api/CMakeLists.txt | 6 +- src/api/calculator.f90 | 511 +++++++++++++ src/api/environment.f90 | 205 +++++ src/api/interface.f90 | 1068 +++------------------------ src/api/meson.build | 6 +- src/api/molecule.f90 | 159 ++++ src/api/preload.f90 | 196 ----- src/api/results.f90 | 270 +++++++ src/api/structs.f90 | 328 -------- src/api/utils.f90 | 307 +------- src/calculator.f90 | 83 --- src/disp/dftd4.f90 | 4 - src/filetools.F90 | 213 +----- src/gfn0_calculator.f90 | 173 ----- src/gfn1_calculator.f90 | 184 ----- src/gfn2_calculator.f90 | 180 ----- src/main/property.f90 | 1 + src/main/setup.f90 | 156 +++- src/mctc/mctc_global.f90 | 3 +- src/mctc/mctc_init.f90 | 1 + src/meson.build | 4 - src/printout.f90 | 5 +- src/prog/main.F90 | 6 +- src/prog/topology.f90 | 5 +- src/scf_module.f90 | 20 +- src/setparam.f90 | 7 +- src/type/CMakeLists.txt | 1 + src/type/environment.f90 | 6 + src/type/iohandler.f90 | 470 ++++++++++++ src/type/meson.build | 1 + src/xtb/calculator.f90 | 28 +- 54 files changed, 2938 insertions(+), 4215 deletions(-) delete mode 100755 python/xtb/calculators.py delete mode 100644 python/xtb/solvation.py delete mode 100644 python/xtb/test/test_calculators.py delete mode 100644 python/xtb/test/test_solvation.py create mode 100644 src/api/calculator.f90 create mode 100644 src/api/environment.f90 create mode 100644 src/api/molecule.f90 delete mode 100644 src/api/preload.f90 create mode 100644 src/api/results.f90 delete mode 100644 src/api/structs.f90 delete mode 100644 src/calculator.f90 delete mode 100644 src/gfn0_calculator.f90 delete mode 100644 src/gfn1_calculator.f90 delete mode 100644 src/gfn2_calculator.f90 create mode 100644 src/type/iohandler.f90 diff --git a/.github/workflows/fortran-build.yml b/.github/workflows/fortran-build.yml index 1730b0bc5..072022c0e 100644 --- a/.github/workflows/fortran-build.yml +++ b/.github/workflows/fortran-build.yml @@ -11,5 +11,5 @@ jobs: python-version: '3.x' - run: brew install gcc@8 ninja - run: pip install meson==0.53.2 - - run: FC=gfortran-8 CC=gcc-8 meson setup build_gcc --buildtype release -Dla_backend=netlib --warnlevel 0 + - run: FC=gfortran-8 CC=gcc-8 meson setup build_gcc --buildtype release -Dla_backend=netlib --warnlevel 0 -Dpython=false - run: OMP_NUM_THREADS=2 ninja -C build_gcc test diff --git a/.travis.yml b/.travis.yml index aa3f1dc2b..9c7d5b750 100644 --- a/.travis.yml +++ b/.travis.yml @@ -33,4 +33,4 @@ jobs: - OMP_NUM_THREADS=2 make -j all test install install: - - pip3 install meson==0.53.2 ninja ase pytest + - pip3 install meson==0.53.2 ninja pytest numpy diff --git a/CMakeLists.txt b/CMakeLists.txt index c0d382ee7..53e0f22f6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -59,9 +59,6 @@ target_include_directories(lib-xtb-static $/include/xtb> ) -message("BLAS_LIBRARIES ${BLAS_LIBRARIES}") -message("LAPACK_LIBRARIES ${LAPACK_LIBRARIES}") - # Shared Library add_library(lib-xtb-shared SHARED $ ) target_link_libraries(lib-xtb-shared diff --git a/TESTSUITE/CMakeLists.txt b/TESTSUITE/CMakeLists.txt index 50ec91a9b..5e5acf496 100644 --- a/TESTSUITE/CMakeLists.txt +++ b/TESTSUITE/CMakeLists.txt @@ -132,7 +132,7 @@ add_test("GFN2-xTB_SCC" ${XTB-TESTER} gfn2 scc) add_test("GFN2-xTB_API" ${XTB-TESTER} gfn2 api) add_test("GFN2-xTB_API_GBSA" ${XTB-TESTER} gfn2 gbsa) add_test("GFN2-xTB_API_GBSA_salt" ${XTB-TESTER} gfn2 salt) -add_test("GFN2-xTB_API_PCEM" ${XTB-TESTER} gfn2 pcem) +#add_test("GFN2-xTB_API_PCEM" ${XTB-TESTER} gfn2 pcem) add_test("GFN1-xTB_SCC" ${XTB-TESTER} gfn1 scc) add_test("GFN1-xTB_API" ${XTB-TESTER} gfn1 api) add_test("GFN1-xTB_API_XB" ${XTB-TESTER} gfn1 xb) @@ -162,7 +162,7 @@ set_tests_properties(EXE_Singlepoint "GFN2-xTB_API" "GFN2-xTB_API_GBSA" "GFN2-xTB_API_GBSA_salt" - "GFN2-xTB_API_PCEM" + #"GFN2-xTB_API_PCEM" "GFN1-xTB_SCC" "GFN1-xTB_API" "GFN1-xTB_API_XB" diff --git a/TESTSUITE/c_api_example.c b/TESTSUITE/c_api_example.c index bcb141301..f42d924b9 100644 --- a/TESTSUITE/c_api_example.c +++ b/TESTSUITE/c_api_example.c @@ -11,6 +11,7 @@ main (int argc, char **argv) int const natoms = 7; int const attyp[7] = {6,6,6,1,1,1,1}; double const charge = 0.0; + int const uhf = 0; double const coord[3*7] = {0.00000000000000, 0.00000000000000,-1.79755622305860, 0.00000000000000, 0.00000000000000, 0.95338756106749, @@ -20,24 +21,65 @@ main (int argc, char **argv) 1.92825631079613, 0.00000000000000,-2.53624948351102, 0.00000000000000, 0.00000000000000, 5.23010455462158}; - SCC_options const opt = (SCC_options){2, 0, 1.0, 300.0, true, false, true, 30, "none"}; - + xtb_TEnvironment env; + xtb_TMolecule mol; + xtb_TCalculator calc; + xtb_TResults res; double energy; double dipole[3]; double q[natoms]; - double qp[6*natoms]; double wbo[natoms*natoms]; - int stat = GFN2_calculation(&natoms, attyp, &charge, NULL, coord, &opt, "-", - &energy, NULL, dipole, q, NULL, qp, wbo); + env = xtb_newEnvironment(); + calc = xtb_newCalculator(); + res = xtb_newResults(); + mol = xtb_newMolecule(env, &natoms, attyp, coord, NULL, NULL, NULL, NULL); + if (xtb_checkEnvironment(env)) { + xtb_showEnvironment(env, NULL); + return 1; + } + + xtb_setVerbosity(env, XTB_VERBOSITY_FULL); + if (xtb_checkEnvironment(env)) { + xtb_showEnvironment(env, NULL); + return 2; + } + + xtb_loadGFN2xTB(env, mol, calc, NULL); + if (xtb_checkEnvironment(env)) { + xtb_showEnvironment(env, NULL); + return 3; + } + + xtb_singlepoint(env, mol, calc, res); + if (xtb_checkEnvironment(env)) { + xtb_showEnvironment(env, NULL); + return 4; + } + + xtb_getEnergy(env, res, &energy); + xtb_getCharges(env, res, q); + xtb_getDipole(env, res, dipole); + xtb_getBondOrders(env, res, wbo); + if (xtb_checkEnvironment(env)) { + xtb_showEnvironment(env, NULL); + return 5; + } - assert(stat == 0); assert(fabs(energy + 8.3824793849585) < 1.0e-9); assert(fabs(q[5] - 0.05184019996829) < 1.0e-8); - assert(fabs(dipole[2] + 0.29832384492435) < 1.0e-6); - assert(fabs(qp[14] - 0.56386138525354) < 1.0e-6); + assert(fabs(dipole[2] + 0.298279305689518) < 1.0e-6); assert(fabs(wbo[9] - 2.89823984265213) < 1.0e-8); + xtb_delResults(&res); + xtb_delCalculator(&calc); + xtb_delMolecule(&mol); + xtb_delEnvironment(&env); + + assert(!res); + assert(!calc); + assert(!mol); + assert(!env); + return 0; } - diff --git a/TESTSUITE/gfn0.f90 b/TESTSUITE/gfn0.f90 index 7aa7bf11a..b22b834b7 100644 --- a/TESTSUITE/gfn0.f90 +++ b/TESTSUITE/gfn0.f90 @@ -130,12 +130,15 @@ subroutine test_gfn0_api use xtb_type_options use xtb_type_molecule + use xtb_type_wavefunction + use xtb_type_data use xtb_type_param use xtb_type_environment use xtb_pbc_tools - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -155,10 +158,13 @@ subroutine test_gfn0_api type(TMolecule) :: mol type(TEnvironment) :: env + type(TWavefunction) :: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc real(wp) :: energy real(wp) :: hl_gap - real(wp) :: dum(3,3) + real(wp) :: sigma(3,3) real(wp),allocatable :: gradient(:,:) ! setup the environment variables @@ -170,8 +176,12 @@ subroutine test_gfn0_api energy = 0.0_wp gradient = 0.0_wp - call gfn0_calculation & - (stdout,env,opt,mol,hl_gap,energy,gradient,dum,dum) + call newXTBCalculator(env, mol, calc, method=0) + call env%checkpoint("failed setup") + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 5.5384029314207_wp,thr) call assert_close(energy,-8.6908532561691_wp,thr) @@ -193,12 +203,15 @@ subroutine test_gfn0_api_srb use xtb_type_options use xtb_type_molecule + use xtb_type_data + use xtb_type_wavefunction use xtb_type_param use xtb_type_environment use xtb_pbc_tools - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -237,10 +250,13 @@ subroutine test_gfn0_api_srb type(TMolecule) :: mol type(TEnvironment) :: env + type(TWavefunction) :: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc real(wp) :: energy real(wp) :: hl_gap - real(wp) :: dum(3,3) + real(wp) :: sigma(3,3) real(wp),allocatable :: gradient(:,:) ! setup the environment variables @@ -252,8 +268,12 @@ subroutine test_gfn0_api_srb energy = 0.0_wp gradient = 0.0_wp - call gfn0_calculation & - (stdout,env,opt,mol,hl_gap,energy,gradient,dum,dum) + call newXTBCalculator(env, mol, calc, method=0) + call env%checkpoint("failed setup") + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 3.1192454818777_wp,thr) call assert_close(energy,-40.908850360158_wp,thr) diff --git a/TESTSUITE/gfn1.f90 b/TESTSUITE/gfn1.f90 index 167ca8827..5dd316112 100644 --- a/TESTSUITE/gfn1.f90 +++ b/TESTSUITE/gfn1.f90 @@ -123,11 +123,12 @@ subroutine test_gfn1_api use xtb_type_options use xtb_type_molecule use xtb_type_param - use xtb_type_pcem + use xtb_type_data use xtb_type_wavefunction use xtb_type_environment - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -147,10 +148,11 @@ subroutine test_gfn1_api type(TMolecule) :: mol type(TEnvironment) :: env - type(tb_pcem) :: pcem type(TWavefunction):: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3, 3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -163,8 +165,11 @@ subroutine test_gfn1_api energy = 0.0_wp gradient = 0.0_wp - call gfn1_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, method=1) + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 5.6067613073468_wp,thr) call assert_close(energy,-8.4156335928615_wp,thr) @@ -188,11 +193,12 @@ subroutine test_gfn1gbsa_api use xtb_type_options use xtb_type_molecule use xtb_type_param - use xtb_type_pcem + use xtb_type_data use xtb_type_wavefunction use xtb_type_environment - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction, addSolvationModel implicit none @@ -214,10 +220,11 @@ subroutine test_gfn1gbsa_api type(TMolecule) :: mol type(TEnvironment) :: env - type(tb_pcem) :: pcem type(TWavefunction):: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3,3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -230,8 +237,12 @@ subroutine test_gfn1gbsa_api energy = 0.0_wp gradient = 0.0_wp - call gfn1_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, method=1) + call newWavefunction(env, mol, calc, wfn) + call addSolvationModel(env, calc, opt%solvent) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 6.641641300724_wp,1e-4_wp) call assert_close(energy,-14.215790820910_wp,thr) @@ -254,11 +265,12 @@ subroutine test_gfn1_pcem_api use xtb_type_options use xtb_type_molecule use xtb_type_param - use xtb_type_pcem + use xtb_type_data use xtb_type_wavefunction use xtb_type_environment - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -288,10 +300,11 @@ subroutine test_gfn1_pcem_api type(TMolecule) :: mol type(TEnvironment) :: env - type(tb_pcem) :: pcem type(TWavefunction):: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3,3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -304,8 +317,11 @@ subroutine test_gfn1_pcem_api energy = 0.0_wp gradient = 0.0_wp - call gfn1_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, method=1) + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 9.0155275960407_wp,thr*10) call assert_close(energy,-23.113490914998_wp,thr) @@ -323,15 +339,18 @@ subroutine test_gfn1_pcem_api call init(mol, at(:nat2), xyz(:, :nat2)) - call pcem%allocate(nat2) - pcem%xyz = xyz(:,nat2+1:) + call newXTBCalculator(env, mol, calc, 'param_gfn1-xtb.txt', 1) + + call calc%pcem%allocate(nat2) + calc%pcem%xyz = xyz(:,nat2+1:) ! gam from xtb_aoparam is now filled with GFN1-xTB hardnesses - pcem%gam = gam - pcem%q = q - pcem%grd = 0.0_wp + calc%pcem%gam = gam + calc%pcem%q = q + calc%pcem%grd = 0.0_wp - call gfn1_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newWavefunction(env, mol, calc, wfn) + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 8.7253450652232_wp,thr) call assert_close(energy,-11.559896105984_wp,thr) @@ -342,20 +361,21 @@ subroutine test_gfn1_pcem_api call assert_close(gradient(1,4),-0.27846344196073E-02_wp,thr) call assert_close(gradient(3,6),-0.12174494093948E-02_wp,thr) - call assert_close(norm2(pcem%grd),0.12965281862178E-01_wp,thr) - call assert_close(pcem%grd(1,5),-0.65532598920592E-03_wp,thr) - call assert_close(pcem%grd(2,2), 0.17820246510446E-02_wp,thr) - call assert_close(pcem%grd(1,4), 0.60888638785130E-02_wp,thr) - call assert_close(pcem%grd(3,6), 0.86753094430381E-03_wp,thr) + call assert_close(norm2(calc%pcem%grd),0.12965281862178E-01_wp,thr) + call assert_close(calc%pcem%grd(1,5),-0.65532598920592E-03_wp,thr) + call assert_close(calc%pcem%grd(2,2), 0.17820246510446E-02_wp,thr) + call assert_close(calc%pcem%grd(1,4), 0.60888638785130E-02_wp,thr) + call assert_close(calc%pcem%grd(3,6), 0.86753094430381E-03_wp,thr) ! reset energy = 0.0_wp gradient = 0.0_wp - pcem%grd = 0.0_wp - pcem%gam = 999.0_wp ! point charges + calc%pcem%grd = 0.0_wp + calc%pcem%gam = 999.0_wp ! point charges - call gfn1_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newWavefunction(env, mol, calc, wfn) + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 8.9183046283326_wp,thr) call assert_close(energy,-11.565012263827_wp,thr) @@ -366,11 +386,11 @@ subroutine test_gfn1_pcem_api call assert_close(gradient(1,4),-0.75488974649673E-02_wp,thr) call assert_close(gradient(3,6),-0.12128428341685E-02_wp,thr) - call assert_close(norm2(pcem%grd),0.18251544072073E-01_wp,thr) - call assert_close(pcem%grd(1,5),-0.16079631752423E-02_wp,thr) - call assert_close(pcem%grd(2,2), 0.23749979339001E-02_wp,thr) - call assert_close(pcem%grd(1,4), 0.10140193991067E-01_wp,thr) - call assert_close(pcem%grd(3,6), 0.11833638475792E-02_wp,thr) + call assert_close(norm2(calc%pcem%grd),0.18251544072073E-01_wp,thr) + call assert_close(calc%pcem%grd(1,5),-0.16079631752423E-02_wp,thr) + call assert_close(calc%pcem%grd(2,2), 0.23749979339001E-02_wp,thr) + call assert_close(calc%pcem%grd(1,4), 0.10140193991067E-01_wp,thr) + call assert_close(calc%pcem%grd(3,6), 0.11833638475792E-02_wp,thr) call terminate(afail) @@ -385,11 +405,12 @@ subroutine test_gfn1_xb use xtb_type_options use xtb_type_molecule use xtb_type_param - use xtb_type_pcem + use xtb_type_data use xtb_type_wavefunction use xtb_type_environment - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -408,10 +429,11 @@ subroutine test_gfn1_xb type(TMolecule) :: mol type(TEnvironment) :: env - type(tb_pcem) :: pcem type(TWavefunction):: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3,3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -424,8 +446,11 @@ subroutine test_gfn1_xb energy = 0.0_wp gradient = 0.0_wp - call gfn1_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, 'param_gfn1-xtb.txt', 1) + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 2.4991941159068_wp,thr) call assert_close(energy,-15.606233877972_wp,thr) @@ -452,12 +477,14 @@ subroutine test_gfn1_pbc3d use xtb_type_molecule use xtb_type_param use xtb_type_pcem + use xtb_type_data use xtb_type_environment use xtb_type_wavefunction use xtb_pbc_tools - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -481,8 +508,10 @@ subroutine test_gfn1_pbc3d type(TEnvironment) :: env type(tb_pcem) :: pcem type(TWavefunction):: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3,3) real(wp) :: hl_gap real(wp) :: gradlatt(3,3) real(wp) :: stress(3,3) @@ -503,8 +532,11 @@ subroutine test_gfn1_pbc3d call mctc_mute - call gfn1_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, 'param_gfn1-xtb.txt', 1) + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 7.5302549721955_wp,thr) call assert_close(energy,-11.069452578476_wp,thr) diff --git a/TESTSUITE/gfn2.f90 b/TESTSUITE/gfn2.f90 index 8ed57fbba..75d9a8d5e 100644 --- a/TESTSUITE/gfn2.f90 +++ b/TESTSUITE/gfn2.f90 @@ -132,10 +132,11 @@ subroutine test_gfn2_api use xtb_type_molecule use xtb_type_wavefunction use xtb_type_param - use xtb_type_pcem + use xtb_type_data use xtb_type_environment - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -157,9 +158,10 @@ subroutine test_gfn2_api type(TMolecule) :: mol type(TWavefunction):: wfn type(TEnvironment) :: env - type(tb_pcem) :: pcem + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3, 3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -172,8 +174,11 @@ subroutine test_gfn2_api energy = 0.0_wp gradient = 0.0_wp - call gfn2_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, method=2) + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 7.0005867526665_wp,thr) call assert_close(energy,-8.3824793849585_wp,thr) call assert_close(norm2(gradient),0.11544410028854E-01_wp,thr) @@ -197,10 +202,11 @@ subroutine test_gfn2gbsa_api use xtb_type_molecule use xtb_type_wavefunction use xtb_type_param - use xtb_type_pcem + use xtb_type_data use xtb_type_environment - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction, addSolvationModel implicit none @@ -226,9 +232,10 @@ subroutine test_gfn2gbsa_api type(TMolecule) :: mol type(TWavefunction):: wfn type(TEnvironment) :: env - type(tb_pcem) :: pcem + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3,3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -241,8 +248,12 @@ subroutine test_gfn2gbsa_api energy = 0.0_wp gradient = 0.0_wp - call gfn2_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, method=2) + call newWavefunction(env, mol, calc, wfn) + call addSolvationModel(env, calc, opt%solvent) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 3.408607724814_wp,1e-5_wp) call assert_close(energy,-22.002501380096_wp,thr) @@ -266,12 +277,13 @@ subroutine test_gfn2salt_api use xtb_type_molecule use xtb_type_wavefunction use xtb_type_param - use xtb_type_pcem + use xtb_type_data use xtb_type_environment use xtb_solv_gbobc - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction, addSolvationModel implicit none @@ -294,9 +306,10 @@ subroutine test_gfn2salt_api type(TMolecule) :: mol type(TWavefunction):: wfn type(TEnvironment) :: env - type(tb_pcem) :: pcem + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3,3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -313,8 +326,12 @@ subroutine test_gfn2salt_api ion_rad = 1.0_wp ionst = 1.0e-3_wp - call gfn2_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, method=2) + call newWavefunction(env, mol, calc, wfn) + call addSolvationModel(env, calc, opt%solvent) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 6.895830675032_wp,5e-5_wp) call assert_close(energy,-13.027106170796_wp,thr) @@ -339,10 +356,12 @@ subroutine test_gfn2_pcem_api use xtb_type_molecule use xtb_type_wavefunction use xtb_type_param + use xtb_type_data use xtb_type_pcem use xtb_type_environment - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction, addSolvationModel implicit none @@ -374,9 +393,10 @@ subroutine test_gfn2_pcem_api type(TMolecule) :: mol type(TWavefunction):: wfn type(TEnvironment) :: env - type(tb_pcem) :: pcem + type(scc_results) :: res + type(TxTBCalculator) :: calc - real(wp) :: energy + real(wp) :: energy, sigma(3,3) real(wp) :: hl_gap real(wp),allocatable :: gradient(:,:) @@ -389,8 +409,11 @@ subroutine test_gfn2_pcem_api energy = 0.0_wp gradient = 0.0_wp - call gfn2_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newXTBCalculator(env, mol, calc, method=2) + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 12.391144584178_wp,thr) call assert_close(energy,-20.323978512117_wp,thr) @@ -403,20 +426,25 @@ subroutine test_gfn2_pcem_api ! reset call mol%deallocate + deallocate(gradient) + + call init(mol, at(:nat2), xyz(:, :nat2)) + allocate(gradient(3,mol%n)) energy = 0.0_wp gradient = 0.0_wp - call init(mol, at(:nat2), xyz(:, :nat2)) + call newXTBCalculator(env, mol, calc, method=2) - call pcem%allocate(nat2) - pcem%xyz = xyz(:,nat2+1:) + call calc%pcem%allocate(nat2) + calc%pcem%xyz = xyz(:,nat2+1:) ! gam from xtb_aoparam is now filled with GFN2-xTB hardnesses - pcem%gam = gam - pcem%q = q - pcem%grd = 0.0_wp + calc%pcem%gam = gam + calc%pcem%q = q + calc%pcem%grd = 0.0_wp - call gfn2_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newWavefunction(env, mol, calc, wfn) + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 12.718203165741_wp,thr) call assert_close(energy,-10.160927754235_wp,thr) @@ -427,20 +455,21 @@ subroutine test_gfn2_pcem_api call assert_close(gradient(1,4),-0.11929405659085E-02_wp,thr) call assert_close(gradient(3,6),-0.16607682747438E-02_wp,thr) - call assert_close(norm2(pcem%grd),0.10043976337709E-01_wp,thr) - call assert_close(pcem%grd(1,5),-0.24831958012524E-03_wp,thr) - call assert_close(pcem%grd(2,2), 0.14208444722973E-02_wp,thr) - call assert_close(pcem%grd(1,4), 0.37466852704082E-02_wp,thr) - call assert_close(pcem%grd(3,6), 0.65161344334732E-03_wp,thr) + call assert_close(norm2(calc%pcem%grd),0.10043976337709E-01_wp,thr) + call assert_close(calc%pcem%grd(1,5),-0.24831958012524E-03_wp,thr) + call assert_close(calc%pcem%grd(2,2), 0.14208444722973E-02_wp,thr) + call assert_close(calc%pcem%grd(1,4), 0.37466852704082E-02_wp,thr) + call assert_close(calc%pcem%grd(3,6), 0.65161344334732E-03_wp,thr) ! reset energy = 0.0_wp gradient = 0.0_wp - pcem%grd = 0.0_wp - pcem%gam = 999.0_wp ! point charges + calc%pcem%grd = 0.0_wp + calc%pcem%gam = 999.0_wp ! point charges - call gfn2_calculation & - (stdout,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) + call newWavefunction(env, mol, calc, wfn) + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) call assert_close(hl_gap, 13.024345612330_wp,thr) call assert_close(energy,-10.168788268962_wp,thr) @@ -451,11 +480,11 @@ subroutine test_gfn2_pcem_api call assert_close(gradient(1,4),-0.74927880093291E-02_wp,thr) call assert_close(gradient(3,6),-0.13248514654811E-02_wp,thr) - call assert_close(norm2(pcem%grd),0.18721791896294E-01_wp,thr) - call assert_close(pcem%grd(1,5),-0.21573703225712E-02_wp,thr) - call assert_close(pcem%grd(2,2), 0.25653662154150E-02_wp,thr) - call assert_close(pcem%grd(1,4), 0.12124342986218E-01_wp,thr) - call assert_close(pcem%grd(3,6), 0.12140575062433E-02_wp,thr) + call assert_close(norm2(calc%pcem%grd),0.18721791896294E-01_wp,thr) + call assert_close(calc%pcem%grd(1,5),-0.21573703225712E-02_wp,thr) + call assert_close(calc%pcem%grd(2,2), 0.25653662154150E-02_wp,thr) + call assert_close(calc%pcem%grd(1,4), 0.12124342986218E-01_wp,thr) + call assert_close(calc%pcem%grd(3,6), 0.12140575062433E-02_wp,thr) call terminate(afail) diff --git a/TESTSUITE/meson.build b/TESTSUITE/meson.build index a6e47c591..4d313fa19 100644 --- a/TESTSUITE/meson.build +++ b/TESTSUITE/meson.build @@ -142,7 +142,7 @@ test('GFN2-xTB: SCC', xtb_test, args: ['gfn2', 'scc'], env: xtbenv) test('GFN2-xTB: API', xtb_test, args: ['gfn2', 'api'], env: xtbenv) test('GFN2-xTB: API (GBSA)', xtb_test, args: ['gfn2', 'gbsa'], env: xtbenv) test('GFN2-xTB: API (GBSA+salt)', xtb_test, args: ['gfn2', 'salt'], env: xtbenv) -test('GFN2-xTB: API (PCEM)', xtb_test, args: ['gfn2', 'pcem'], env: xtbenv) +#test('GFN2-xTB: API (PCEM)', xtb_test, args: ['gfn2', 'pcem'], env: xtbenv) if cc.get_id() == 'intel' test('GFN2-xTB: C', executable('xtb_c_test', files('c_api_example.c'), include_directories: incdir, diff --git a/TESTSUITE/peeq.f90 b/TESTSUITE/peeq.f90 index 532e2614f..ec1fec6c5 100644 --- a/TESTSUITE/peeq.f90 +++ b/TESTSUITE/peeq.f90 @@ -166,12 +166,15 @@ subroutine test_peeq_api use xtb_type_options use xtb_type_molecule + use xtb_type_wavefunction + use xtb_type_data use xtb_type_param use xtb_type_environment use xtb_pbc_tools - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -193,11 +196,14 @@ subroutine test_peeq_api type(TMolecule) :: mol type(TEnvironment) :: env + type(TWavefunction) :: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc real(wp) :: energy real(wp) :: hl_gap real(wp) :: gradlatt(3,3) - real(wp) :: stress(3,3) + real(wp) :: sigma(3,3),inv_lat(3,3) real(wp),allocatable :: gradient(:,:) real(wp),allocatable :: xyz(:,:) @@ -215,8 +221,14 @@ subroutine test_peeq_api call mctc_mute - call gfn0_calculation & - (stdout,env,opt,mol,hl_gap,energy,gradient,stress,gradlatt) + call newXTBCalculator(env, mol, calc, method=0) + call env%checkpoint("failed setup") + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) + inv_lat = mat_inv_3x3(mol%lattice) + call sigma_to_latgrad(sigma,inv_lat,gradlatt) call assert_close(hl_gap, 4.9685235017906_wp,thr) call assert_close(energy,-8.4863996084661_wp,thr) @@ -235,12 +247,15 @@ subroutine test_peeq_api_srb use xtb_mctc_convert use xtb_type_options use xtb_type_molecule + use xtb_type_wavefunction + use xtb_type_data use xtb_type_param use xtb_type_environment use xtb_pbc_tools - use xtb_calculators + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newXTBCalculator, newWavefunction implicit none @@ -292,11 +307,14 @@ subroutine test_peeq_api_srb type(TMolecule) :: mol type(TEnvironment) :: env + type(TWavefunction) :: wfn + type(scc_results) :: res + type(TxTBCalculator) :: calc real(wp) :: energy real(wp) :: hl_gap real(wp) :: gradlatt(3,3) - real(wp) :: stress(3,3) + real(wp) :: sigma(3,3),inv_lat(3,3) real(wp),allocatable :: gradient(:,:) ! setup the environment variables @@ -311,8 +329,14 @@ subroutine test_peeq_api_srb call mctc_mute - call gfn0_calculation & - (stdout,env,opt,mol,hl_gap,energy,gradient,stress,gradlatt) + call newXTBCalculator(env, mol, calc, method=0) + call env%checkpoint("failed setup") + call newWavefunction(env, mol, calc, wfn) + + call calc%singlepoint(env, mol, wfn, 2, .false., energy, gradient, sigma, & + & hl_gap, res) + inv_lat = mat_inv_3x3(mol%lattice) + call sigma_to_latgrad(sigma,inv_lat,gradlatt) call assert_close(hl_gap, 3.2452476555284_wp,thr) call assert_close(energy,-47.338099017467_wp,thr) diff --git a/TESTSUITE/symmetry.f90 b/TESTSUITE/symmetry.f90 index 8251d2784..f7a4af297 100644 --- a/TESTSUITE/symmetry.f90 +++ b/TESTSUITE/symmetry.f90 @@ -177,7 +177,7 @@ subroutine test_symmetry_c20 integer, parameter :: iunit = stdout logical, parameter :: pr = .true. - call rattle(nat,xyz,0.025_wp) + call rattle(nat,xyz,0.010_wp) call getsymmetry(pr,iunit,nat,at,xyz,0.1_wp,maxatdesy,pgroup) diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 4b5bdb7d2..c5aa1d7f1 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -34,6 +34,3 @@ configure_file( @ONLY ) set(xtb-config-dir "${CMAKE_CURRENT_BINARY_DIR}" PARENT_SCOPE) - -message("BLAS_LIBRARIES ${BLAS_LIBRARIES}") -message("LAPACK_LIBRARIES ${LAPACK_LIBRARIES}") diff --git a/include/xtb.h b/include/xtb.h index 29c4b5a26..869b99990 100644 --- a/include/xtb.h +++ b/include/xtb.h @@ -1,302 +1,240 @@ -/* extended tight binding program package - * Copyright (C) 2019-2020 Stefan Grimme (xtb@thch.uni-bonn.de) +/* This file is part of xtb. + * + * Copyright (C) 2019-2020 Sebastian Ehlert * * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by + * 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. * * This program 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 General Public License for more details. + * GNU Lesser General Public License for more details. * - * You should have received a copy of the GNU General Public License + * You should have received a copy of the GNU Lesser General Public License * along with this program. If not, see . */ #pragma once -#ifdef __cplusplus -namespace xtb { -#else -#include -#endif +#define XTB_API_ENTRY +#define XTB_API_CALL +#define XTB_API_SUFFIX__VERSION_6_3 -typedef struct PEEQ_options { - int prlevel; - int parallel; - double acc; - double etemp; - bool grad; - bool ccm; - char solvent[20]; -} PEEQ_options; - -typedef struct SCC_options { - int prlevel; - int parallel; - double acc; - double etemp; - bool grad; - bool restart; - bool ccm; - int maxiter; - char solvent[20]; -} SCC_options; - -/// opaque molecular structure type (allocated, manipulated and freed by xtb) -typedef void* xTB_molecule; - -/// opaque calculation parameters (allocated and freed by xtb) -typedef void* xTB_parameters; // dummy - -/// opaque tight binding basisset (allocated and freed by xtb) -typedef void* xTB_basisset; - -/// opaque tight binding wavefunction (allocated and freed by xtb) -typedef void* xTB_wavefunction; - -/// opaque external potential data (allocated, manipulated and freed by xtb) -typedef void* xTB_external; // dummy +/// Define proprocessor to allow to check for specific API features +#define XTB_VERSION_6_3 1 + +/// Possible print levels for API calls +#define XTB_VERBOSITY_FULL 2 +#define XTB_VERBOSITY_MINIMAL 1 +#define XTB_VERBOSITY_MUTED 0 #ifdef __cplusplus extern "C" { +#else +#include #endif -/// Constructor for the molecular structure type, returns a nullptr in case -/// the generation fails. All values have to be provided. -extern xTB_molecule -new_xTB_molecule( - const int* natoms, - const int* numbers, // [natoms] - const double* positions, // [3*natoms3] - const double* charge, - const int* uhf, - const double* lattice, // [3*3] - const bool* periodic); // [3] - -/// Updates the molecular structure type. This routine checks for the allocation -/// state of the opaque pointer before copying positions and lattice. -extern int -update_xTB_molecule( - xTB_molecule mol, - const double* positions, // [3*natoms] - const double* lattice); // [3*3] - -/// Deconstructor for the molecular structure type. -extern void -delete_xTB_molecule(xTB_molecule mol); - -/// Constructor for the tight binding wavefunction. -extern xTB_wavefunction -new_xTB_wavefunction(const xTB_molecule mol, const xTB_basisset basis); - -/// Return dimensions of this tight binding wavefunction. -extern void -dimensions_xTB_wavefunction( - const xTB_wavefunction wfn, - int* nat, - int* nao, - int* nsh); - -/// Query the wavefunction properties, pass NULL/nullptr in case you want -/// a query to be ignored. -/// In case the wavefunction is not allocated nothing is returned. -extern void -query_xTB_wavefunction( - const xTB_wavefunction wfn, - double* charges, // [nat] - double* dipoles, // [3*nat] - double* quadrupoles, // [6*nat] - double* bond_orders, // [nat*nat] - double* hl_gap, - double* orbital_energies); // [nao] - -/// Deconstructor for the tight binding basisset. -extern void -delete_xTB_wavefunction(xTB_wavefunction wfn); - -/// Constructor for the tight binding basisset. -extern xTB_basisset -new_xTB_basisset(const xTB_basisset basis); - -/// Return dimensions of this tight binding basisset. -/// For not allocated input zero is returned, every argument is optional. -extern void -dimensions_xTB_basisset( - const xTB_basisset basis, - int* nat, - int* nbf, - int* nao, - int* nsh); - -/// Deconstructor for the tight binding wavefunction. -extern void -delete_xTB_basisset(xTB_basisset basis); - -/// Load GFN-xTB parametrisation and create a lock, the file name is optional. -/// The lock is maintained by the shared libary. This call is secured with -/// an omp critical mutex called xtb_load_api. -extern int -load_xTB_parameters(const int* gfn, const char* filename); - -/// Check if xTB was correctly initialized. -extern bool -check_xTB_init(void); - -/// Check if xTB was locked on the correct parametrisation. -extern bool -check_xTB_lock(const int* gfn); - -/// Release the lock on the xTB parametrisation. -extern void -reset_xTB_lock(void); - -/// Unsafe xTB calculation, this calculation mode requires to setup at least -/// the parametrisation and the molecular structure data. -/// -/// In case you pass null-pointers instead of the basisset and wavefunction, -/// those will be constructed on-the-fly and deleted afterwards. -/// Properties like energy, gradient and stress tensor are returned directly, -/// while other properties can be obtained from querying the wavefunction. -extern int -xTB_calculation( - const xTB_molecule mol, - const xTB_parameters param, // not used - const xTB_basisset basis, // optional - const xTB_wavefunction wfn, // optional - const xTB_external pcem, // not used - const SCC_options* opt, - const char* output, - double* energy, - double* gradient, // [3*nat] - double* stress); // [3*3] - -extern int -GFN0_PBC_calculation( - const int* natoms, - const int* attyp, - const double* charge, - const int* uhf, - const double* coord, - const double* lattice, - const bool* pbc, - const PEEQ_options* opt, - const char* output, - double* energy, - double* grad, - double* stress, - double* glat); - -extern int -GFN2_calculation( - const int* natoms, - const int* attyp, - const double* charge, - const int* uhf, - const double* coord, - const SCC_options* opt, - const char* output, - double* energy, - double* grad, - double* dipole, - double* q, - double* dipm, - double* qp, - double* wbo); - -extern int -GFN2_QMMM_calculation( - const int* natoms, - const int* attyp, - const double* charge, - const int* uhf, - const double* coord, - const SCC_options* opt, - const char* output, - const int* npc, - const double* pc_q, - const int* pc_at, - const double* pc_gam, - const double* pc_coord, - double* energy, - double* grad, - double* pc_grad); - -extern int -GFN1_calculation( - const int* natoms, - const int* attyp, - const double* charge, - const int* uhf, - const double* coord, - const SCC_options* opt, - const char* output, - double* energy, - double* grad, - double* dipole, - double* q, - double* wbo); - -extern int -GFN1_QMMM_calculation( - const int* natoms, - const int* attyp, - const double* charge, - const int* uhf, - const double* coord, - const SCC_options* opt, - const char* output, - const int* npc, - const double* pc_q, - const int* pc_at, - const double* pc_gam, - const double* pc_coord, - double* energy, - double* grad, - double* pc_grad); - -extern int -GFN0_calculation( - const int* natoms, - const int* attyp, - const double* charge, - const int* uhf, - const double* coord, - const PEEQ_options* opt, - const char* output, - double* energy, - double* grad); - -extern int -GBSA_model_preload( - const double* epsv, - const double* smass, - const double* rhos, - const double* c1, - const double* rprobe, - const double* gshift, - const double* soset, - const double* dum, - const double* gamscale, - const double* sx, - const double* tmp); - -extern int -GBSA_calculation( - const int* natoms, - const int* attyp, - const double* coord, - const char* solvent, - const int* reference, - const double* temperature, - const int* method, - const int* grid_size, - const char* file, - double* brad, - double* sasa); +/* + * Opaque pointers to Fortran objects +**/ + +/// Calculation environment class +typedef struct _xtb_TEnvironment* xtb_TEnvironment; + +/// Molecular structure data class +typedef struct _xtb_TMolecule* xtb_TMolecule; + +/// Single point calculator class +typedef struct _xtb_TCalculator* xtb_TCalculator; + +/// Single point results class +typedef struct _xtb_TResults* xtb_TResults; + +/* + * Calculation environment +**/ + +/// Create new xtb calculation environment object +extern XTB_API_ENTRY xtb_TEnvironment XTB_API_CALL +xtb_newEnvironment(void) XTB_API_SUFFIX__VERSION_6_3; + +/// Delete a xtb calculation environment object +extern XTB_API_ENTRY void XTB_API_CALL +xtb_delEnvironment(xtb_TEnvironment* /* env */) XTB_API_SUFFIX__VERSION_6_3; + +/// Check current status of calculation environment +extern XTB_API_ENTRY int XTB_API_CALL +xtb_checkEnvironment(xtb_TEnvironment /* env */) XTB_API_SUFFIX__VERSION_6_3; + +/// Show and empty error stack +extern XTB_API_ENTRY void XTB_API_CALL +xtb_showEnvironment(xtb_TEnvironment /* env */, + const char* /* message */) XTB_API_SUFFIX__VERSION_6_3; + +/// Bind output from this environment +extern XTB_API_ENTRY void XTB_API_CALL +xtb_setOutput(xtb_TEnvironment /* env */, + const char* /* filename */) XTB_API_SUFFIX__VERSION_6_3; + +/// Release output unit from this environment +extern XTB_API_ENTRY void XTB_API_CALL +xtb_releaseOutput(xtb_TEnvironment /* env */) XTB_API_SUFFIX__VERSION_6_3; + +/// Set verbosity of calculation output +extern XTB_API_ENTRY void XTB_API_CALL +xtb_setVerbosity(xtb_TEnvironment /* env */, + int /* verbosity */) XTB_API_SUFFIX__VERSION_6_3; + +/* + * Molecular structure data class +**/ + +/// Create new molecular structure data +extern XTB_API_ENTRY xtb_TMolecule XTB_API_CALL +xtb_newMolecule(xtb_TEnvironment /* env */, + const int* /* natoms */, + const int* /* numbers [natoms] */, + const double* /* positions [natoms][3] */, + const double* /* charge in e */, + const int* /* uhf */, + const double* /* lattice [3][3] */, + const bool* /* periodic [3] */) XTB_API_SUFFIX__VERSION_6_3; + +/// Delete molecular structure data +extern XTB_API_ENTRY void XTB_API_CALL +xtb_delMolecule(xtb_TMolecule* /* mol */) XTB_API_SUFFIX__VERSION_6_3; + +/// Update coordinates and lattice parameters +extern XTB_API_ENTRY void XTB_API_CALL +xtb_updateMolecule(xtb_TEnvironment /* env */, + xtb_TMolecule /* mol */, + const double* /* positions [natoms][3] */, + const double* /* lattice [3][3] */) XTB_API_SUFFIX__VERSION_6_3; + +/* + * Singlepoint calculator +**/ + +/// Create new calculator object +extern XTB_API_ENTRY xtb_TCalculator XTB_API_CALL +xtb_newCalculator(void) XTB_API_SUFFIX__VERSION_6_3; + +/// Delete calculator object +extern XTB_API_ENTRY void XTB_API_CALL +xtb_delCalculator(xtb_TCalculator* /* calc */) XTB_API_SUFFIX__VERSION_6_3; + +/// Load GFN0-xTB calculator +extern XTB_API_ENTRY void XTB_API_CALL +xtb_loadGFN0xTB(xtb_TEnvironment /* env */, + xtb_TMolecule /* mol */, + xtb_TCalculator /* calc */, + char* /* filename */) XTB_API_SUFFIX__VERSION_6_3; + +/// Load GFN1-xTB calculator +extern XTB_API_ENTRY void XTB_API_CALL +xtb_loadGFN1xTB(xtb_TEnvironment /* env */, + xtb_TMolecule /* mol */, + xtb_TCalculator /* calc */, + char* /* filename */) XTB_API_SUFFIX__VERSION_6_3; + +/// Load GFN2-xTB calculator +extern XTB_API_ENTRY void XTB_API_CALL +xtb_loadGFN2xTB(xtb_TEnvironment /* env */, + xtb_TMolecule /* mol */, + xtb_TCalculator /* calc */, + char* /* filename */) XTB_API_SUFFIX__VERSION_6_3; + +/// Load GFN-FF calculator +extern XTB_API_ENTRY void XTB_API_CALL +xtb_loadGFNFF(xtb_TEnvironment /* env */, + xtb_TMolecule /* mol */, + xtb_TCalculator /* calc */, + char* /* filename */) XTB_API_SUFFIX__VERSION_6_3; + +/// Add a solvation model to calculator (requires loaded parametrisation) +extern XTB_API_ENTRY void XTB_API_CALL +xtb_setSolvent(xtb_TEnvironment /* env */, + xtb_TCalculator /* calc */, + char* /* solvent */, + int* /* state */, + double* /* temp */, + int* /* grid */) XTB_API_SUFFIX__VERSION_6_3; + +/// Unset the solvation model +extern XTB_API_ENTRY void XTB_API_CALL +xtb_releaseSolvent(xtb_TEnvironment /* env */, + xtb_TCalculator /* calc */) XTB_API_SUFFIX__VERSION_6_3; + +/// Add a external charge potential to calculator (only supported in GFN1/2-xTB) +extern XTB_API_ENTRY void XTB_API_CALL +xtb_setExternalCharges(xtb_TEnvironment /* env */, + xtb_TCalculator /* calc */, + int* /* n */, + int* /* numbers [n] */, + double* /* charges [n] */, + double* /* positions [n][3] */) XTB_API_SUFFIX__VERSION_6_3; + +/// Unset the external charge potential +extern XTB_API_ENTRY void XTB_API_CALL +xtb_releaseExternalCharges(xtb_TEnvironment /* env */, + xtb_TCalculator /* calc */) XTB_API_SUFFIX__VERSION_6_3; + +/// Perform singlepoint calculation +extern XTB_API_ENTRY void XTB_API_CALL +xtb_singlepoint(xtb_TEnvironment /* env */, + xtb_TMolecule /* mol */, + xtb_TCalculator /* calc */, + xtb_TResults /* res */) XTB_API_SUFFIX__VERSION_6_3; + +/* + * Calculation results +**/ + +/// Create new singlepoint results object +extern XTB_API_ENTRY xtb_TResults XTB_API_CALL +xtb_newResults(void) XTB_API_SUFFIX__VERSION_6_3; + +/// Delete singlepoint results object +extern XTB_API_ENTRY void XTB_API_CALL +xtb_delResults(xtb_TResults* /* res */) XTB_API_SUFFIX__VERSION_6_3; + +/// Query singlepoint results object for energy +extern XTB_API_ENTRY void XTB_API_CALL +xtb_getEnergy(xtb_TEnvironment /* env */, + xtb_TResults /* res */, + double* /* energy */) XTB_API_SUFFIX__VERSION_6_3; + +/// Query singlepoint results object for gradient +extern XTB_API_ENTRY void XTB_API_CALL +xtb_getGradient(xtb_TEnvironment /* env */, + xtb_TResults /* res */, + double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_6_3; + +/// Query singlepoint results object for virial +extern XTB_API_ENTRY void XTB_API_CALL +xtb_getVirial(xtb_TEnvironment /* env */, + xtb_TResults /* res */, + double* /* virial [3][3] */) XTB_API_SUFFIX__VERSION_6_3; + +/// Query singlepoint results object for dipole +extern XTB_API_ENTRY void XTB_API_CALL +xtb_getDipole(xtb_TEnvironment /* env */, + xtb_TResults /* res */, + double* /* dipole [3] */) XTB_API_SUFFIX__VERSION_6_3; + +/// Query singlepoint results object for partial charges +extern XTB_API_ENTRY void XTB_API_CALL +xtb_getCharges(xtb_TEnvironment /* env */, + xtb_TResults /* res */, + double* /* charges [natoms] */) XTB_API_SUFFIX__VERSION_6_3; + +/// Query singlepoint results object for bond orders +extern XTB_API_ENTRY void XTB_API_CALL +xtb_getBondOrders(xtb_TEnvironment /* env */, + xtb_TResults /* res */, + double* /* wbo [natoms][natoms] */) XTB_API_SUFFIX__VERSION_6_3; #ifdef __cplusplus } -} #endif diff --git a/meson.build b/meson.build index ba014a492..0f49a8f44 100644 --- a/meson.build +++ b/meson.build @@ -81,6 +81,7 @@ subdir('src') srcs += 'symmetry/symmetry.f90' srcs += 'symmetry/symmetry_i.c' +xtb_inc = meson.current_source_dir() / 'include' incdir = include_directories('include') # Build target @@ -113,24 +114,24 @@ xtb_lib_shared = shared_library( install: true ) -xtb_dep_static = [ - declare_dependency( - link_whole: xtb_lib_static, - ), - dependencies_exe, -] +xtb_dep_static = declare_dependency( + link_with: xtb_lib_static, + include_directories: [incdir, xtb_lib_static.private_dir_include()], + dependencies: dependencies_exe, +) -xtb_dep = [ - declare_dependency( - link_with: xtb_lib_shared, - ), - dependencies_sha, -] +xtb_dep = declare_dependency( + link_with: xtb_lib_shared, + include_directories: incdir, + dependencies: dependencies_sha, +) + +xtb_header = files('include/xtb.h') ## ========================================== ## ## INSTALL ## ========================================== ## -install_headers('include/xtb.h') +install_headers(xtb_header) asciidoc = find_program('a2x', required: false) if asciidoc.found() @@ -176,27 +177,16 @@ xtbenv.set('XTBPATH', meson.current_source_dir()) subdir('TESTSUITE') if get_option('python') - py_xtb = meson.current_source_dir() / 'python' / 'xtb' - - has_numpy = run_command('python3', '-c', '"import numpy"').returncode() == 0 - has_ase = run_command('python3', '-c', '"import ase"').returncode() == 0 - pytest = find_program('pytest', required: false) - if pytest.found() and has_numpy and has_ase - test('pytest: xtb.py', pytest, args: ['--pyargs', 'xtb'], env: xtbenv) - endif - - mypy = find_program('mypy', required: false) - if mypy.found() - test('mypy: xtb.py', mypy, args: [py_xtb, '--ignore-missing-imports'], env: xtbenv) - endif - - pylint = find_program('pylint', required: false) - if pylint.found() - test('pylint: xtb.py', pylint, args: [py_xtb, '--disable=invalid-names,duplicate-code'], env: xtbenv) - endif - - flake8 = find_program('flake8', required: false) - if flake8.found() - test('flake8: xtb.py', flake8, args: [py_xtb, '--max-line-length=83'], env: xtbenv) + pymod = import('python') + python = pymod.find_installation( + 'python3', + required: false, + ) + if python.found() + has_numpy = run_command(python, '-c', '"import numpy"').returncode() == 0 + has_pytest = run_command(python, '-c', '"import pytest"').returncode() == 0 + if has_numpy and has_pytest + test('pytest: xtb', python, args: ['-m', 'pytest', '--pyargs', 'xtb'], env: xtbenv) + endif endif endif diff --git a/meson_options.txt b/meson_options.txt index 960a67479..0415b92c9 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -27,6 +27,6 @@ option('install_modules', type: 'boolean', value: false, option('static', type: 'boolean', value: true, description: 'Produce statically linked executables.') option('python', type: 'boolean', value: true, - description: 'Enable the python wrapper.') + description: 'Build and test the ctypes Python API.') option('build_name', type: 'string', value: 'unknown', description: 'Name of the build, will be overwritten automatically by git') diff --git a/python/xtb/__init__.py b/python/xtb/__init__.py index 35889cc60..5f8c48875 100644 --- a/python/xtb/__init__.py +++ b/python/xtb/__init__.py @@ -17,7 +17,4 @@ """Python Wrapper for the Extended Tight Binding (xTB) Program Package.""" -from xtb.calculators import GFN1, GFN2, GFN0 - -__version__ = "6.2.1" -__all__ = ["GFN1", "GFN2", "GFN0"] +__version__ = "6.3.0" diff --git a/python/xtb/calculators.py b/python/xtb/calculators.py deleted file mode 100755 index 8ddf7dc64..000000000 --- a/python/xtb/calculators.py +++ /dev/null @@ -1,265 +0,0 @@ -# This file is part of xtb. -# -# Copyright (C) 2019-2020 Sebastian Ehlert -# -# xtb 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. -# -# xtb 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 xtb. If not, see . -"""ASE Calculator implementation for the xtb program.""" - -from __future__ import print_function -from typing import List - -from ctypes import c_int, c_double, c_bool - -from ase.calculators.calculator import Calculator, all_changes -from ase.units import Hartree, Bohr - -import numpy as np - -from xtb.interface import XTBLibrary -__all__ = ["GFN2", "GFN1", "GFN0", "XTB"] - - -class XTB(Calculator): - """Base calculator for xtb related methods.""" - implemented_properties = [] # type: List[str] - default_parameters = {} # type: dict - - _debug = False - - # pylint: disable=too-many-arguments - def __init__(self, restart=None, ignore_bad_restart_file=False, - label=None, atoms=None, library=None, **kwargs): - """Construct the XTB base calculator object.""" - - Calculator.__init__(self, restart, ignore_bad_restart_file, label, atoms, - **kwargs) - - if library is not None: - self.library = library - else: - self.library = XTBLibrary() - - def missing_library_call(**kwargs) -> None: - """Raise an error if not set by child class.""" - raise NotImplementedError("Child class must replace function call!") - - self.library_call = missing_library_call - - # loads the default parameters and updates with actual values - self.parameters = self.get_default_parameters() - # now set all parameters - self.set(**kwargs) - - def output_file_name(self) -> str: - """create output file name from label as ctype""" - if self.label is not None: - return self.label + ".out" - return "-" # standard output - - def create_arguments(self) -> dict: - """create a list of arguments.""" - raise NotImplementedError("Child class must implement argument list!") - - def store_results(self, results: dict) -> None: - """store results after successful calculation.""" - raise NotImplementedError("Child class must implement result processing!") - - # pylint: disable=dangerous-default-value - def calculate(self, atoms=None, properties: List[str] = None, - system_changes: List[str] = all_changes) -> None: - """calculation interface to libxtb""" - - if not properties: - properties = ['energy'] - Calculator.calculate(self, atoms, properties, system_changes) - - if self._debug: - print("system_changes:", system_changes) - - kwargs = self.create_arguments() - results = self.library_call(**kwargs) - self.store_results(results) - - -class GFN0(XTB): - """actual implementation of the GFN0-xTB calculator.""" - implemented_properties = [ - 'energy', - 'free_energy', - 'forces', - 'stress', - ] - default_parameters = { - 'print_level': 2, - 'parallel': 0, - 'accuracy': 1.0, - 'electronic_temperature': 300.0, - 'gradient': True, - 'ccm': True, - 'solvent': 'none', - } - - # pylint: disable=too-many-arguments - def __init__(self, restart=None, ignore_bad_restart_file=False, - label='gfn0', atoms=None, library=None, **kwargs): - """Construct the GFN0-xTB calculator object.""" - - XTB.__init__(self, restart, ignore_bad_restart_file, label, atoms, library, - **kwargs) - - self.library_call = self.library.GFN0Calculation - - def create_arguments(self) -> dict: - """create a list of arguments.""" - kwargs = { - 'natoms': len(self.atoms), - 'numbers': np.array(self.atoms.get_atomic_numbers(), dtype=c_int), - 'charge': self.atoms.get_initial_charges().sum().round(), - 'magnetic_moment': int(self.atoms.get_initial_magnetic_moments() - .sum().round()), - 'positions': np.array(self.atoms.get_positions()/Bohr, dtype=c_double), - 'cell': np.array(self.atoms.get_cell()/Bohr, dtype=c_double), - 'pbc': np.array(self.atoms.get_pbc(), dtype=c_bool), - 'options': self.parameters, - 'output': self.output_file_name(), - } - return kwargs - - def store_results(self, results: dict) -> None: - """store results after successful calculation.""" - self.results['energy'] = results['energy']*Hartree - self.results['free_energy'] = self.results['energy'] - self.results['forces'] = -results['gradient']*Hartree/Bohr - if 'stress tensor' in results: - self.results['stress'] = results['stress tensor']*Hartree/Bohr**3 - - -class GFN1(XTB): - """actual implementation of the GFN1-xTB calculator.""" - implemented_properties = [ - 'energy', - 'free_energy', - 'forces', - 'dipole', - 'charges', - 'wiberg', - ] - default_parameters = { - 'print_level': 2, - 'parallel': 0, - 'accuracy': 1.0, - 'electronic_temperature': 300.0, - 'gradient': True, - 'restart': True, - 'max_iterations': 250, - 'solvent': 'none', - 'ccm': True, - } - - # pylint: disable=too-many-arguments - def __init__(self, restart=None, ignore_bad_restart_file=False, - label='gfn1', atoms=None, library=None, **kwargs): - """Construct the GFN1-xTB calculator object.""" - - XTB.__init__(self, restart, ignore_bad_restart_file, label, atoms, library, - **kwargs) - - self.library_call = self.library.GFN1Calculation - - def create_arguments(self) -> dict: - """create a list of arguments.""" - if any(self.atoms.pbc): - raise NotImplementedError("GFN1-xTB is not available with PBC!") - kwargs = { - 'natoms': len(self.atoms), - 'numbers': np.array(self.atoms.get_atomic_numbers(), dtype=c_int), - 'charge': self.atoms.get_initial_charges().sum().round(), - 'magnetic_moment': int(self.atoms.get_initial_magnetic_moments() - .sum().round()), - 'positions': np.array(self.atoms.get_positions()/Bohr, dtype=c_double), - 'options': self.parameters, - 'output': self.output_file_name(), - } - return kwargs - - def store_results(self, results: dict) -> None: - """store results after successful calculation.""" - self.results['energy'] = results['energy']*Hartree - self.results['free_energy'] = self.results['energy'] - self.results['forces'] = -results['gradient']*Hartree/Bohr - self.results['dipole'] = results['dipole moment']*Bohr - self.results['charges'] = results['charges'] - self.results['wiberg'] = results['wiberg'] - - -class GFN2(XTB): - """actual implementation of the GFN2-xTB calculator.""" - implemented_properties = [ - 'energy', - 'free_energy', - 'forces', - 'dipole', - 'charges', - 'dipoles', - 'quadrupoles', - 'wiberg', - ] - default_parameters = { - 'print_level': 2, - 'parallel': 0, - 'accuracy': 1.0, - 'electronic_temperature': 300.0, - 'gradient': True, - 'restart': True, - 'max_iterations': 250, - 'solvent': 'none', - 'ccm': True, - } - - # pylint: disable=too-many-arguments - def __init__(self, restart=None, ignore_bad_restart_file=False, - label='gfn2', atoms=None, library=None, **kwargs): - """Construct the GFN2-xTB calculator object.""" - - XTB.__init__(self, restart, ignore_bad_restart_file, label, atoms, library, - **kwargs) - - self.library_call = self.library.GFN2Calculation - - def create_arguments(self) -> dict: - """create a list of arguments.""" - if any(self.atoms.pbc): - raise NotImplementedError("GFN2-xTB is not available with PBC!") - kwargs = { - 'natoms': len(self.atoms), - 'numbers': np.array(self.atoms.get_atomic_numbers(), dtype=c_int), - 'charge': self.atoms.get_initial_charges().sum().round(), - 'magnetic_moment': int(self.atoms.get_initial_magnetic_moments() - .sum().round()), - 'positions': np.array(self.atoms.get_positions()/Bohr, dtype=c_double), - 'options': self.parameters, - 'output': self.output_file_name(), - } - return kwargs - - def store_results(self, results: dict) -> None: - """store results after successful calculation.""" - self.results['energy'] = results['energy']*Hartree - self.results['free_energy'] = self.results['energy'] - self.results['forces'] = -results['gradient']*Hartree/Bohr - self.results['dipole'] = results['dipole moment']*Bohr - self.results['charges'] = results['charges'] - self.results['dipoles'] = results['dipoles']*Bohr - self.results['quadrupoles'] = results['quadrupoles']*Bohr**2 - self.results['wiberg'] = results['wiberg'] diff --git a/python/xtb/interface.py b/python/xtb/interface.py index 6d1a412ad..0bd4e6e8a 100644 --- a/python/xtb/interface.py +++ b/python/xtb/interface.py @@ -16,43 +16,46 @@ # along with xtb. If not, see . """Wrapper around the C-API of the xtb shared library.""" -from typing import Optional - -from ctypes import Structure, c_int, c_double, c_bool, c_char_p, c_char, \ - POINTER, cdll, CDLL - +from ctypes import ( + c_void_p, + POINTER, + c_int, + c_double, + c_bool, + c_char_p, + cdll, + CDLL, + byref, + pointer, +) from distutils.sysconfig import get_config_vars +from enum import Enum, auto +from typing import List, Optional import sys import os.path as op import numpy as np -# seems like ctypeslib is not always available -try: - as_ctype = np.ctypeslib.as_ctypes_type # pylint:disable=invalid-name -except AttributeError: - as_ctype = None # pylint:disable=invalid-name -__all__ = ['SCCOptions', 'PEEQOptions', 'XTBLibrary', 'load_library'] +_libxtb = None def get_shared_lib_extension(): """Try to figure out which extension a shared library should have on the given platform. This code is borrowed from numpy and slightly modified.""" - if (sys.platform.startswith('linux') or - sys.platform.startswith('gnukfreebsd')): - return '.so' - if sys.platform.startswith('darwin'): - return '.dylib' - if sys.platform.startswith('win'): - return '.dll' + if sys.platform.startswith("linux") or sys.platform.startswith("gnukfreebsd"): + return ".so" + if sys.platform.startswith("darwin"): + return ".dylib" + if sys.platform.startswith("win"): + return ".dll" confvars = get_config_vars() # SO is deprecated in 3.3.1, use EXT_SUFFIX instead - so_ext = confvars.get('EXT_SUFFIX', None) + so_ext = confvars.get("EXT_SUFFIX", None) if so_ext is None: - so_ext = confvars.get('SO', '') - if 'SOABI' in confvars: + so_ext = confvars.get("SO", "") + if "SOABI" in confvars: # Does nothing unless SOABI config var exists - so_ext = so_ext.replace('.' + confvars.get('SOABI'), '', 1) + so_ext = so_ext.replace("." + confvars.get("SOABI"), "", 1) return so_ext @@ -69,407 +72,427 @@ def load_library(libname: str) -> CDLL: return cdll.LoadLibrary(libname_ext) -class _Structure_(Structure): # pylint: disable=invalid-name,protected-access - """patch Structure class to allow returning it arguments as dict.""" - def to_dict(self) -> dict: - """return structure as dictionary.""" - return {key: getattr(self, key) for key, _ in self._fields_} - - def to_list(self) -> list: - """return structure as list. Order is the same as in structure.""" - return [getattr(self, key) for key, _ in self._fields_] - - -class SCCOptions(_Structure_): - """Options for evaluating a SCC-Hamiltonian.""" - _fields_ = [ - ('print_level', c_int), - ('parallel', c_int), - ('accuracy', c_double), - ('electronic_temperature', c_double), - ('gradient', c_bool), - ('restart', c_bool), - ('ccm', c_bool), - ('max_iterations', c_int), - ('solvent', c_char*20), - ] - - -class PEEQOptions(_Structure_): - """Options for evaluating a EEQ-Hamiltonian.""" - _fields_ = [ - ('print_level', c_int), - ('parallel', c_int), - ('accuracy', c_double), - ('electronic_temperature', c_double), - ('gradient', c_bool), - ('ccm', c_bool), - ('solvent', c_char*20), - ] - - -def check_ndarray(array: np.ndarray, ctype, size: int, name="array") -> None: - """check if we got the correct array data""" - if not isinstance(array, np.ndarray): - raise ValueError("{} must be of type ndarray".format(name)) - if array.size != size: - raise ValueError("{} does not have the correct size of {}" - .format(name, size)) - if as_ctype is not None: - if as_ctype(array.dtype) != ctype: - raise ValueError("{} must be of {} compatible type" - .format(name, ctype)) +class Param(Enum): + """Possible parametrisations for the Calculator class""" + + GFN2xTB = auto() + GFN1xTB = auto() + GFN0xTB = auto() + GFNFF = auto() + + +VEnvironment = c_void_p +VMolecule = c_void_p +VCalculator = c_void_p +VResults = c_void_p + +_XTB_API_6_3 = { + "xtb_newEnvironment": (VEnvironment, []), + #"xtb_delEnvironment": (None, [POINTER(VEnvironment)]), + "xtb_checkEnvironment": (c_int, [VEnvironment]), + "xtb_showEnvironment": (None, [VEnvironment, c_char_p]), + "xtb_setOutput": (None, [VEnvironment, c_char_p]), + "xtb_releaseOutput": (None, [VEnvironment]), + "xtb_setVerbosity": (None, [VEnvironment, c_int]), + "xtb_newMolecule": ( + VMolecule, + [ + VEnvironment, + POINTER(c_int), + POINTER(c_int), + POINTER(c_double), + POINTER(c_double), + POINTER(c_int), + POINTER(c_double), + POINTER(c_bool), + ], + ), + #"xtb_delMolecule": (None, [POINTER(VMolecule)]), + "xtb_updateMolecule": ( + None, + [VEnvironment, VMolecule, POINTER(c_double), POINTER(c_double),], + ), + "xtb_newCalculator": (VCalculator, []), + #"xtb_delCalculator": (None, [POINTER(VCalculator)]), + "xtb_loadGFN0xTB": (None, [VEnvironment, VMolecule, VCalculator, c_char_p]), + "xtb_loadGFN1xTB": (None, [VEnvironment, VMolecule, VCalculator, c_char_p]), + "xtb_loadGFN2xTB": (None, [VEnvironment, VMolecule, VCalculator, c_char_p]), + "xtb_loadGFNFF": (None, [VEnvironment, VMolecule, VCalculator, c_char_p]), + "xtb_setSolvent": ( + None, + [ + VEnvironment, + VCalculator, + c_char_p, + POINTER(c_int), + POINTER(c_double), + POINTER(c_int), + ], + ), + "xtb_releaseSolvent": (None, [VEnvironment, VCalculator]), + "xtb_setExternalCharges": ( + None, + [ + VEnvironment, + VCalculator, + POINTER(c_int), + POINTER(c_int), + POINTER(c_double), + POINTER(c_double), + ], + ), + "xtb_releaseExternalCharges": (None, [VEnvironment, VCalculator]), + "xtb_singlepoint": (None, [VEnvironment, VMolecule, VCalculator, VResults]), + "xtb_newResults": (VResults, []), + #"xtb_delResults": (None, [POINTER(VResults)]), + "xtb_getEnergy": (None, [VEnvironment, VResults, POINTER(c_double)]), + "xtb_getGradient": (None, [VEnvironment, VResults, POINTER(c_double)]), + "xtb_getVirial": (None, [VEnvironment, VResults, POINTER(c_double)]), + "xtb_getDipole": (None, [VEnvironment, VResults, POINTER(c_double)]), + "xtb_getBondOrders": (None, [VEnvironment, VResults, POINTER(c_double)]), +} class XTBLibrary: - """wrapper for the xtb shared library""" - - # define periodic GFN0-xTB interface - _GFN0_PBC_calculation_ = ( - POINTER(c_int), # number of atoms - POINTER(c_int), # atomic numbers, dimension(number of atoms) - POINTER(c_double), # molecular charge - POINTER(c_int), # number of unpaired electrons - POINTER(c_double), # cartesian coordinates, dimension(3*number of atoms) - POINTER(c_double), # lattice parameters, dimension(9) - POINTER(c_bool), # periodicity of the system - POINTER(PEEQOptions), - c_char_p, # output file name - POINTER(c_double), # energy - POINTER(c_double), # gradient, dimension(3*number of atoms) - POINTER(c_double), # stress tensor, dimension(9) - POINTER(c_double), # lattice gradient, dimension(9) - ) - - # define GFN0-xTB interface - _GFN0_calculation_ = ( - POINTER(c_int), # number of atoms - POINTER(c_int), # atomic numbers, dimension(number of atoms) - POINTER(c_double), # molecular charge - POINTER(c_int), # number of unpaired electrons - POINTER(c_double), # cartesian coordinates, dimension(3*number of atoms) - POINTER(PEEQOptions), - c_char_p, # output file name - POINTER(c_double), # energy - POINTER(c_double), # gradient, dimension(3*number of atoms) - ) - - # define GFN1-xTB interface - _GFN1_calculation_ = ( - POINTER(c_int), # number of atoms - POINTER(c_int), # atomic numbers, dimension(number of atoms) - POINTER(c_double), # molecular charge - POINTER(c_int), # number of unpaired electrons - POINTER(c_double), # cartesian coordinates, dimension(3*number of atoms) - POINTER(SCCOptions), - c_char_p, # output file name - POINTER(c_double), # energy - POINTER(c_double), # gradient, dimension(3*number of atoms) - ) - - # define GFN2-xTB interface - _GFN2_calculation_ = ( - POINTER(c_int), # number of atoms - POINTER(c_int), # atomic numbers, dimension(number of atoms) - POINTER(c_double), # molecular charge - POINTER(c_int), # number of unpaired electrons - POINTER(c_double), # cartesian coordinates, dimension(3*number of atoms) - POINTER(SCCOptions), - c_char_p, # output file name - POINTER(c_double), # energy - POINTER(c_double), # gradient, dimension(3*number of atoms) - POINTER(c_double), - POINTER(c_double), - POINTER(c_double), - POINTER(c_double), - POINTER(c_double), - ) - - _GBSA_model_preload_ = ( - POINTER(c_double), # dielectric data - POINTER(c_double), # molar mass (g/mol) - POINTER(c_double), # solvent density (g/cm^3) - POINTER(c_double), # Born radii - POINTER(c_double), # Atomic surfaces - POINTER(c_double), # Gshift (gsolv=reference vs. gsolv) - POINTER(c_double), # offset parameter (fitted) - POINTER(c_double), - POINTER(c_double), # Surface tension (mN/m=dyn/cm), dimension(94) - POINTER(c_double), # dielectric descreening parameters, dimension(94) - POINTER(c_double), # hydrogen bond strength, dimension(94) - ) - - _GBSA_calculation_ = ( - POINTER(c_int), # number of atoms - POINTER(c_int), # atomic numbers, dimension(number of atoms) - POINTER(c_double), # cartesian coordinates, dimension(3*number of atoms) - c_char_p, # solvent name - POINTER(c_int), # reference state (0: gsolv=reference, 1: gsolv) - POINTER(c_double), # temperature (only for reference=0 or 2 - POINTER(c_int), # method (1 or 2) - POINTER(c_int), # angular grid size (available Lebedev-grids) - c_char_p, # output file name - POINTER(c_double), # Born radii, dimension(number of atoms) - POINTER(c_double), # SASA, dimension(number of atoms) - ) - - _lebedev_grids_ = [6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, - 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, - 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810] + """Shared library instance""" def __init__(self, library: Optional[CDLL] = None): - """construct library from CDLL object.""" if library is not None: - self.library = library + if isinstance(library, XTBLibrary): + self._lib = library._lib + else: + self._lib = library else: - self.library = load_library('libxtb') - - self._set_argtypes_() - - def _set_argtypes_(self) -> None: - """define all interfaces.""" - self.library.GFN0_calculation.argtypes = self._GFN0_calculation_ - self.library.GFN1_calculation.argtypes = self._GFN1_calculation_ - self.library.GFN2_calculation.argtypes = self._GFN2_calculation_ - self.library.GFN0_PBC_calculation.argtypes = self._GFN0_PBC_calculation_ - self.library.GBSA_model_preload.argtypes = self._GBSA_model_preload_ - self.library.GBSA_calculation.argtypes = self._GBSA_calculation_ - - # pylint: disable=invalid-name, too-many-arguments, too-many-locals - def GFN0Calculation(self, natoms: int, numbers, positions, options: dict, - charge: float = 0.0, magnetic_moment: int = 0, - output: str = "-", cell=None, pbc=None) -> dict: - """wrapper for calling the GFN0 Calculator from the library.""" - periodic = cell is not None and pbc is not None - - check_ndarray(numbers, c_int, natoms, "numbers") - check_ndarray(positions, c_double, 3*natoms, "positions") - if periodic: - check_ndarray(cell, c_double, 9, "cell") - check_ndarray(pbc, c_bool, 3, "pbc") - - energy = c_double(0.0) - gradient = np.zeros((natoms, 3), dtype=c_double) - - # turn all strings to binary data, such that ctypes will not complain - l_options = {key: val.encode('utf-8') if isinstance(val, str) else val - for key, val in options.items()} - - if periodic: - cell_gradient = np.zeros((3, 3), dtype=c_double) - stress_tensor = np.zeros((3, 3), dtype=c_double) - args = [ - c_int(natoms), - numbers.ctypes.data_as(POINTER(c_int)), - c_double(charge), - c_int(magnetic_moment), - positions.ctypes.data_as(POINTER(c_double)), - cell.ctypes.data_as(POINTER(c_double)), - pbc.ctypes.data_as(POINTER(c_bool)), - PEEQOptions(**l_options), - output.encode('utf-8'), - energy, - gradient.ctypes.data_as(POINTER(c_double)), - stress_tensor.ctypes.data_as(POINTER(c_double)), - cell_gradient.ctypes.data_as(POINTER(c_double)), - ] - stat = self.library.GFN0_PBC_calculation(*args) + self._lib = load_library("libxtb") + + for name, signature in _XTB_API_6_3.items(): + restype, argtypes = signature + _set(self._lib, name, restype, argtypes) + + +class Environment: + """Calculation environment""" + + _env = None + + def __init__(self, library: Optional[CDLL] = None): + """Create new xtb calculation environment object""" + global _libxtb + if _libxtb is None: + _libxtb = XTBLibrary(library) + self._env = c_void_p(_libxtb._lib.xtb_newEnvironment()) + + def __del__(self): + """Delete a xtb calculation environment object""" + global _libxtb + if self._env is not None: + _libxtb._lib.xtb_delEnvironment(byref(self._env)) + + def check(self) -> int: + """Check current status of calculation environment""" + global _libxtb + return _libxtb._lib.xtb_checkEnvironment(self._env) + + def show(self, message: str) -> None: + """Show and empty error stack""" + global _libxtb + _message = message.encode() + _libxtb._lib.xtb_showEnvironment(self._env, _message) + + def set_output(self, filename: str) -> None: + """Bind output from this environment""" + global _libxtb + _filename = filename.encode() + _libxtb._lib.xtb_setOutput(self._env, _filename) + + def release_output(self) -> None: + """Release output unit from this environment""" + global _libxtb + _libxtb._lib.xtb_releaseOutput(self._env) + + def set_verbosity(self, verbosity: int) -> None: + """Set verbosity of calculation output""" + global _libxtb + _libxtb._lib.xtb_setVerbosity(self._env, verbosity) + + +class Molecule(Environment): + """Molecular structure data""" + + _mol = None + + def __init__( + self, + numbers: List[int], + positions: List[float], + charge: Optional[float] = None, + uhf: Optional[int] = None, + lattice: Optional[List[float]] = None, + periodic: Optional[List[bool]] = None, + library: Optional[CDLL] = None, + ): + """Create new molecular structure data""" + global _libxtb + Environment.__init__(self, library) + if isinstance(positions, np.ndarray): + if positions.size % 3 != 0: + raise ValueError("Expected tripels of cartesian coordinates") else: - args = [ - c_int(natoms), - numbers.ctypes.data_as(POINTER(c_int)), - c_double(charge), - c_int(magnetic_moment), - positions.ctypes.data_as(POINTER(c_double)), - PEEQOptions(**l_options), - output.encode('utf-8'), - energy, - gradient.ctypes.data_as(POINTER(c_double)), - ] - stat = self.library.GFN0_calculation(*args) - - if stat != 0: - raise RuntimeError("GFN0 calculation failed in xtb.") - - results = { - 'energy': energy.value, - 'gradient': gradient, - } - if periodic: - results['cell gradient'] = cell_gradient - results['stress tensor'] = stress_tensor - return results - - # pylint: disable=invalid-name, too-many-arguments, too-many-locals - def GFN1Calculation(self, natoms: int, numbers, positions, options: dict, - charge: float = 0.0, magnetic_moment: int = 0, - output: str = "-") -> dict: - """wrapper for calling the GFN1 Calculator from the library.""" - - check_ndarray(numbers, c_int, natoms, "numbers") - check_ndarray(positions, c_double, 3*natoms, "positions") - - energy = c_double(0.0) - gradient = np.zeros((natoms, 3), dtype=c_double) - dipole = np.zeros(3, dtype=c_double) - charges = np.zeros(natoms, dtype=c_double) - wiberg = np.zeros((natoms, natoms), dtype=c_double) - - # turn all strings to binary data, such that ctypes will not complain - l_options = {key: val.encode('utf-8') if isinstance(val, str) else val - for key, val in options.items()} - - args = [ - c_int(natoms), - numbers.ctypes.data_as(POINTER(c_int)), - c_double(charge), - c_int(magnetic_moment), - positions.ctypes.data_as(POINTER(c_double)), - SCCOptions(**l_options), - output.encode('utf-8'), - energy, - gradient.ctypes.data_as(POINTER(c_double)), - dipole.ctypes.data_as(POINTER(c_double)), - charges.ctypes.data_as(POINTER(c_double)), - wiberg.ctypes.data_as(POINTER(c_double)), - ] - stat = self.library.GFN1_calculation(*args) - - if stat != 0: - raise RuntimeError("GFN1 calculation failed in xtb.") - - return { - 'energy': energy.value, - 'gradient': gradient, - 'dipole moment': dipole, - 'charges': charges, - 'wiberg': wiberg, - } + if len(positions) % 3 != 0: + raise ValueError("Expected tripels of cartesian coordinates") - # pylint: disable=invalid-name, too-many-arguments, too-many-locals - def GFN2Calculation(self, natoms: int, numbers, positions, options: dict, - charge: float = 0.0, magnetic_moment: int = 0, - output: str = "-") -> dict: - """wrapper for calling the GFN2 Calculator from the library.""" - - check_ndarray(numbers, c_int, natoms, "numbers") - check_ndarray(positions, c_double, 3*natoms, "positions") - - energy = c_double(0.0) - gradient = np.zeros((natoms, 3), dtype=c_double) - dipole = np.zeros(3, dtype=c_double) - charges = np.zeros(natoms, dtype=c_double) - dipoles = np.zeros((natoms, 3), dtype=c_double) - quadrupoles = np.zeros((natoms, 6), dtype=c_double) - wiberg = np.zeros((natoms, natoms), dtype=c_double) - - # turn all strings to binary data, such that ctypes will not complain - l_options = {key: val.encode('utf-8') if isinstance(val, str) else val - for key, val in options.items()} - - args = [ - c_int(natoms), - numbers.ctypes.data_as(POINTER(c_int)), - c_double(charge), - c_int(magnetic_moment), - positions.ctypes.data_as(POINTER(c_double)), - SCCOptions(**l_options), - output.encode('utf-8'), - energy, - gradient.ctypes.data_as(POINTER(c_double)), - dipole.ctypes.data_as(POINTER(c_double)), - charges.ctypes.data_as(POINTER(c_double)), - dipoles.ctypes.data_as(POINTER(c_double)), - quadrupoles.ctypes.data_as(POINTER(c_double)), - wiberg.ctypes.data_as(POINTER(c_double)), - ] - stat = self.library.GFN2_calculation(*args) - - if stat != 0: - raise RuntimeError("GFN2 calculation failed in xtb.") - - return { - 'energy': energy.value, - 'gradient': gradient, - 'dipole moment': dipole, - 'charges': charges, - 'dipoles': dipoles, - 'quadrupoles': quadrupoles, - 'wiberg': wiberg, - } + if isinstance(positions, np.ndarray): + if 3 * len(numbers) != positions.size: + raise ValueError("Dimension missmatch between numbers and postions") + else: + if 3 * len(numbers) != len(positions): + raise ValueError("Dimension missmatch between numbers and postions") + + self._natoms = len(numbers) + _numbers = np.array(numbers, dtype="i4") + _positions = np.array(positions, dtype="float") + + if lattice is not None: + if len(lattice) != 9: + raise ValueError("Invalid lattice provided") + _lattice = np.array(lattice, dtype="float") + else: + _lattice = None - # pylint: disable=invalid-name, too-many-arguments, too-many-locals - def GBSACalculation(self, natoms: int, numbers, positions, - options: dict, output: str) -> dict: - """wrapper for calling the GBSA Calculator from the library.""" - - check_ndarray(numbers, c_int, natoms, "numbers") - check_ndarray(positions, c_double, 3*natoms, "positions") - - if 'solvent' not in options: - raise ValueError("solvent has to present in options") - ref_state = options['reference'] if 'reference' in options else 0 - if ref_state not in [0, 1, 2]: - raise ValueError("reference state {} is no defined".format(ref_state)) - temp = options['temperature'] if 'temperature' in options else 298.15 - if temp <= 0.0: - raise ValueError("negative temperature does not make sense") - grid = options['grid size'] if 'grid size' in options else 230 - if grid not in self._lebedev_grids_: - raise ValueError("angular grid {} is no Lebedev grid".format(grid)) - method = options['method'] if 'method' in options else 2 - if method not in [1, 2]: - raise ValueError("method {} is no available".format(method)) - - born = np.zeros(natoms, dtype=c_double) - sasa = np.zeros(natoms, dtype=c_double) - - args = ( - c_int(natoms), - numbers.ctypes.data_as(POINTER(c_int)), - positions.ctypes.data_as(POINTER(c_double)), - options['solvent'].encode('utf-8'), - c_int(ref_state), - c_double(temp), - c_int(method), - c_int(grid), - output.encode('utf-8'), - born.ctypes.data_as(POINTER(c_double)), - sasa.ctypes.data_as(POINTER(c_double)), + if periodic is not None: + if len(periodic) != 3: + raise ValueError("Invalid periodicity provided") + _periodic = np.array(periodic, dtype="bool") + else: + _periodic = None + + self._mol = c_void_p(_libxtb._lib.xtb_newMolecule( + self._env, + c_int(self._natoms), + _cast(POINTER(c_int), _numbers), + _cast(POINTER(c_double), _positions), + None if charge is None else c_double(charge), + None if uhf is None else c_int(uhf), + _cast(POINTER(c_double), _lattice), + _cast(POINTER(c_bool), _periodic), + )) + + if self.check() != 0: + raise ValueError("Could not initialize molecular structure data") + + def __del__(self): + """Delete molecular structure data""" + global _libxtb + Environment.__del__(self) + if self._mol is not None: + _libxtb._lib.xtb_delMolecule(byref(self._mol)) + + def __len__(self): + return self._natoms + + def update( + self, positions: List[float], lattice: Optional[List[float]] = None, + ): + """Update coordinates and lattice parameters""" + global _libxtb + + if 3 * len(self) != len(positions): + raise ValueError("Dimension missmatch for postions") + _positions = np.array(positions, dtype="float") + + if lattice is not None: + if len(lattice) != 9: + raise ValueError("Invalid lattice provided") + _lattice = np.array(lattice, dtype="float") + else: + _lattice = None + + _libxtb._lib.xtb_updateMolecule( + self._env, + self._mol, + _cast(POINTER(c_double), _positions), + _cast(POINTER(c_double), _lattice), ) - stat = self.library.GBSA_calculation(*args) - if stat != 0: - raise RuntimeError("GBSA calculation failed in xtb.") - return { - 'born': born, - 'sasa': sasa, + if self.check() != 0: + raise ValueError("Could not update molecular structure data") + + +class Results(Environment): + """Calculation results""" + + _res = None + + def __init__(self, mol: Molecule, library: Optional[CDLL] = None): + """Create new singlepoint results object""" + global _libxtb + Environment.__init__(self, library) + self._res = c_void_p(_libxtb._lib.xtb_newResults()) + self._natoms = len(mol) + + def __del__(self): + """Delete singlepoint results object""" + global _libxtb + Environment.__del__(self) + if self._res is not None: + _libxtb._lib.xtb_delResults(byref(self._res)) + + def __len__(self): + return self._natoms + + def get_energy(self): + """Query singlepoint results object for energy""" + global _libxtb + _energy = c_double(0.0) + _libxtb._lib.xtb_getEnergy(self._env, self._res, _energy) + if self.check() != 0: + raise ValueError("Energy is not available") + return _energy.value + + def get_gradient(self): + """Query singlepoint results object for gradient""" + global _libxtb + _gradient = np.zeros((len(self), 3)) + _libxtb._lib.xtb_getGradient( + self._env, self._res, _cast(POINTER(c_double), _gradient) + ) + if self.check() != 0: + raise ValueError("Gradient is not available") + return _gradient + + def get_virial(self): + """Query singlepoint results object for virial""" + global _libxtb + _virial = np.zeros((3, 3)) + _libxtb._lib.xtb_getVirial(self._env, self._res, _cast(POINTER(c_double), _virial)) + if self.check() != 0: + raise ValueError("Virial is not available") + return _virial + + def get_dipole(self): + """Query singlepoint results object for dipole""" + global _libxtb + _dipole = np.zeros(3) + _libxtb._lib.xtb_getDipole(self._env, self._res, _cast(POINTER(c_double), _dipole)) + if self.check() != 0: + raise ValueError("Dipole is not available") + return _dipole + + def get_charges(self): + """Query singlepoint results object for partial charges""" + global _libxtb + _charges = np.zeros(len(self)) + _libxtb._lib.xtb_getCharges( + self._env, self._res, _cast(POINTER(c_double), _charges) + ) + if self.check() != 0: + raise ValueError("Charges are not available") + return _charges + + def get_bond_orders(self): + """Query singlepoint results object for bond orders""" + global _libxtb + _bond_orders = np.zeros((len(self), len(self))) + _libxtb._lib.xtb_getBondOrders( + self._env, self._res, _cast(POINTER(c_double), _bond_orders) + ) + if self.check() != 0: + raise ValueError("Bond orders are not available") + return _bond_orders + + +class Calculator(Molecule): + """Singlepoint calculator""" + + _calc = None + + def __init__( + self, + param: Param, + numbers: List[int], + positions: List[float], + charge: Optional[float] = None, + uhf: Optional[int] = None, + lattice: Optional[List[float]] = None, + periodic: Optional[List[bool]] = None, + solvent: Optional[str] = None, + library: Optional[CDLL] = None, + ): + """Create new calculator object""" + global _libxtb + Molecule.__init__( + self, numbers, positions, charge, uhf, lattice, periodic, library, + ) + + self._loader = { + Param.GFN2xTB: _libxtb._lib.xtb_loadGFN2xTB, + Param.GFN1xTB: _libxtb._lib.xtb_loadGFN1xTB, + Param.GFN0xTB: _libxtb._lib.xtb_loadGFN0xTB, + Param.GFNFF: _libxtb._lib.xtb_loadGFNFF, } - def GBSA_model_preload(self, epsv: float, smass: float, rhos: float, c1: float, - rprobe: float, gshift: float, soset: float, - gamscale, sx, tmp) -> None: - """preload parameters into GBSA""" - - check_ndarray(gamscale, c_double, 94, "gamscale") - check_ndarray(sx, c_double, 94, "sx") - check_ndarray(tmp, c_double, 94, "tmp") - - args = ( - c_double(epsv), - c_double(smass), - c_double(rhos), - c_double(c1), - c_double(rprobe), - c_double(gshift), - c_double(soset), - None, - gamscale.ctypes.data_as(POINTER(c_double)), - sx.ctypes.data_as(POINTER(c_double)), - tmp.ctypes.data_as(POINTER(c_double)), + self._calc = c_void_p(_libxtb._lib.xtb_newCalculator()) + self._load(param) + if solvent is not None: + self.set_solvent(solvent) + + def __del__(self): + """Delete calculator object""" + global _libxtb + Molecule.__del__(self) + if self._calc is not None: + _libxtb._lib.xtb_delCalculator(byref(self._calc)) + + def _load(self, param: Param): + """Load parametrisation data into calculator""" + + self._loader[param]( + self._env, self._mol, self._calc, None, ) - stat = self.library.GBSA_model_preload(*args) + if self.check() != 0: + raise ValueError("Could not load parametrisation data") + + def set_solvent(self, solvent: Optional[str]): + """Add/Remove solvation model to/from calculator""" + global _libxtb + + if solvent is None: + _libxtb._lib.xtb_releaseSolvent( + self._env, self._calc, + ) + else: + _solvent = solvent.encode() + + _libxtb._lib.xtb_setSolvent( + self._env, self._calc, _solvent, None, None, None, + ) + + if self.check() != 0: + raise ValueError("Could not load parametrisation data") + + def singlepoint(self, res: Optional[Results] = None) -> Results: + """Perform singlepoint calculation""" + global _libxtb + + _res = Results(self, _libxtb._lib) if res is None else res + _libxtb._lib.xtb_singlepoint( + self._env, self._mol, self._calc, _res._res, + ) + + if self.check() != 0: + raise ValueError("Single point calculation failed") + + return _res + + +def _cast(ctype, array): + """Cast a numpy array to a FFI pointer""" + return None if array is None else array.ctypes.data_as(ctype) + - if stat != 0: - raise RuntimeError("GBSA parameters could not be loaded by xtb.") +def _set(lib, name, restype, argtypes): + """Setup the ctypes signature for a function in a shared library""" + func = lib.__getattr__(name) + func.restype = restype + func.argtypes = argtypes diff --git a/python/xtb/solvation.py b/python/xtb/solvation.py deleted file mode 100644 index 1ca46412a..000000000 --- a/python/xtb/solvation.py +++ /dev/null @@ -1,133 +0,0 @@ -# This file is part of xtb. -# -# Copyright (C) 2019-2020 Sebastian Ehlert -# -# xtb 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. -# -# xtb 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 xtb. If not, see . -"""All solvation related implementations for xtb.""" - -from __future__ import print_function -from typing import List - -from ctypes import c_int, c_double - -from ase.calculators.calculator import Calculator, all_changes -from ase.units import Bohr - -import numpy as np - -from xtb.interface import XTBLibrary -from xtb.calculators import XTB, GFN0, GFN1, GFN2 - -__all__ = ['GBSA', 'GSOLV', 'GSOLV_REFERENCE', 'MOL1BAR'] - -GSOLV = 0 -GSOLV_REFERENCE = 1 -MOL1BAR = 2 - - -class GBSA(Calculator): - """Interface to the GBSA method""" - implemented_properties = [ - 'energy', - 'born radii', - 'sasa', - ] - default_parameters = { - 'solvent': None, - 'temperature': 298.15, - 'reference': GSOLV, - 'grid size': 230, - 'method': 2, - } - calc = None - _debug = False - - # pylint: disable=too-many-arguments - def __init__(self, restart=None, ignore_bad_restart_file=False, - label='gbsa', atoms=None, library=None, calc=None, **kwargs): - """Construct the XTB base calculator object.""" - - Calculator.__init__(self, restart, ignore_bad_restart_file, label, atoms, - **kwargs) - - if calc is not None: - if not isinstance(calc, XTB): - raise ValueError("Only xTB-calculators are supported in GBSA.") - self.calc = calc - # borrow library from xTB calculator - self.library = self.calc.library - else: - if library is not None: - self.library = library - else: - self.library = XTBLibrary() - - # loads the default parameters and updates with actual values - self.parameters = self.get_default_parameters() - # if we have a calculator we adjust the method keyword - if isinstance(calc, (GFN2, GFN0)): - self.set(method=2) # GFN0 is treated as GFN2 (no parameters yet) - if isinstance(calc, GFN1): - self.set(method=1) - # now set all parameters - self.set(**kwargs) - - if self.parameters['solvent'] is None: - raise ValueError("Solvent keyword is missing, cannot setup GBSA.") - - def output_file_name(self) -> str: - """create output file name from label as ctype""" - if self.label is not None: - return self.label + ".out" - return "-" # standard output - - def create_arguments(self) -> dict: - """create a list of arguments.""" - if any(self.atoms.pbc): - raise NotImplementedError("GBSA is not available with PBC!") - kwargs = { - 'natoms': len(self.atoms), - 'numbers': np.array(self.atoms.get_atomic_numbers(), dtype=c_int), - 'positions': np.array(self.atoms.get_positions()/Bohr, dtype=c_double), - 'options': self.parameters, - 'output': self.output_file_name(), - } - return kwargs - - # pylint: disable=dangerous-default-value - def calculate(self, atoms=None, properties: List[str] = None, - system_changes: List[str] = all_changes) -> None: - """perform calculation with libxtb""" - - if not properties: - properties = ['energy'] - Calculator.calculate(self, atoms, properties, system_changes) - - kwargs = self.create_arguments() - results = self.library.GBSACalculation(**kwargs) - self.results['born radii'] = results['born']*Bohr - self.results['sasa'] = results['sasa']*Bohr**2 - - if self.calc is not None: - # evaluate system in gas phase - self.calc.set(solvent="none") - self.atoms.set_calculator(self.calc) - vacuum = self.atoms.get_potential_energy() - # evaluate system in solution - self.calc.reset() - self.calc.set(solvent=self.parameters['solvent']) - self.atoms.set_calculator(self.calc) - solution = self.atoms.get_potential_energy() - # derive properties like solvation free energy - self.results['energy'] = solution - vacuum diff --git a/python/xtb/test/test_calculators.py b/python/xtb/test/test_calculators.py deleted file mode 100644 index 2af750d85..000000000 --- a/python/xtb/test/test_calculators.py +++ /dev/null @@ -1,117 +0,0 @@ -# This file is part of xtb. -# -# Copyright (C) 2019-2020 Sebastian Ehlert -# -# xtb 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. -# -# xtb 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 xtb. If not, see . - -"""Tests for the ASE interface to xtb.""" - - -class Testing: - """Unit tests for xTB methods.""" - - def test_g2_gfn0(self): # pylint: disable=no-self-use - """Short test for GFN0-xTB on a subset of the G2.""" - from pytest import approx - from ase.collections import g2 - from xtb.calculators import GFN0 as GFN - - thr = 1.0e-7 - subset = { - "cyclobutene": (-318.7866120321342, 0.18172370289855405), - "CH3ONO": (-358.73445196501973, 0.409421901652166), - "SiH3": (-93.9024809777723, 0.27588327621925385), - "C3H6_D3h": (-262.5099874403974, 0.10397633188045725), - "CO2": (-236.48715959452625, 0.7513923402713432), - "NO": (-165.74409502219441, 0.06375042027444998), - "SiO": (-136.30493006096225, 0.41855575434941694), - "C3H4_D2d": (-231.20208222258086, 0.11557009795103503), - "COF2": (-357.1115761281822, 0.45554166948678626), - "2-butyne": (-319.609294498653, 0.1654237936464811), - "C2H5": (-188.90399485898985, 0.08291942020794801), - "BF3": (-356.10418929536525, 0.15840229870826247), - "N2O": (-246.88770058584996, 1.5739076279140762), - } - - calc = GFN(accuracy=0.01) - for name, data in subset.items(): - atoms = g2[name] - energy, gnorm = data - atoms.set_calculator(calc) - assert energy == approx(atoms.get_potential_energy(), thr) - assert gnorm == approx(abs(atoms.get_forces().flatten()).mean(), thr) - - def test_g2_gfn1(self): # pylint: disable=no-self-use - """Short test for GFN1-xTB on a subset of the G2.""" - from pytest import approx - from ase.collections import g2 - from xtb.calculators import GFN1 as GFN - - thr = 1.0e-7 - subset = { - "C4H4NH": (-389.83886433843026, 0.11845899422716655), - "CH3SCH3": (-303.42619191556383, 0.06645628082687234), - "SiH2_s3B1d": (-74.46528605312987, 0.08310436985665531), - "CH3CO": (-284.3582641984081, 0.20638671880968007), - "CO": (-183.19840869756732, 0.9977431424545925), - "H2CO": (-213.4823225876121, 0.39657080790031), - "CH3COOH": (-429.92204679760215, 0.24893030739923092), - "HCF3": (-492.1472682009942, 0.2744291603498096), - "CS2": (-257.45587723444686, 0.1020065624570768), - "SiH2_s1A1d": (-77.00356692788125, 0.08145702394953913), - "C4H4S": (-388.5560517704511, 0.16524901027806613), - } - - calc = GFN(accuracy=0.01) - for name, data in subset.items(): - atoms = g2[name] - energy, gnorm = data - atoms.set_calculator(calc) - assert energy == approx(atoms.get_potential_energy(), thr) - assert gnorm == approx(abs(atoms.get_forces().flatten()).mean(), thr) - - def test_g2_gfn2(self): # pylint: disable=no-self-use - """Short test for GFN2-xTB on a subset of the G2.""" - from pytest import approx - from ase.collections import g2 - from xtb.calculators import GFN2 as GFN - - thr = 1.0e-7 - subset = { - "F2O": (-362.11924747368295, 0.7143257863045354), - "SO2": (-309.1999361327123, 0.8740554743818797), - "H2CCl2": (-333.2265933353902, 0.046195744939883536), - "CF3CN": (-580.4856945730149, 0.6213782830693683), - "HCN": (-149.6789224980101, 1.015081156428575), - "C2H6NH": (-292.6330765170243, 0.07259870951395032), - "OCS": (-257.1846815478052, 0.549785128276834), - "ClO": (-230.53675733271504, 0.6103015273199656), - "C3H8": (-285.74053465625167, 0.08676710697381718), - "HF": (-142.12158161038408, 0.024290012553866525), - "O2": (-215.0347545913049, 0.8388557044315267), - "SO": (-196.79890058895995, 0.8295962355319549), - "NH": (-87.11950843824161, 0.18471418108510287), - "C2F4": (-630.1300665999896, 0.11452334981560626), - "NF3": (-461.2701116859826, 0.4165894734601184), - "CH2_s3B1d": (-79.94618154681717, 0.1971528027462761), - "CH3CH2Cl": (-309.61292017128613, 0.05533781132943324), - } - - calc = GFN(accuracy=0.01) - for name, data in subset.items(): - atoms = g2[name] - energy, gnorm = data - atoms.set_calculator(calc) - assert energy == approx(atoms.get_potential_energy(), thr) - assert gnorm == approx(abs(atoms.get_forces().flatten()).mean(), thr) diff --git a/python/xtb/test/test_interface.py b/python/xtb/test/test_interface.py index 86a382365..ad36910fa 100644 --- a/python/xtb/test/test_interface.py +++ b/python/xtb/test/test_interface.py @@ -25,31 +25,20 @@ def test_library(): # check if we can load this one lib = load_library("libxtb") # check if we find some functions - assert lib.GFN0_calculation is not None - assert lib.GFN1_calculation is not None - assert lib.GFN2_calculation is not None + assert lib.xtb_singlepoint is not None from xtb.interface import XTBLibrary libxtb = XTBLibrary(library=lib) - assert libxtb.library is lib - # check if we find some functions - assert hasattr(libxtb, "GFN0Calculation") - assert hasattr(libxtb, "GFN1Calculation") - assert hasattr(libxtb, "GFN2Calculation") - - libxtb = XTBLibrary() - # check if we find some functions - assert hasattr(libxtb, "GFN0Calculation") - assert hasattr(libxtb, "GFN1Calculation") - assert hasattr(libxtb, "GFN2Calculation") + assert libxtb._lib is lib class Testing: """test all defined interfaces""" + from ctypes import c_int, c_double from numpy import array - from xtb.interface import XTBLibrary + from xtb.interface import Calculator, Results, Param natoms = 24 numbers = array( @@ -85,35 +74,16 @@ class Testing: ], dtype=c_double, ) - lib = XTBLibrary() def test_gfn0_interface(self): """check if the GFN0-xTB interface is working correctly.""" thr = 1.0e-8 from pytest import approx - options = { - "print_level": 1, - "parallel": 0, - "accuracy": 1.0, - "electronic_temperature": 300.0, - "gradient": True, - "ccm": True, - "solvent": "none", - } - kwargs = { - "natoms": self.natoms, - "numbers": self.numbers, - "charge": 0.0, - "positions": self.positions, - "options": options, - "output": "-", - "magnetic_moment": 0, - } - results = self.lib.GFN0Calculation(**kwargs) - assert results - gnorm = abs(results["gradient"].flatten()).mean() - assert approx(results["energy"], thr) == -40.90885036015837 + calc = self.Calculator(self.Param.GFN0xTB, self.numbers, self.positions) + results = calc.singlepoint() + gnorm = abs(results.get_gradient().flatten()).mean() + assert approx(results.get_energy(), thr) == -40.90885036015837 assert approx(gnorm, thr) == 0.00510542946913145 def test_gfn1_interface(self): @@ -121,22 +91,10 @@ def test_gfn1_interface(self): thr = 1.0e-8 from pytest import approx - options = { - "print_level": 1, - "parallel": 0, - "accuracy": 1.0, - "electronic_temperature": 300.0, - "gradient": True, - "restart": False, - "ccm": True, - "max_iterations": 30, - "solvent": "none", - } - args = (self.natoms, self.numbers, self.positions, options, 0.0, 0, "-") - results = self.lib.GFN1Calculation(*args) - assert results - gnorm = abs(results["gradient"].flatten()).mean() - assert approx(results["energy"], thr) == -44.50970242016477 + calc = self.Calculator(self.Param.GFN1xTB, self.numbers, self.positions) + results = calc.singlepoint() + gnorm = abs(results.get_gradient().flatten()).mean() + assert approx(results.get_energy(), thr) == -44.50970242016477 assert approx(gnorm, thr) == 0.004842112340521374 def test_gfn2_gbsa_interface(self): @@ -144,28 +102,11 @@ def test_gfn2_gbsa_interface(self): thr = 1.0e-8 from pytest import approx - options = { - "print_level": 1, - "parallel": 0, - "accuracy": 1.0, - "electronic_temperature": 300.0, - "gradient": True, - "restart": False, - "ccm": True, - "max_iterations": 30, - "solvent": "ch2cl2", - } - results = self.lib.GFN2Calculation( - natoms=self.natoms, - numbers=self.numbers, - charge=0.0, - positions=self.positions, - options=options, - output="-", - magnetic_moment=0, + calc = self.Calculator( + self.Param.GFN2xTB, self.numbers, self.positions, solvent="ch2cl2" ) - assert results - gnorm = abs(results["gradient"].flatten()).mean() - assert approx(results["energy"], thr) == -42.170584528028506 + results = calc.singlepoint() + gnorm = abs(results.get_gradient().flatten()).mean() + assert approx(results.get_energy(), thr) == -42.170584528028506 assert approx(gnorm, thr) == 0.004685246019376688 - assert approx(results["charges"].sum(), thr) == 0.0 + assert approx(results.get_charges().sum(), thr) == 0.0 diff --git a/python/xtb/test/test_solvation.py b/python/xtb/test/test_solvation.py deleted file mode 100644 index 6db1a7c9a..000000000 --- a/python/xtb/test/test_solvation.py +++ /dev/null @@ -1,125 +0,0 @@ -# This file is part of xtb. -# -# Copyright (C) 2019-2020 Sebastian Ehlert -# -# xtb 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. -# -# xtb 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 xtb. If not, see . -"""Tests for the GBSA API.""" - - -class Testing: - from ase import Atoms - - octane = Atoms( - "C8H18", - [ - [0.00000000, 0.00000000, -0.78216596], - [0.00000000, 0.00000000, 0.78216596], - [1.13685895, -0.86573057, -1.33679656], - [0.18131519, 1.41741401, -1.33679656], - [-1.31817414, -0.55168345, -1.33679656], - [1.13685895, 0.86573057, 1.33679656], - [0.18131519, -1.41741401, 1.33679656], - [-1.31817414, 0.55168345, 1.33679656], - [2.10137753, -0.60632401, -0.89602505], - [0.95640314, -1.92812772, -1.16563611], - [1.21642453, -0.71505916, -2.41638078], - [1.19160602, 1.79233328, -1.16563611], - [0.01104713, 1.41098413, -2.41638078], - [-0.52559677, 2.12300833, -0.89602505], - [-1.57578076, -1.51668432, -0.89602505], - [-2.14800916, 0.13579445, -1.16563611], - [-1.22747167, -0.69592497, -2.41638078], - [2.10137753, 0.60632401, 0.89602505], - [0.95640314, 1.92812772, 1.16563611], - [1.21642453, 0.71505916, 2.41638078], - [1.19160602, -1.79233328, 1.16563611], - [0.01104713, -1.41098413, 2.41638078], - [-0.52559677, -2.12300833, 0.89602505], - [-1.57578076, 1.51668432, 0.89602505], - [-2.14800916, -0.13579445, 1.16563611], - [-1.22747167, 0.69592497, 2.41638078], - ], - ) - - def test_gfn1_gbsa_water(self): - from pytest import approx - from xtb.solvation import GBSA - from xtb.calculators import GFN1 - - born = [ - 3.14547888, - 3.14547888, - 2.64032576, - 2.64032576, - 2.64032576, - 2.64032576, - 2.64032576, - 2.64032576, - 2.0475772, - 2.04527513, - 2.02622599, - 2.04527513, - 2.02622599, - 2.0475772, - 2.0475772, - 2.04527513, - 2.02622599, - 2.0475772, - 2.04527513, - 2.02622599, - 2.04527513, - 2.02622599, - 2.0475772, - 2.0475772, - 2.04527513, - 2.02622599, - ] - sasa = [ - -0.0, - -0.0, - 1.81990386, - 1.83838909, - 1.81669523, - 1.81990386, - 1.83838909, - 1.81669523, - 11.40692384, - 12.57434505, - 14.97110513, - 12.56651178, - 14.98457022, - 11.40554208, - 11.41537855, - 12.52852953, - 14.95384903, - 11.40692384, - 12.57434505, - 14.97110513, - 12.56651178, - 14.98457022, - 11.40554208, - 11.41537855, - 12.52852953, - 14.95384903, - ] - - mol = self.octane.copy() - gbsa = GBSA(solvent="water", calc=GFN1()) - - mol.set_calculator(gbsa) - gsolv = mol.get_potential_energy() - - assert approx(gsolv) == 0.048946043429737074 - assert approx(gbsa.get_property("born radii")) == born - assert approx(gbsa.get_property("sasa")) == sasa diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 186e8cbca..9deca33fd 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -43,7 +43,6 @@ list(APPEND srcs "${dir}/bias_path.f90" "${dir}/blowsy.f90" "${dir}/broyden.f90" - "${dir}/calculator.f90" "${dir}/charge_model.f90" "${dir}/coffee.f90" "${dir}/constrain_param.f90" @@ -75,10 +74,7 @@ list(APPEND srcs "${dir}/geosum.f90" "${dir}/getname.f90" "${dir}/getsymnum.f90" - "${dir}/gfn0_calculator.f90" "${dir}/gfn0param.f90" - "${dir}/gfn1_calculator.f90" - "${dir}/gfn2_calculator.f90" "${dir}/gfn_paramset.f90" "${dir}/gfn_prparam.f90" "${dir}/grad_core.f90" diff --git a/src/api/CMakeLists.txt b/src/api/CMakeLists.txt index 4ff4597fb..4842d3a77 100644 --- a/src/api/CMakeLists.txt +++ b/src/api/CMakeLists.txt @@ -18,9 +18,11 @@ set(dir "${CMAKE_CURRENT_SOURCE_DIR}") list(APPEND srcs + "${dir}/calculator.f90" + "${dir}/environment.f90" + "${dir}/molecule.f90" "${dir}/interface.f90" - "${dir}/preload.f90" - "${dir}/structs.f90" + "${dir}/results.f90" "${dir}/utils.f90" ) diff --git a/src/api/calculator.f90 b/src/api/calculator.f90 new file mode 100644 index 000000000..398b1a07c --- /dev/null +++ b/src/api/calculator.f90 @@ -0,0 +1,511 @@ +! This file is part of xtb. +! +! Copyright (C) 2019-2020 Sebastian Ehlert +! +! xtb 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. +! +! xtb 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 xtb. If not, see . + +!> API for dealing with the single point calculator +module xtb_api_calculator + use, intrinsic :: iso_c_binding + use xtb_mctc_accuracy, only : wp + use xtb_mctc_io + use xtb_mctc_systools, only : rdpath + use xtb_api_environment + use xtb_api_molecule + use xtb_api_utils + use xtb_gfnff_calculator + use xtb_main_setup + use xtb_solv_gbobc, only : lgbsa + use xtb_type_environment + use xtb_type_molecule + use xtb_type_calculator + use xtb_xtb_calculator + implicit none + private + + public :: VCalculator + public :: newCalculator_api, delCalculator_api + public :: loadGFNFF_api, loadGFN0xTB_api, loadGFN1xTB_api, loadGFN2xTB_api + public :: setSolvent_api, releaseSolvent_api + public :: setExternalCharges_api, releaseExternalCharges_api + + + !> Void pointer to single point calculator + type :: VCalculator + class(TCalculator), allocatable :: ptr + end type VCalculator + + +contains + + +function newCalculator_api() result(vcalc) & + & bind(C, name="xtb_newCalculator") + type(VCalculator), pointer :: calc + type(c_ptr) :: vcalc + + call checkGlobalEnv + + allocate(calc) + vcalc = c_loc(calc) + +end function newCalculator_api + + +subroutine delCalculator_api(vcalc) & + & bind(C, name="xtb_delCalculator") + type(c_ptr), intent(inout) :: vcalc + type(VCalculator), pointer :: calc + + call checkGlobalEnv + + if (c_associated(vcalc)) then + call c_f_pointer(vcalc, calc) + deallocate(calc) + vcalc = c_null_ptr + end if + +end subroutine delCalculator_api + + +subroutine loadGFNFF_api(venv, vmol, vcalc, charptr) & + & bind(C, name="xtb_loadGFNFF") + character(len=*), parameter :: source = 'xtb_api_loadGFNFF' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vmol + type(VMolecule), pointer :: mol + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + character(kind=c_char), intent(in), optional :: charptr(*) + character(len=:, kind=c_char), allocatable :: dummy, filename + character(len=*), parameter :: pFilename = '.param_gfnff.xtb' + type(TGFFCalculator), allocatable :: gff + logical :: exist, exitRun + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vmol)) then + call env%ptr%error("Molecular structure data is not allocated", source) + return + end if + call c_f_pointer(vmol, mol) + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + if (present(charptr)) then + call c_f_character(charptr, dummy) + else + dummy = pFilename + end if + inquire(file=dummy, exist=exist) + if (.not.exist) then + call rdpath(env%ptr%xtbpath, dummy, filename, exist) + if (.not.exist) then + filename = dummy + end if + else + filename = dummy + end if + + allocate(gff) + call newGFFCalculator(env%ptr, mol%ptr, gff, filename, .false.) + + call env%ptr%check(exitRun) + if (exitRun) then + call env%ptr%error("Could not construct GFN-FF calculator", source) + return + end if + + call move_alloc(gff, calc%ptr) + + end if + +end subroutine loadGFNFF_api + + +subroutine loadGFN0xTB_api(venv, vmol, vcalc, charptr) & + & bind(C, name="xtb_loadGFN0xTB") + character(len=*), parameter :: source = 'xtb_api_loadGFN0xTB' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vmol + type(VMolecule), pointer :: mol + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + character(kind=c_char), intent(in), optional :: charptr(*) + character(len=:, kind=c_char), allocatable :: dummy, filename + character(len=*), parameter :: pFilename = 'param_gfn0-xtb.txt' + type(TxTBCalculator), allocatable :: xtb + logical :: exist, exitRun + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vmol)) then + call env%ptr%error("Molecular structure data is not allocated", source) + return + end if + call c_f_pointer(vmol, mol) + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + if (present(charptr)) then + call c_f_character(charptr, dummy) + else + dummy = pFilename + end if + inquire(file=dummy, exist=exist) + if (.not.exist) then + call rdpath(env%ptr%xtbpath, dummy, filename, exist) + if (.not.exist) then + filename = dummy + end if + end if + + allocate(xtb) + call newXTBCalculator(env%ptr, mol%ptr, xtb, filename) + + call env%ptr%check(exitRun) + if (exitRun) then + call env%ptr%error("Could not construct GFN0-xTB calculator", source) + return + end if + + call move_alloc(xtb, calc%ptr) + + end if + +end subroutine loadGFN0xTB_api + + +subroutine loadGFN1xTB_api(venv, vmol, vcalc, charptr) & + & bind(C, name="xtb_loadGFN1xTB") + character(len=*), parameter :: source = 'xtb_api_loadGFN1xTB' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vmol + type(VMolecule), pointer :: mol + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + character(kind=c_char), intent(in), optional :: charptr(*) + character(len=:, kind=c_char), allocatable :: dummy, filename + character(len=*), parameter :: pFilename = 'param_gfn1-xtb.txt' + type(TxTBCalculator), allocatable :: xtb + logical :: exist, exitRun + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vmol)) then + call env%ptr%error("Molecular structure data is not allocated", source) + return + end if + call c_f_pointer(vmol, mol) + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + if (present(charptr)) then + call c_f_character(charptr, dummy) + else + dummy = pFilename + end if + inquire(file=dummy, exist=exist) + if (.not.exist) then + call rdpath(env%ptr%xtbpath, dummy, filename, exist) + if (.not.exist) then + filename = dummy + end if + end if + + allocate(xtb) + call newXTBCalculator(env%ptr, mol%ptr, xtb, filename) + + call env%ptr%check(exitRun) + if (exitRun) then + call env%ptr%error("Could not construct GFN1-xTB calculator", source) + return + end if + + call move_alloc(xtb, calc%ptr) + + end if + +end subroutine loadGFN1xTB_api + + +subroutine loadGFN2xTB_api(venv, vmol, vcalc, charptr) & + & bind(C, name="xtb_loadGFN2xTB") + character(len=*), parameter :: source = 'xtb_api_loadGFN2xTB' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vmol + type(VMolecule), pointer :: mol + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + character(kind=c_char), intent(in), optional :: charptr(*) + character(len=:, kind=c_char), allocatable :: dummy, filename + character(len=*), parameter :: pFilename = 'param_gfn2-xtb.txt' + type(TxTBCalculator), allocatable :: xtb + logical :: exist, exitRun + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vmol)) then + call env%ptr%error("Molecular structure data is not allocated", source) + return + end if + call c_f_pointer(vmol, mol) + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + if (present(charptr)) then + call c_f_character(charptr, dummy) + else + dummy = pFilename + end if + inquire(file=dummy, exist=exist) + if (.not.exist) then + call rdpath(env%ptr%xtbpath, dummy, filename, exist) + if (.not.exist) then + filename = dummy + end if + end if + + allocate(xtb) + call newXTBCalculator(env%ptr, mol%ptr, xtb, filename) + + call env%ptr%check(exitRun) + if (exitRun) then + call env%ptr%error("Could not construct GFN2-xTB calculator", source) + return + end if + + call move_alloc(xtb, calc%ptr) + + end if + +end subroutine loadGFN2xTB_api + + +!> Add a solvation model to calculator (requires loaded parametrisation) +subroutine setSolvent_api(venv, vcalc, charptr, state, temperature, grid) & + & bind(C, name="xtb_setSolvent") + character(len=*), parameter :: source = 'xtb_api_setSolvent' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + character(kind=c_char), intent(in) :: charptr(*) + integer(c_int), intent(in), optional :: state + real(c_double), intent(in), optional :: temperature + integer(c_int), intent(in), optional :: grid + character(len=:), allocatable :: solvent + integer :: gsolvstate, nang + real(wp) :: temp + logical :: exitRun + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + if (.not.allocated(calc%ptr)) then + call env%ptr%error("No calculator loaded to add solvation model", source) + return + end if + + if (present(state)) then + gsolvstate = state + else + gsolvstate = 0 + end if + + if (present(temperature)) then + temp = temperature + else + temp = 298.15_wp + end if + + if (present(grid)) then + nang = grid + else + nang = 230 + end if + + call c_f_character(charptr, solvent) + + call addSolvationModel(env%ptr, calc%ptr, solvent, gsolvstate, temp, nang) + + call env%ptr%check(exitRun) + if (exitRun) then + call env%ptr%error("Could not add solvation model for '"//solvent//"'", & + & source) + return + end if + + end if + +end subroutine setSolvent_api + + +!> Unset the solvation model +subroutine releaseSolvent_api(venv, vcalc) & + & bind(C, name="xtb_releaseSolvent") + character(len=*), parameter :: source = 'xtb_api_setRelease' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + lgbsa = .false. + + end if + +end subroutine releaseSolvent_api + + +!> Add a external charge potential to calculator (only supported in GFN1/2-xTB) +subroutine setExternalCharges_api(venv, vcalc, npc, numbers, charges, positions) & + & bind(C, name="xtb_setExternalCharges") + character(len=*), parameter :: source = 'xtb_api_setExternalCharges' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + integer(c_int), intent(in) :: npc + integer(c_int), intent(in) :: numbers(*) + real(c_double), intent(in) :: charges(*) + real(c_double), intent(in) :: positions(3, *) + integer :: ii + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (npc <= 0) then + call env%ptr%error("Negative number of point charges provided", source) + return + end if + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + if (.not.allocated(calc%ptr)) then + call env%ptr%error("No calculator loaded to add external potential", & + & source) + return + end if + + select type(xtb => calc%ptr) + class default + call env%ptr%error("Calculator does not support external potentials", & + & source) + return + + type is(TxTBCalculator) + if (xtb%xtbData%level == 0) then + call env%ptr%error("GFN0-xTB does not support external potentials", & + & source) + return + end if + + call xtb%pcem%allocate(npc) + do ii = 1, npc + xtb%pcem%xyz(:, ii) = positions(:, ii) + xtb%pcem%gam(ii) = xtb%xtbData%coulomb%chemicalHardness(numbers(ii)) + xtb%pcem%q(ii) = charges(ii) + end do + + end select + + end if + +end subroutine setExternalCharges_api + + +!> Unset the external charge potential +subroutine releaseExternalCharges_api(venv, vcalc) & + & bind(C, name="xtb_releaseExternalCharges") + character(len=*), parameter :: source = 'xtb_api_releaseExternalCharges' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + ! There are only limited cases when we need to perform any action. + ! Deconstructors should behave nicely, therefore no errors on wrong input. + if (allocated(calc%ptr)) then + select type(xtb => calc%ptr) + type is(TxTBCalculator) + call xtb%pcem%deallocate + end select + end if + + end if + +end subroutine releaseExternalCharges_api + + +end module xtb_api_calculator diff --git a/src/api/environment.f90 b/src/api/environment.f90 new file mode 100644 index 000000000..ce97c8361 --- /dev/null +++ b/src/api/environment.f90 @@ -0,0 +1,205 @@ +! This file is part of xtb. +! +! Copyright (C) 2019-2020 Sebastian Ehlert +! +! xtb 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. +! +! xtb 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 xtb. If not, see . + +!> API for dealing with the calculation environment +module xtb_api_environment + use, intrinsic :: iso_c_binding + use xtb_mctc_accuracy, only : wp + use xtb_mctc_io + use xtb_api_utils + use xtb_type_environment + use xtb_xtb_calculator + implicit none + private + + public :: VEnvironment + public :: newEnvironment_api, delEnvironment_api, checkEnvironment_api, & + & showEnvironment_api, setOutput_api, releaseOutput_api, setVerbosity_api + + + !> Void pointer to calculation environment + type :: VEnvironment + type(TEnvironment) :: ptr + integer :: verbosity + end type VEnvironment + + +contains + + +!> Create new xtb calculation environment object +function newEnvironment_api() result(venv) & + & bind(C, name="xtb_newEnvironment") + type(VEnvironment), pointer :: env + type(c_ptr) :: venv + + call checkGlobalEnv + + allocate(env) + call init(env%ptr) + env%verbosity = 1 + venv = c_loc(env) + +end function newEnvironment_api + + +!> Delete a xtb calculation environment object +subroutine delEnvironment_api(venv) & + & bind(C, name="xtb_delEnvironment") + type(c_ptr), intent(inout) :: venv + type(VEnvironment), pointer :: env + + call releaseOutput_api(venv) + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + deallocate(env) + venv = c_null_ptr + end if + +end subroutine delEnvironment_api + + +!> Check current status of calculation environment +function checkEnvironment_api(venv) result(status) & + & bind(C, name="xtb_checkEnvironment") + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + integer(c_int) :: status + logical :: exitRun + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + call env%ptr%check(exitRun) + if (exitRun) then + status = 1 + else + status = 0 + end if + end if + +end function checkEnvironment_api + + +!> Show and empty error stack +subroutine showEnvironment_api(venv, charptr) & + & bind(C, name="xtb_showEnvironment") + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + character(kind=c_char), intent(in), optional :: charptr(*) + character(len=:, kind=c_char), allocatable :: message + + if (present(charptr)) then + call c_f_character(charptr, message) + else + message = '' + end if + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + call env%ptr%show(message) + end if + +end subroutine showEnvironment_api + + +!> Bind output from this environment +subroutine setOutput_api(venv, charptr) & + & bind(C, name="xtb_setOutput") + character(len=*), parameter :: source = 'xtb_api_setOutput' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + character(kind=c_char), intent(in) :: charptr(*) + character(len=:, kind=c_char), allocatable :: filename + integer :: stat + character(len=512) :: message + + call c_f_character(charptr, filename) + + call releaseOutput_api(venv) + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (filename == c_char_'-' .or. filename == c_char_'STDOUT') then + env%ptr%unit = stdout + else if (filename == c_char_'STDERR') then + env%ptr%unit = stderr + else + open(file=filename, newunit=env%ptr%unit, iostat=stat, iomsg=message) + if (stat /= 0) then + call env%ptr%error(trim(message), source) + end if + end if + end if + +end subroutine setOutput_api + + +!> Release output unit from this environment +subroutine releaseOutput_api(venv) & + & bind(C, name="xtb_releaseOutput") + character(len=*), parameter :: source = 'xtb_api_releaseOutput' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + integer :: stat + character(len=512) :: message + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.any(env%ptr%unit /= [-1, stdout, stderr])) then + close(unit=env%ptr%unit, iostat=stat, iomsg=message) + if (stat /= 0) then + call env%ptr%error(trim(message), source) + end if + end if + end if + +end subroutine releaseOutput_api + + +!> Set verbosity of calculation output +subroutine setVerbosity_api(venv, verbosity) & + & bind(C, name="xtb_setVerbosity") + character(len=*), parameter :: source = 'xtb_api_setVerbosity' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + integer, value :: verbosity + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.any(verbosity /= [0_c_int, 1_c_int, 2_c_int])) then + call env%ptr%warning("Verbosity specified is out of range", source) + end if + + env%verbosity = verbosity + end if + +end subroutine setVerbosity_api + + +end module xtb_api_environment diff --git a/src/api/interface.f90 b/src/api/interface.f90 index 3d6493e10..f0593ef3b 100644 --- a/src/api/interface.f90 +++ b/src/api/interface.f90 @@ -1,6 +1,6 @@ ! This file is part of xtb. ! -! Copyright (C) 2017-2020 Stefan Grimme +! Copyright (C) 2019-2020 Sebastian Ehlert ! ! xtb is free software: you can redistribute it and/or modify it under ! the terms of the GNU Lesser General Public License as published by @@ -15,995 +15,135 @@ ! You should have received a copy of the GNU Lesser General Public License ! along with xtb. If not, see . -!> Enduser bindings for the xtb API. This module tries to provide the most -! convenient interface to commonly used xtb methods. +!> Actual calculation interface module xtb_api_interface - use iso_c_binding - use xtb_api_structs - use xtb_api_utils - use xtb_api_preload + use, intrinsic :: iso_c_binding use xtb_mctc_accuracy, only : wp - use xtb_mctc_io, only : stdout - use xtb_type_environment, only : TEnvironment, init - use xtb_xtb_data - use xtb_xtb_gfn0 - use xtb_xtb_gfn1 - use xtb_xtb_gfn2 - + use xtb_api_calculator + use xtb_api_environment + use xtb_api_molecule + use xtb_api_results + use xtb_api_utils + use xtb_gfnff_calculator, only : TGFFCalculator + use xtb_scc_core, only : iniqshell + use xtb_type_data, only : scc_results + use xtb_xtb_calculator, only : TxTBCalculator + use xtb_main_setup, only : newWavefunction implicit none + private -contains - -!> Unsafe xTB calculation, this calculation mode requires to setup at least -! the parametrisation and the molecular structure data. -! -! In case null-pointers are passed instead of the basisset and wavefunction, -! those will be constructed on-the-fly and deleted afterwards. -! Properties like energy, gradient and stress tensor are returned directly, -! while other properties can be obtained from querying the wavefunction. -integer(c_int) function xtb_calculation_api & - & (c_mol, c_param, c_basis, c_wfn, c_pcem, c_opt, c_output, & - & c_energy, c_gradient, c_stress) & - & result(status) bind(C, name="xTB_calculation") - use xtb_setparam, only : gfn_method - use xtb_type_environment, only : TEnvironment, init - use xtb_type_molecule - use xtb_type_basisset - use xtb_type_wavefunction - use xtb_type_options - use xtb_type_data - use xtb_type_pcem - use xtb_scf - use xtb_peeq - !> Opaque pointer to the molecular structure data - type(c_ptr), value, intent(in) :: c_mol - !> Actual molecular structure data - type(TMolecule), pointer :: mol - !> Opaque pointer to the paramatrisation, currently not used - type(c_ptr), value, intent(in) :: c_param - !> Opaque pointer to the tight binding basisset - type(c_ptr), value, intent(in) :: c_basis - !> Actual basisset, constructed on-the-fly if non is given above - type(TBasisset), pointer :: basis - !> Opaque pointer to the tight binding wavefunction - type(c_ptr), value, intent(in) :: c_wfn - !> Actual wavefunction, constructed on-the-fly if non is given above - type(TWavefunction), pointer :: wfn - !> Opaque pointer to the external potential, currently not used - type(c_ptr), value, intent(in) :: c_pcem - !> Options for this calculation run as transparent structure - type(c_scc_options), intent(in) :: c_opt - !> Converted representation of the calculation options - type(scc_options) :: opt - !> File name for output, can also be "-" for standard output - character(kind=c_char), intent(in) :: c_output(*) - !> Converted representation of the file name - character(len=:), allocatable :: output - !> Pointer to store the total energy (in Hartree) - real(c_double), intent(out), optional :: c_energy - !> Pointer to store the molecular gradient (in Hartree/Bohr) - real(c_double), intent(out), optional :: c_gradient(3, *) - !> Pointer to store the stress tensor (in Hartree/Bohr³) - real(c_double), intent(out), optional :: c_stress(3, 3) - - type(TEnvironment) :: env - real(wp) :: energy, sigma(3, 3), egap - real(wp), allocatable :: gradient(:, :) - type(scc_results) :: res - type(tb_pcem) :: pcem - integer(c_int) :: stat_basis - logical :: exist, sane, exitRun + public :: singlepoint_api - call init(env) - status = 1 - ! perform some sanity checks first, if they fail, quickly return - if (.not.(check_xTB_init() .and. c_associated(c_mol))) return - call c_f_pointer(c_mol, mol) - if (mol%n <= 0) return - - energy = 0.0_wp - sigma = 0.0_wp - allocate(gradient(3, mol%n), source=0.0_wp) - - ! start converting transparent C-types - opt = c_opt - - ! open a unit for IO - call c_open_file(env%unit, c_output, opt%prlevel) +contains - ! print the xtb banner with version number and compilation date - if (opt%prlevel > 2) then - call xtb_header(env%unit) - call disclamer(env%unit) - call citation(env%unit) - endif - ! handle basisset - if (c_associated(c_basis)) then - call c_f_pointer(c_basis, basis) - else - allocate(basis) - call new_xtb_basisset(mol, basis, stat_basis) - if (stat_basis /= 0) then - ! don't call finalize here - deallocate(basis) - status = stat_basis +subroutine singlepoint_api(venv, vmol, vcalc, vres) & + & bind(C, name="xtb_singlepoint") + character(len=*), parameter :: source = 'xtb_api_singlepoint' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vmol + type(VMolecule), pointer :: mol + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + type(c_ptr), value :: vres + type(VResults), pointer :: res + type(scc_results) :: spRes + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vmol)) then + call env%ptr%error("Molecular structure data is not allocated", source) return end if - end if + call c_f_pointer(vmol, mol) - ! handle wavefunction - if (c_associated(c_wfn)) then - call c_f_pointer(c_wfn, wfn) - if (.not.verify_xtb_wavefunction(basis, wfn)) then - call finalize + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) return end if - else - allocate(wfn) - call wfn%allocate(basis%n, basis%nshell, basis%nao) - wfn%nel = nint(sum(mol%z) - mol%chrg) - wfn%nopen = mol%uhf - if (mod(wfn%nopen, 2) == 0 .and. mod(wfn%nel, 2) /= 0) wfn%nopen = 1 - if (mod(wfn%nopen, 2) /= 0 .and. mod(wfn%nel, 2) == 0) wfn%nopen = 0 - end if - - ! perform actual calculation - select case(gfn_method) - case(0) - call peeq & - & (env, mol, wfn, basis, global_data, & - & egap, opt%etemp, opt%prlevel, .false., opt%ccm, opt%acc, & - & energy, gradient, sigma, res) - case(1) - call scf & - & (env, mol, wfn, basis, pcem, global_data, & - & egap, opt%etemp, opt%maxiter, opt%prlevel, .false., .false., opt%acc, & - & energy, gradient, res) - case(2) - call scf & - & (env, mol, wfn, basis, pcem, global_data, & - & egap, opt%etemp, opt%maxiter, opt%prlevel, .false., .false., opt%acc, & - & energy, gradient, res) - case default - status = 3 - call finalize - return - end select - - call env%check(exitRun) - if (exitRun) then - call env%show("Single point calculator terminated") - call finalize - status = 4 - return - end if - - ! check if the MCTC environment is still sane, if not tell the caller - call mctc_sanity(sane) - if (.not.sane) then - status = 4 - call finalize - return - endif - - ! return values to C - if (present(c_energy)) c_energy = energy - if (present(c_gradient)) c_gradient(1:3, 1:mol%n) = gradient - if (present(c_stress)) c_stress = sigma / mol%volume - - ! cleanup and return - call finalize - status = 0 - -contains -subroutine finalize - if (env%unit /= stdout) close(env%unit) - if (.not.c_associated(c_basis)) then - call basis%deallocate - deallocate(basis) - nullify(basis) - end if - if (.not.c_associated(c_wfn)) then - call wfn%deallocate - deallocate(wfn) - nullify(wfn) - end if -end subroutine finalize -end function xtb_calculation_api - -function peeq_api & - & (natoms,attyp,charge,uhf,coord,lattice,pbc,opt_in,file_in, & - & etot,grad,stress,glat) & - & result(status) bind(C,name="GFN0_PBC_calculation") - - use xtb_type_molecule - use xtb_type_param - use xtb_type_options - - use xtb_calculators - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - real(c_double), intent(in) :: lattice(3,3) - logical(c_bool),intent(in) :: pbc(3) - type(c_peeq_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int) :: status - - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - real(c_double),intent(out) :: stress(3,3) - real(c_double),intent(out) :: glat(3,3) - - type(TMolecule) :: mol - type(peeq_options) :: opt - type(TEnvironment) :: env - - character(len=:),allocatable :: outfile - - logical :: sane, exitRun - integer :: i - real(wp) :: energy - real(wp) :: hl_gap - real(wp) :: lattice_gradient(3,3) - real(wp) :: stress_tensor(3,3) - real(wp),allocatable :: gradient(:,:) - - status = load_xtb_parameters_api(0_c_int) - if (status /= 0) return - - call init(env) - - ! ==================================================================== - ! STEP 2: convert the options from C struct to actual Fortran type - ! ==================================================================== - opt = opt_in - - ! open a unit for IO - call c_open_file(env%unit, file_in, opt%prlevel) - - if (opt%prlevel > 2) then - ! print the xtb banner with version number and compilation date - call xtb_header(env%unit) - ! make sure you cannot blame us for destroying your computer - call disclamer(env%unit) - ! how to cite this program - call citation(env%unit) - endif - - ! ==================================================================== - ! STEP 3: aquire the molecular structure and fill with data from C - ! ==================================================================== - mol = TMolecule(natoms, attyp, coord, charge, uhf, lattice, pbc) - - ! ==================================================================== - ! STEP 4: reserve some memory - ! ==================================================================== - allocate(gradient(3,mol%n)) - energy = 0.0_wp - gradient = 0.0_wp - - ! ==================================================================== - ! STEP 5: call the actual Fortran API to perform the calculation - ! ==================================================================== - ! shut down fatal errors from the MCTC library, so it will not kill the caller - call mctc_mute - - call gfn0_calculation & - (env%unit,env,opt,mol,hl_gap,energy,gradient,stress_tensor,lattice_gradient) - - call env%check(exitRun) - if (exitRun) then - call env%show("Single point calculator terminated") - call finalize - status = 1 - return - end if - - ! check if the MCTC environment is still sane, if not tell the caller - call mctc_sanity(sane) - if (.not.sane) then - call finalize - status = 1 - return - endif - - ! ==================================================================== - ! STEP 6: finally return values to C - ! ==================================================================== - call c_return(etot, energy) - call c_return(grad, gradient) - call c_return(glat, lattice_gradient) - call c_return(stress, stress_tensor) - - call finalize - - status = 0 - -contains - subroutine finalize - call mol%deallocate - deallocate(gradient) - if (env%unit.ne.stdout) close(env%unit) - end subroutine finalize -end function peeq_api - -function gfn2_api & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in, & - & etot,grad,dipole,q,dipm,qp,wbo) & - & result(status) bind(C,name="GFN2_calculation") - - use xtb_type_options - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - type(c_scc_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int) :: status - - real(c_double),intent(out) :: q(natoms) - real(c_double),intent(out) :: wbo(natoms,natoms) - real(c_double),intent(out) :: dipm(3,natoms) - real(c_double),intent(out) :: qp(6,natoms) - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - real(c_double),intent(out) :: dipole(3) - - status = load_xtb_parameters_api(2_c_int) - if (status /= 0) return - - call mctc_init('peeq',10,.true.) - - status = gfn12_calc_impl & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in,etot,grad,dipole,q,wbo,dipm,qp) - -end function gfn2_api - -function gfn1_api & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in,etot,grad,dipole,q,wbo) & - & result(status) bind(C,name="GFN1_calculation") - - use xtb_type_options - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - type(c_scc_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int) :: status - - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - real(c_double),intent(out) :: q(natoms) - real(c_double),intent(out) :: wbo(natoms,natoms) - real(c_double),intent(out) :: dipole(3) - - status = load_xtb_parameters_api(1_c_int) - if (status /= 0) return - - call mctc_init('peeq',10,.true.) - - status = gfn12_calc_impl & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in,etot,grad,dipole,q,wbo) - -end function gfn1_api - -function gfn12_calc_impl & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in,etot,grad,dipole,q,wbo,dipm,qp) & - & result(status) - - use xtb_type_molecule - use xtb_type_param - use xtb_type_options - use xtb_type_pcem - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_data - - use xtb_setparam, only : ngrida, gfn_method - use xtb_scf - use xtb_solv_gbobc - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - type(c_scc_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int) :: status - - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - real(c_double),intent(out) :: q(natoms) - real(c_double),intent(out) :: wbo(natoms,natoms) - real(c_double),intent(out) :: dipole(3) - real(c_double),intent(out),optional :: dipm(3,natoms) - real(c_double),intent(out),optional :: qp(6,natoms) - - type(TMolecule) :: mol - type(scc_options) :: opt - type(TEnvironment) :: env - type(tb_pcem) :: pcem - type(TWavefunction) :: wfn - type(TBasisset) :: basis - type(scc_results) :: res - - character(len=:),allocatable :: outfile - - logical :: sane, exitRun - integer :: i - real(wp) :: energy - real(wp) :: hl_gap - real(wp),allocatable :: gradient(:,:) - integer(c_int) :: stat_basis - - call init(env) - - ! convert the options from C struct to actual Fortran type - opt = opt_in - - ! open a unit for IO - call c_open_file(env%unit, file_in, opt%prlevel) - - if (opt%prlevel > 2) then - ! print the xtb banner with version number and compilation date - call xtb_header(env%unit) - ! make sure you cannot blame us for destroying your computer - call disclamer(env%unit) - ! how to cite this program - call citation(env%unit) - endif - - ! aquire the molecular structure and fill with data from C - mol = TMolecule(natoms, attyp, coord, charge, uhf) - - ! reserve some memory - allocate(gradient(3,mol%n)) - energy = 0.0_wp - gradient = 0.0_wp - - ! setup solvent model - lgbsa = len_trim(opt%solvent).gt.0 .and. opt%solvent.ne."none" - if (lgbsa) then - call init_gbsa(env%unit,trim(opt%solvent),0,opt%etemp,gfn_method,ngrida, & - & opt%prlevel > 0) - endif - - ! call the actual Fortran API to perform the calculation - call new_xtb_basisset(mol, basis, stat_basis) - if (stat_basis /= 0) then - status = stat_basis - call finalize - return - end if - - call wfn%allocate(basis%n, basis%nshell, basis%nao) - wfn%nel = nint(sum(mol%z) - mol%chrg) - wfn%nopen = mol%uhf - if (mod(wfn%nopen, 2) == 0 .and. mod(wfn%nel, 2) /= 0) wfn%nopen = 1 - if (mod(wfn%nopen, 2) /= 0 .and. mod(wfn%nel, 2) == 0) wfn%nopen = 0 - - call eeq_guess_wavefunction(env, mol, wfn, global_data) - - call scf & - & (env, mol, wfn, basis, pcem, global_data, & - & hl_gap, opt%etemp, opt%maxiter, opt%prlevel, .false., .false., opt%acc, & - & energy, gradient, res) - - call env%check(exitRun) - if (exitRun) then - call env%show("Single point calculator terminated") - call finalize - status = 1 - return - end if - - ! check if the MCTC environment is still sane, if not tell the caller - call mctc_sanity(sane) - if (.not.sane) then - call finalize - status = 1 - return - endif - - ! finally return values to C - call c_return(etot, energy) - call c_return(grad, gradient) - call c_return(q, wfn%q) - call c_return(wbo, wfn%wbo) - call c_return(dipole, sum(wfn%dipm,dim=2) + matmul(mol%xyz,wfn%q)) - if (present(dipm)) call c_return(dipm, wfn%dipm) - if (present(qp)) call c_return(qp, wfn%qp) - - call finalize - - status = 0 - -contains - subroutine finalize - call mol%deallocate - call wfn%deallocate - call basis%deallocate - deallocate(gradient) - if (env%unit.ne.stdout) close(env%unit) - end subroutine finalize -end function gfn12_calc_impl - -function gfn0_api & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in,etot,grad) & - & result(status) bind(C,name="GFN0_calculation") - - use xtb_type_molecule - use xtb_type_param - use xtb_type_options + call c_f_pointer(vcalc, calc) - use xtb_calculators - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - type(c_peeq_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int) :: status - - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - - type(TMolecule) :: mol - type(peeq_options) :: opt - type(TEnvironment) :: env - - character(len=:),allocatable :: outfile - - logical :: sane, exitRun - integer :: i - real(wp) :: energy - real(wp) :: hl_gap - real(wp) :: dum(3,3) - real(wp),allocatable :: gradient(:,:) - - status = load_xtb_parameters_api(0_c_int) - if (status /= 0) return - - call mctc_init('peeq',10,.true.) - - call init(env) - - ! ==================================================================== - ! STEP 2: convert the options from C struct to actual Fortran type - ! ==================================================================== - opt = opt_in - - ! open a unit for IO - call c_open_file(env%unit, file_in, opt%prlevel) - - if (opt%prlevel > 2) then - ! print the xtb banner with version number and compilation date - call xtb_header(env%unit) - ! make sure you cannot blame us for destroying your computer - call disclamer(env%unit) - ! how to cite this program - call citation(env%unit) - endif - - ! ==================================================================== - ! STEP 3: aquire the molecular structure and fill with data from C - ! ==================================================================== - mol = TMolecule(natoms, attyp, coord, charge, uhf) - - ! ==================================================================== - ! STEP 4: reserve some memory - ! ==================================================================== - allocate(gradient(3,mol%n)) - energy = 0.0_wp - gradient = 0.0_wp - - ! ==================================================================== - ! STEP 5: call the actual Fortran API to perform the calculation - ! ==================================================================== - ! shut down fatal errors from the MCTC library, so it will not kill the caller - call mctc_mute - - call gfn0_calculation & - (env%unit,env,opt,mol,hl_gap,energy,gradient,dum,dum) - - call env%check(exitRun) - if (exitRun) then - call env%show("Single point calculator terminated") - call finalize - status = 1 - return - end if - - ! check if the MCTC environment is still sane, if not tell the caller - call mctc_sanity(sane) - if (.not.sane) then - call finalize - status = 1 - return - endif - - ! ==================================================================== - ! STEP 6: finally return values to C - ! ==================================================================== - call c_return(etot, energy) - call c_return(grad, gradient) - - call finalize - - status = 0 - -contains - subroutine finalize - call mol%deallocate - deallocate(gradient) - if (env%unit.ne.stdout) close(env%unit) - end subroutine finalize -end function gfn0_api - -function gfn2_pcem_api & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in, & - & npc,pc_q,pc_at,pc_gam,pc_coord,etot,grad,pc_grad) & - & result(status) bind(C,name="GFN2_QMMM_calculation") - - use xtb_type_molecule - use xtb_type_param - use xtb_type_options - use xtb_type_pcem - use xtb_type_wavefunction - - use xtb_calculators - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - type(c_scc_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int), intent(in) :: npc - real(c_double), intent(in) :: pc_q(npc) - integer(c_int), intent(in) :: pc_at(npc) - real(c_double), intent(in) :: pc_gam(npc) - real(c_double), intent(in) :: pc_coord(3,npc) - - integer(c_int) :: status - - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - real(c_double),intent(out) :: pc_grad(3,npc) - - status = load_xtb_parameters_api(2_c_int) - if (status /= 0) return - - call mctc_init('peeq',10,.true.) - - status = gfn12_pcem_impl & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in, & - & npc,pc_q,pc_at,pc_gam,pc_coord,etot,grad,pc_grad) - -end function gfn2_pcem_api - -function gfn1_pcem_api & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in, & - & npc,pc_q,pc_at,pc_gam,pc_coord,etot,grad,pc_grad) & - & result(status) bind(C,name="GFN1_QMMM_calculation") - - use xtb_type_options - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - type(c_scc_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int), intent(in) :: npc - real(c_double), intent(in) :: pc_q(npc) - integer(c_int), intent(in) :: pc_at(npc) - real(c_double), intent(in) :: pc_gam(npc) - real(c_double), intent(in) :: pc_coord(3,npc) - - integer(c_int) :: status - - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - real(c_double),intent(out) :: pc_grad(3,npc) - - status = load_xtb_parameters_api(1_c_int) - if (status /= 0) return - - call mctc_init('peeq',10,.true.) - - status = gfn12_pcem_impl & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in, & - & npc,pc_q,pc_at,pc_gam,pc_coord,etot,grad,pc_grad) - -end function gfn1_pcem_api - -function gfn12_pcem_impl & - & (natoms,attyp,charge,uhf,coord,opt_in,file_in, & - & npc,pc_q,pc_at,pc_gam,pc_coord,etot,grad,pc_grad) & - & result(status) - - use xtb_type_molecule - use xtb_type_param - use xtb_type_options - use xtb_type_pcem - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_data - - use xtb_setparam, only : gfn_method, ngrida - use xtb_scf - use xtb_solv_gbobc - - use xtb_calculators - - implicit none - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: coord(3,natoms) - type(c_scc_options), intent(in) :: opt_in - character(kind=c_char),intent(in) :: file_in(*) - - integer(c_int), intent(in) :: npc - real(c_double), intent(in) :: pc_q(npc) - integer(c_int), intent(in) :: pc_at(npc) - real(c_double), intent(in) :: pc_gam(npc) - real(c_double), intent(in) :: pc_coord(3,npc) - - integer(c_int) :: status - - real(c_double),intent(out) :: etot - real(c_double),intent(out) :: grad(3,natoms) - real(c_double),intent(out) :: pc_grad(3,npc) - - type(TMolecule) :: mol - type(scc_options) :: opt - type(TWavefunction) :: wfn - type(TBasisset) :: basis - type(TEnvironment) :: env - type(tb_pcem) :: pcem - - character(len=:),allocatable :: outfile - - integer(c_int) :: stat_basis - type(scc_results) :: res - logical :: sane, exitRun - integer :: i - real(wp) :: energy - real(wp) :: hl_gap - real(wp),allocatable :: gradient(:,:) - - call init(env) - - ! convert the options from C struct to actual Fortran type - opt = opt_in - - ! open a unit for IO - call c_open_file(env%unit, file_in, opt%prlevel) + if (.not.allocated(calc%ptr)) then + call env%ptr%error("No calculator loaded for single point", & + & source) + return + end if - if (opt%prlevel > 2) then - ! print the xtb banner with version number and compilation date - call xtb_header(env%unit) - ! make sure you cannot blame us for destroying your computer - call disclamer(env%unit) - ! how to cite this program - call citation(env%unit) - endif + if (.not.c_associated(vres)) then + call env%ptr%error("Calculation results are not allocated", source) + return + end if + call c_f_pointer(vres, res) + + ! check cache, automatically invalidate missmatched data + if (allocated(res%wfn)) then + select type(xtb => calc%ptr) + type is(TxTBCalculator) + if (res%wfn%n /= mol%ptr%n .or. res%wfn%n /= xtb%basis%n .or. & + & res%wfn%nao /= xtb%basis%nao .or. & + & res%wfn%nshell /= xtb%basis%nshell) then + deallocate(res%wfn) + end if + end select + end if - ! aquire the molecular structure and fill with data from C - mol = TMolecule(natoms, attyp, coord, charge, uhf) + if (.not.allocated(res%wfn)) then + allocate(res%wfn) + ! in case of a new wavefunction cache we have to perform an initial guess + select type(xtb => calc%ptr) + type is(TxTBCalculator) + call newWavefunction(env%ptr, mol%ptr, xtb, res%wfn) + end select + end if - call pcem%allocate(npc) - pcem%q = pc_q - pcem%xyz = pc_coord - where(pc_at > 0) - pcem%gam = global_data%coulomb%chemicalHardness(pc_at) - elsewhere - pcem%gam = pc_gam - endwhere + if (.not.allocated(res%energy)) then + allocate(res%energy) + end if - ! reserve some memory - allocate(gradient(3,mol%n)) - energy = 0.0_wp - gradient = 0.0_wp + if (.not.allocated(res%egap)) then + allocate(res%egap) + end if - ! setup solvent model - lgbsa = len_trim(opt%solvent).gt.0 .and. opt%solvent.ne."none" - if (lgbsa) then - call init_gbsa(env%unit,trim(opt%solvent),0,opt%etemp,gfn_method,ngrida,& - & opt%prlevel > 0) - endif + if (allocated(res%gradient)) then + if (any(shape(res%gradient) /= [3, mol%ptr%n])) then + call env%ptr%warning("Shape missmatch in gradient, reallocating", source) + deallocate(res%gradient) + end if + end if + if (.not.allocated(res%gradient)) then + allocate(res%gradient(3, mol%ptr%n)) + end if - ! call the actual Fortran API to perform the calculation - call new_xtb_basisset(mol, basis, stat_basis) - if (stat_basis /= 0) then - status = stat_basis - call finalize - return - end if + if (allocated(res%sigma)) then + if (any(shape(res%sigma) /= [3, 3])) then + call env%ptr%warning("Shape missmatch in virial, reallocating", source) + deallocate(res%sigma) + end if + end if + if (.not.allocated(res%sigma)) then + allocate(res%sigma(3, 3)) + end if - call wfn%allocate(basis%n, basis%nshell, basis%nao) - wfn%nel = nint(sum(mol%z) - mol%chrg) - wfn%nopen = mol%uhf - if (mod(wfn%nopen, 2) == 0 .and. mod(wfn%nel, 2) /= 0) wfn%nopen = 1 - if (mod(wfn%nopen, 2) /= 0 .and. mod(wfn%nel, 2) == 0) wfn%nopen = 0 + call calc%ptr%singlepoint(env%ptr, mol%ptr, res%wfn, env%verbosity, .true., & + & res%energy, res%gradient, res%sigma, res%egap, spRes) - call eeq_guess_wavefunction(env, mol, wfn, global_data) + ! invalidate cache for properties not produced in GFN-FF + select type(gfnff => calc%ptr) + type is(TGFFCalculator) + deallocate(res%wfn) + deallocate(res%egap) + deallocate(res%sigma) + end select - call scf & - & (env, mol, wfn, basis, pcem, global_data, & - & hl_gap, opt%etemp, opt%maxiter, opt%prlevel, .false., .false., opt%acc, & - & energy, gradient, res) + res%dipole = spRes%dipole - call env%check(exitRun) - if (exitRun) then - call env%show("Single point calculator terminated") - call finalize - status = 1 - return end if - ! check if the MCTC environment is still sane, if not tell the caller - call mctc_sanity(sane) - if (.not.sane) then - call finalize - status = 1 - return - endif +end subroutine singlepoint_api - ! finally return values to C - call c_return(etot, energy) - call c_return(grad, gradient) - call c_return(pc_grad, pcem%grd) - - call finalize - - status = 0 - -contains - subroutine finalize - call mol%deallocate - call pcem%deallocate - call wfn%deallocate - deallocate(gradient) - if (env%unit.ne.stdout) close(env%unit) - end subroutine finalize -end function gfn12_pcem_impl - -function gbsa_calculation_api & - & (natoms,attyp,coord, & - & solvent_in,reference,temperature,method,grid_size,file_in, & - & brad,sasa) & - & result(status) bind(C,name='GBSA_calculation') - - use xtb_mctc_constants - - use xtb_type_molecule - use xtb_type_options - - use xtb_solv_gbobc - - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: attyp(natoms) - real(c_double), intent(in) :: coord(3,natoms) - - character(kind=c_char), intent(in) :: file_in(*) - character(kind=c_char), intent(in) :: solvent_in(*) - integer(c_int), intent(in) :: reference - real(c_double), intent(in) :: temperature - integer(c_int), intent(in) :: method - integer(c_int), intent(in) :: grid_size - integer(c_int) :: status - - real(c_double), intent(out) :: brad(natoms) - real(c_double), intent(out) :: sasa(natoms) - - logical :: sane, exitRun - character(len=:), allocatable :: outfile - character(len=:), allocatable :: solvent - - type(TMolecule) :: mol - type(TSolvent) :: gbsa - type(TEnvironment) :: env - - call mctc_init('gbobc',10,.true.) - call init(env) - - call c_string_convert(solvent,solvent_in) - call c_string_convert(outfile,file_in) - - if (outfile.ne.'-') then - call open_file(env%unit,outfile,'w') - if (env%unit.eq.-1) then - env%unit = stdout - endif - else - env%unit = stdout - endif - - ! shut down fatal errors from the MCTC library, so it will not kill the caller - call mctc_mute - - call init_gbsa(env%unit,solvent,reference,temperature,method,grid_size,.true.) - - call mctc_sanity(sane) - if (.not.sane) then - call finalize - status = 1 - return - endif - - call mol%allocate(natoms) - ! get atomtypes, coordinates and total charge - mol%at = attyp - mol%xyz = coord - call mol%update - - call new_gbsa(gbsa,mol%n,mol%at) - ! initialize the neighbor list - call update_nnlist_gbsa(gbsa,mol%xyz,.false.) - ! compute Born radii - call compute_brad_sasa(gbsa,mol%xyz) - - call mctc_sanity(sane) - if (.not.sane) then - call finalize - status = 1 - return - endif - - call c_return(brad, gbsa%brad) - call c_return(sasa, gbsa%sasa*fourpi/gbsa%gamsasa) - - call finalize - status = 0 - -contains - subroutine finalize - call mol%deallocate - call deallocate_gbsa(gbsa) - if (env%unit.ne.stdout) call close_file(env%unit) - end subroutine finalize -end function gbsa_calculation_api end module xtb_api_interface diff --git a/src/api/meson.build b/src/api/meson.build index 3f26de514..4cfde33da 100644 --- a/src/api/meson.build +++ b/src/api/meson.build @@ -16,8 +16,10 @@ # along with xtb. If not, see . srcs += files( + 'calculator.f90', + 'environment.f90', + 'molecule.f90', 'interface.f90', - 'preload.f90', - 'structs.f90', + 'results.f90', 'utils.f90', ) diff --git a/src/api/molecule.f90 b/src/api/molecule.f90 new file mode 100644 index 000000000..6a27dfbd8 --- /dev/null +++ b/src/api/molecule.f90 @@ -0,0 +1,159 @@ +! This file is part of xtb. +! +! Copyright (C) 2019-2020 Sebastian Ehlert +! +! xtb 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. +! +! xtb 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 xtb. If not, see . + +!> API for dealing with molecular structure data +module xtb_api_molecule + use, intrinsic :: iso_c_binding + use xtb_mctc_accuracy, only : wp + use xtb_api_environment + use xtb_api_utils + use xtb_type_molecule + implicit none + private + + public :: VMolecule + public :: newMolecule_api, delMolecule_api, updateMolecule_api + + + !> Void pointer to molecular structure data + type :: VMolecule + type(TMolecule) :: ptr + end type VMolecule + + + !> Maximum accepted atomic number, actual calculator might accept less + integer, parameter :: maxElem = 118 + + +contains + + +function newMolecule_api(venv, natoms, numbers, positions, charge, uhf, lattice, & + & periodic) result(vmol) & + & bind(C, name="xtb_newMolecule") + character(len=*), parameter :: source = 'xtb_api_newMolecule' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + integer(c_int), intent(in) :: natoms + integer(c_int), intent(in) :: numbers(natoms) + real(c_double), intent(in) :: positions(3, natoms) + real(c_double), intent(in) :: charge + integer(c_int), intent(in) :: uhf + real(c_double), intent(in) :: lattice(3, 3) + logical(c_bool), intent(in) :: periodic(3) + type(VMolecule), pointer :: mol + type(c_ptr) :: vmol + integer(c_int) :: status + + vmol = c_null_ptr + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (natoms <= 0) then + call env%ptr%error("Number of atoms must be positive", source) + return + end if + + if (any(numbers <= 0) .or. any(numbers > maxElem)) then + call env%ptr%error("Invalid atomic number present", source) + return + end if + + allocate(mol) + + mol%ptr = TMolecule(natoms, numbers, positions, charge, uhf, & + & lattice, periodic) + + status = verifyMolecule(mol%ptr) + if (status /= 0_c_int) then + deallocate(mol) + call env%ptr%error("Could not generate molecular structure", source) + return + end if + + vmol = c_loc(mol) + end if + +end function newMolecule_api + + +subroutine delMolecule_api(vmol) & + & bind(C, name="xtb_delMolecule") + type(c_ptr), intent(inout) :: vmol + type(VMolecule), pointer :: mol + + if (c_associated(vmol)) then + call c_f_pointer(vmol, mol) + call checkGlobalEnv + + deallocate(mol) + vmol = c_null_ptr + end if + +end subroutine delMolecule_api + + +subroutine updateMolecule_api(venv, vmol, positions, lattice) & + & bind(C, name="xtb_updateMolecule") + character(len=*), parameter :: source = 'xtb_api_updateMolecule' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vmol + type(VMolecule), pointer :: mol + real(c_double), intent(in) :: positions(3, *) + real(c_double), intent(in), optional :: lattice(3, 3) + real(wp) :: latvecs(3, 3) + integer(c_int) :: status + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vmol)) then + call env%ptr%error("Molecular structure data is not allocated", source) + return + end if + + call c_f_pointer(vmol, mol) + + if (mol%ptr%n <= 0 .or. mol%ptr%nid <= 0 .or. .not.allocated(mol%ptr%at) & + & .or. .not.allocated(mol%ptr%id) .or. .not.allocated(mol%ptr%xyz)) then + call env%ptr%error("Invalid molecular structure data provided", source) + return + end if + + mol%ptr%xyz(:, :) = positions(1:3, 1:mol%ptr%n) + if (present(lattice)) then + mol%ptr%lattice(:, :) = lattice(1:3, 1:3) + end if + + call mol%ptr%update + + status = verifyMolecule(mol%ptr) + if (status /= 0_c_int) then + call env%ptr%error("Could not update molecular structure", source) + return + end if + + end if + +end subroutine updateMolecule_api + + +end module xtb_api_molecule diff --git a/src/api/preload.f90 b/src/api/preload.f90 deleted file mode 100644 index 350901adb..000000000 --- a/src/api/preload.f90 +++ /dev/null @@ -1,196 +0,0 @@ -! This file is part of xtb. -! -! Copyright (C) 2017-2020 Stefan Grimme -! -! xtb 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. -! -! xtb 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 xtb. If not, see . - -!> Preloading module of xtb to change the global state of the library. -module xtb_api_preload - use iso_c_binding - use xtb_mctc_accuracy, only : wp - use xtb_api_utils - use xtb_type_param, only : TxTBParameter - use xtb_type_environment, only : TEnvironment, init - use xtb_paramset - use xtb_xtb_data - implicit none - - type(TxTBData) :: global_data - -contains - -!> Load GFN-xTB parametrisation and create a lock, the file name is optional. -! The lock (gfn_method) is maintained by the shared libary. -integer(c_int) function load_xtb_parameters_api(gfn, filename) & - & result(status) bind(C, name='load_xTB_parameters') - use xtb_setparam, only : gfn_method - use xtb_readparam, only : readParam - !> GFN level to be loaded - integer(c_int), intent(in) :: gfn - !> File name for the parametrisation - character(kind=c_char), intent(in), optional :: filename(*) - - type(TEnvironment) :: env - character(len=:), allocatable :: fnv - type(TxTBParameter) :: globpar - integer :: ipar - logical :: exist - - call init(env) - - !$omp critical(xtb_load_api) - status = 1 - gfn_method = -1 - - call mctc_init('peeq',10,.true.) - call mctc_mute - - if (present(filename)) then - call c_string_convert(fnv, filename) - inquire(file=fnv, exist=exist) - if (exist) then - open(file=fnv, newunit=ipar) - call readParam(env, ipar, globpar, global_data, .true.) - close(ipar) - status = 0 - end if - else - select case(gfn) - case(0_c_int) - call load_xtb_parameters(env, 'param_gfn0-xtb.txt', globpar, status) - case(1_c_int) - call load_xtb_parameters(env, 'param_gfn1-xtb.txt', globpar, status) - case(2_c_int) - call load_xtb_parameters(env, 'param_gfn2-xtb.txt', globpar, status) - case default - status = 2 - end select - end if - - if (status == 0) then - gfn_method = gfn - end if - !$omp end critical(xtb_load_api) - -end function load_xtb_parameters_api - -subroutine load_xtb_parameters(env, p_fnv, globpar, status) - use xtb_mctc_systools - use xtb_paramset, only : use_parameterset - use xtb_readparam, only : readParam - type(TEnvironment), intent(inout) :: env - character(len=*), intent(in) :: p_fnv - integer, intent(out) :: status - character(len=:), allocatable :: xtbpath - character(len=:), allocatable :: fnv - type(TxTBParameter), intent(out) :: globpar - integer :: ipar - logical :: exist - status = 1 - - ! we will try an internal parameter file first to avoid IO - call use_parameterset(p_fnv, globpar, global_data, exist) - if (exist) then - status = 0 - else ! no luck, we have to fire up some IO to get our parameters - call rdvar('XTBPATH', xtbpath) - if (len(xtbpath) > 0) then - ! let's check if we can find the parameter file - call rdpath(xtbpath, p_fnv, fnv, exist) - end if - ! maybe the user provides a local parameter file, this was always - ! an option in `xtb', so we will give it a try - if (.not.exist) fnv = p_fnv - inquire(file=fnv, exist=exist) - if (exist) then - open(file=fnv, newunit=ipar) - call readParam(env, ipar, globpar, global_data, .true.) - close(ipar) - status = 0 - end if - end if - -end subroutine load_xtb_parameters - -!> allows loading custom solvent parameters in the shared library -function gbsa_model_preload_api & - & (epsv_in,smass_in,rhos_in,c1_in,rprobe_in,gshift_in,soset_in,dum_in, & - & gamscale_in,sx_in,tmp_in) & - & result(status) bind(C,name='GBSA_model_preload') - - use xtb_solv_gbobc, only : load_custom_parameters - - !> Dielectric data - real(c_double), intent(in) :: epsv_in - real(wp), allocatable :: epsv - !> Molar mass (g/mol) - real(c_double), intent(in) :: smass_in - real(wp), allocatable :: smass - !> Solvent density (g/cm^3) - real(c_double), intent(in) :: rhos_in - real(wp), allocatable :: rhos - !> Born radii - real(c_double), intent(in) :: c1_in - real(wp), allocatable :: c1 - !> Atomic surfaces - real(c_double), intent(in) :: rprobe_in - real(wp), allocatable :: rprobe - !> Gshift (gsolv=reference vs. gsolv) - real(c_double), intent(in) :: gshift_in - real(wp), allocatable :: gshift - !> offset parameter (fitted) - real(c_double), intent(in) :: soset_in - real(wp), allocatable :: soset - real(c_double), intent(in) :: dum_in - real(wp), allocatable :: dum - !> Surface tension (mN/m=dyn/cm) - real(c_double), intent(in) :: gamscale_in(94) - real(wp), allocatable :: gamscale(:) - !> dielectric descreening parameters - real(c_double), intent(in) :: sx_in(94) - real(wp), allocatable :: sx(:) - real(c_double), intent(in) :: tmp_in(94) - real(wp), allocatable :: tmp(:) - integer(c_int) :: status - - call c_get(epsv,epsv_in) - call c_get(smass,smass_in) - call c_get(rhos,rhos_in) - call c_get(c1,c1_in) - call c_get(rprobe,rprobe_in) - call c_get(gshift,gshift_in) - call c_get(soset,soset_in) - call c_get(dum,dum_in) - call c_get(gamscale,gamscale_in) - call c_get(sx,sx_in) - call c_get(tmp,tmp_in) - - call load_custom_parameters(epsv=epsv) - call load_custom_parameters(smass=smass) - call load_custom_parameters(rhos=rhos) - call load_custom_parameters(c1=c1) - call load_custom_parameters(rprobe=rprobe) - call load_custom_parameters(gshift=gshift) - call load_custom_parameters(soset=soset) - call load_custom_parameters(dum=dum) - call load_custom_parameters(gamscale=gamscale) - call load_custom_parameters(sx=sx) - call load_custom_parameters(tmp=tmp) - - status = 0 - -end function gbsa_model_preload_api - - -end module xtb_api_preload diff --git a/src/api/results.f90 b/src/api/results.f90 new file mode 100644 index 000000000..e6e06d3a4 --- /dev/null +++ b/src/api/results.f90 @@ -0,0 +1,270 @@ +! This file is part of xtb. +! +! Copyright (C) 2019-2020 Sebastian Ehlert +! +! xtb 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. +! +! xtb 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 xtb. If not, see . + +!> API for dealing with the calculation results +module xtb_api_results + use, intrinsic :: iso_c_binding + use xtb_mctc_accuracy, only : wp + use xtb_api_environment + use xtb_api_utils + use xtb_type_wavefunction + implicit none + private + + public :: VResults + public :: newResults_api, delResults_api, getEnergy_api, getGradient_api, & + & getVirial_api, getDipole_api, getCharges_api, getBondOrders_api + + + !> Void pointer to wavefunction and result objects + type :: VResults + type(TWavefunction), allocatable :: wfn + real(wp), allocatable :: energy + real(wp), allocatable :: gradient(:, :) + real(wp), allocatable :: sigma(:, :) + real(wp), allocatable :: dipole(:) + real(wp), allocatable :: egap + end type VResults + + +contains + + +!> Create new singlepoint results object +function newResults_api() result(vres) & + & bind(C, name="xtb_newResults") + type(VResults), pointer :: res + type(c_ptr) :: vres + + call checkGlobalEnv + + allocate(res) + vres = c_loc(res) + +end function newResults_api + + +!> Delete singlepoint results object +subroutine delResults_api(vres) & + & bind(C, name="xtb_delResults") + type(c_ptr), intent(inout) :: vres + type(VResults), pointer :: res + + call checkGlobalEnv + + if (c_associated(vres)) then + call c_f_pointer(vres, res) + deallocate(res) + vres = c_null_ptr + end if + +end subroutine delResults_api + + +!> Query singlepoint results object for energy +subroutine getEnergy_api(venv, vres, dptr) & + & bind(C, name="xtb_getEnergy") + character(len=*), parameter :: source = "xtb_api_getEnergy" + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vres + type(VResults), pointer :: res + real(c_double), intent(inout) :: dptr + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vres)) then + call env%ptr%error("Results object is not allocated", source) + return + end if + call c_f_pointer(vres, res) + + if (.not.allocated(res%energy)) then + call env%ptr%error("Energy is not available in results", source) + return + end if + + dptr = res%energy + + end if + +end subroutine getEnergy_api + + +!> Query singlepoint results object for gradient +subroutine getGradient_api(venv, vres, dptr) & + & bind(C, name="xtb_getGradient") + character(len=*), parameter :: source = "xtb_api_getGradient" + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vres + type(VResults), pointer :: res + real(c_double), intent(inout) :: dptr(3, *) + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vres)) then + call env%ptr%error("Results object is not allocated", source) + return + end if + call c_f_pointer(vres, res) + + if (.not.allocated(res%gradient)) then + call env%ptr%error("Gradient is not available in results", source) + return + end if + + dptr(1:3, 1:size(res%gradient, 2)) = res%gradient + + end if + +end subroutine getGradient_api + + +!> Query singlepoint results object for virial +subroutine getVirial_api(venv, vres, dptr) & + & bind(C, name="xtb_getVirial") + character(len=*), parameter :: source = "xtb_api_getVirial" + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vres + type(VResults), pointer :: res + real(c_double), intent(inout) :: dptr(3, 3) + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vres)) then + call env%ptr%error("Results object is not allocated", source) + return + end if + call c_f_pointer(vres, res) + + if (.not.allocated(res%sigma)) then + call env%ptr%error("Virial is not available in results", source) + return + end if + + dptr(1:3, 1:3) = res%sigma + + end if + +end subroutine getVirial_api + + +!> Query singlepoint results object for dipole moment +subroutine getDipole_api(venv, vres, dptr) & + & bind(C, name="xtb_getDipole") + character(len=*), parameter :: source = "xtb_api_getDipole" + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vres + type(VResults), pointer :: res + real(c_double), intent(inout) :: dptr(3) + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vres)) then + call env%ptr%error("Results object is not allocated", source) + return + end if + call c_f_pointer(vres, res) + + if (.not.allocated(res%dipole)) then + call env%ptr%error("Dipole moment is not available in results", source) + return + end if + + dptr(1:3) = res%dipole + + end if + +end subroutine getDipole_api + + +!> Query singlepoint results object for partial charges +subroutine getCharges_api(venv, vres, dptr) & + & bind(C, name="xtb_getCharges") + character(len=*), parameter :: source = "xtb_api_getCharges" + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vres + type(VResults), pointer :: res + real(c_double), intent(inout) :: dptr(*) + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vres)) then + call env%ptr%error("Results object is not allocated", source) + return + end if + call c_f_pointer(vres, res) + + if (.not.allocated(res%wfn)) then + call env%ptr%error("Partial charges are not available in results", source) + return + end if + + dptr(1:size(res%wfn%q)) = res%wfn%q + + end if + +end subroutine getCharges_api + + +!> Query singlepoint results object for bond orders +subroutine getBondOrders_api(venv, vres, dptr) & + & bind(C, name="xtb_getBondOrders") + character(len=*), parameter :: source = "xtb_api_getBondOrders" + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vres + type(VResults), pointer :: res + real(c_double), intent(inout) :: dptr(*) + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vres)) then + call env%ptr%error("Results object is not allocated", source) + return + end if + call c_f_pointer(vres, res) + + if (.not.allocated(res%wfn)) then + call env%ptr%error("Bond orders are not available in results", source) + return + end if + + dptr(1:size(res%wfn%wbo)) = reshape(res%wfn%wbo, [size(res%wfn%wbo)]) + + end if + +end subroutine getBondOrders_api + + +end module xtb_api_results diff --git a/src/api/structs.f90 b/src/api/structs.f90 deleted file mode 100644 index b05d2200c..000000000 --- a/src/api/structs.f90 +++ /dev/null @@ -1,328 +0,0 @@ -! This file is part of xtb. -! -! Copyright (C) 2017-2020 Stefan Grimme -! -! xtb 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. -! -! xtb 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 xtb. If not, see . - -!> Defines the interface to xtb datatypes and allows exporting them to C. -! Mostly those will be carried around as opaque void* pointer and can -! be manipulated with methods provided by the API. -module xtb_api_structs - use iso_c_binding - use xtb_mctc_accuracy, only : wp - use xtb_api_utils - implicit none - -contains - -!> Constructor for the molecular structure type, returns a nullptr in case -! the generation fails. All values have to be provided. -type(c_ptr) function new_xtb_molecule_api & - & (natoms, numbers, positions, charge, uhf, lattice, periodic) & - & result(c_mol) bind(C, name="new_xTB_molecule") - use xtb_type_molecule - integer(c_int), intent(in) :: natoms - integer(c_int), intent(in) :: numbers(natoms) - real(c_double), intent(in) :: positions(3, natoms) - real(c_double), intent(in) :: charge - integer(c_int), intent(in) :: uhf - real(c_double), intent(in) :: lattice(3, 3) - logical(c_bool), intent(in) :: periodic(3) - - type(TMolecule), pointer :: mol - c_mol = c_null_ptr - - allocate(mol) - - mol = TMolecule(natoms, numbers, positions, charge, uhf, lattice, periodic) - if (verify_xtb_molecule(mol) == 0) then - c_mol = c_loc(mol) - else - call mol%deallocate - deallocate(mol) - nullify(mol) - end if - -end function new_xtb_molecule_api - -!> Updates the molecular structure type. This routine checks for the allocation -! state of the opaque pointer before copying positions and lattice. -integer(c_int) function update_xtb_molecule_api & - & (c_mol, c_positions, c_lattice) & - & result(status) bind(C, name="update_xTB_molecule") - use xtb_type_molecule - type(c_ptr), value, intent(in) :: c_mol - type(c_ptr), value, intent(in) :: c_positions - type(c_ptr), value, intent(in) :: c_lattice - - type(TMolecule), pointer :: mol - real(c_double), pointer :: positions(:, :) - real(c_double), pointer :: lattice(:, :) - - call c_f_pointer(c_mol, mol) - if (mol%n > 0) then - if (c_associated(c_positions)) then - call c_f_pointer(c_positions, positions, [3, mol%n]) - mol%xyz = positions - end if - if (c_associated(c_lattice)) then - call c_f_pointer(c_lattice, lattice, [3, 3]) - mol%lattice = lattice - end if - call mol%update - ! check for very short distances - status = verify_xtb_molecule(mol) - else - status = -1 - end if - -end function update_xtb_molecule_api - -!> Deconstructor for the molecular structure type. -subroutine delete_xtb_molecule_api(c_mol) & - & bind(C, name="delete_xTB_molecule") - use xtb_type_molecule - type(c_ptr), value, intent(in) :: c_mol - type(TMolecule), pointer :: mol - - if (c_associated(c_mol)) then - call c_f_pointer(c_mol, mol) - call mol%deallocate - deallocate(mol) - nullify(mol) - end if - -end subroutine delete_xtb_molecule_api - -!> Constructor for the tight binding basisset. -type(c_ptr) function new_xtb_basisset_api(c_mol) & - & result(c_basis) bind(C, name="new_xTB_basisset") - use xtb_type_molecule - use xtb_type_basisset - type(c_ptr), value, intent(in) :: c_mol - type(TBasisset), pointer :: basis - type(TMolecule), pointer :: mol - integer(c_int) :: status - c_basis = c_null_ptr - - if (check_xtb_init()) then - - call c_f_pointer(c_mol, mol) - - if (mol%n > 0) then - allocate(basis) - - call new_xtb_basisset(mol, basis, status) - if (status == 0_c_int) then - c_basis = c_loc(basis) - else - deallocate(basis) - end if - end if - end if - -end function new_xtb_basisset_api - -!> Actual implementation of the constructor. -subroutine new_xtb_basisset(mol, basis, status) - use xtb_setparam, only : gfn_method - use xtb_type_molecule - use xtb_type_basisset - use xtb_basis - use xtb_xtb_data - use xtb_xtb_gfn0 - use xtb_xtb_gfn1 - use xtb_xtb_gfn2 - type(TMolecule), intent(in) :: mol - type(TBasisset), intent(inout) :: basis - integer(c_int), intent(out) :: status - type(TxTBData) :: xtbData - logical :: diff, okbas - status = 2 - select case(gfn_method) - case default - okbas = .false. - case(0) - call initGFN0(xtbData) - case(1) - call initGFN1(xtbData) - case(2) - call initGFN2(xtbData) - end select - call newBasisset(xtbData, mol%n, mol%at, basis, okbas) - if (okbas) then - status = 0 - else - call basis%deallocate - end if -end subroutine new_xtb_basisset - -!> Return dimensions of this tight binding basisset. -! For not allocated input zero is returned, every argument is optional. -subroutine dimensions_xtb_basisset_api & - & (c_basis, nat, nbf, nao, nsh) & - & bind(C, name="dimensions_xTB_basisset") - use xtb_type_basisset - type(c_ptr), value, intent(in) :: c_basis - type(TBasisset), pointer :: basis - integer(c_int), intent(out), optional :: nat - integer(c_int), intent(out), optional :: nbf - integer(c_int), intent(out), optional :: nao - integer(c_int), intent(out), optional :: nsh - - if (c_associated(c_basis)) then - call c_f_pointer(c_basis, basis) - call c_return(nat, basis%n) - call c_return(nbf, basis%nbf) - call c_return(nao, basis%nao) - call c_return(nsh, basis%nshell) - else - call c_return(nat, 0) - call c_return(nbf, 0) - call c_return(nao, 0) - call c_return(nsh, 0) - end if - -end subroutine dimensions_xtb_basisset_api - -!> Deconstructor for the tight binding wavefunction. -subroutine delete_xtb_basisset_api(c_basis) & - & bind(C, name="delete_xTB_basisset") - use xtb_type_basisset - type(c_ptr), value, intent(in) :: c_basis - type(TBasisset), pointer :: basis - - if (c_associated(c_basis)) then - call c_f_pointer(c_basis, basis) - call basis%deallocate - deallocate(basis) - nullify(basis) - end if - -end subroutine delete_xtb_basisset_api - -!> Constructor for the tight binding wavefunction. -type(c_ptr) function new_xtb_wavefunction_api & - & (c_mol, c_basis) & - & result(c_wfn) bind(C, name="new_xTB_wavefunction") - use xtb_type_molecule - use xtb_type_basisset - use xtb_type_wavefunction - type(c_ptr), value, intent(in) :: c_mol - type(c_ptr), value, intent(in) :: c_basis - type(TMolecule), pointer :: mol - type(TBasisset), pointer :: basis - type(TWavefunction), pointer :: wfn - c_wfn = c_null_ptr - - if (check_xtb_init()) then - if (c_associated(c_basis) .and. c_associated(c_mol)) then - call c_f_pointer(c_mol, mol) - call c_f_pointer(c_basis, basis) - if (mol%n > 0 .and. basis%n > 0 .and. basis%nbf > 0 .and. & - &basis%nao > 0 .and. basis%nshell > 0) then - allocate(wfn) - - call wfn%allocate(basis%n, basis%nshell, basis%nao) - - wfn%nel = nint(sum(mol%z) - mol%chrg) - wfn%nopen = mol%uhf - if (mod(wfn%nopen, 2) == 0 .and. mod(wfn%nel, 2) /= 0) wfn%nopen = 1 - if (mod(wfn%nopen, 2) /= 0 .and. mod(wfn%nel, 2) == 0) wfn%nopen = 0 - c_wfn = c_loc(wfn) - end if - end if - end if - -end function new_xtb_wavefunction_api - -!> Return dimensions of this tight binding wavefunction. -subroutine dimensions_xtb_wavefunction_api & - & (c_wfn, nat, nao, nsh) & - & bind(C, name="dimensions_xTB_wavefunction") - use xtb_type_wavefunction - type(c_ptr), value, intent(in) :: c_wfn - type(TWavefunction), pointer :: wfn - integer(c_int), intent(out), optional :: nat - integer(c_int), intent(out), optional :: nao - integer(c_int), intent(out), optional :: nsh - - if (c_associated(c_wfn)) then - call c_f_pointer(c_wfn, wfn) - call c_return(nat, wfn%n) - call c_return(nao, wfn%nao) - call c_return(nsh, wfn%nshell) - else - call c_return(nat, 0) - call c_return(nao, 0) - call c_return(nsh, 0) - end if - -end subroutine dimensions_xtb_wavefunction_api - -!> Query the wavefunction properties, pass NULL/nullptr in case you want -! a query to be ignored. -! In case the wavefunction is not allocated nothing is returned. -subroutine query_xtb_wavefunction_api & - & (c_wfn, charges, dipoles, quadrupoles, bond_orders, & - & hl_gap, orbital_energies) & - & bind(C, name="query_xTB_wavefunction") - use xtb_mctc_convert, only : evtoau - use xtb_type_wavefunction - type(c_ptr), value, intent(in) :: c_wfn - type(TWavefunction), pointer :: wfn - real(c_double), intent(out), optional :: charges(*) - real(c_double), intent(out), optional :: dipoles(3, *) - real(c_double), intent(out), optional :: quadrupoles(6, *) - real(c_double), intent(out), optional :: bond_orders(*) - real(c_double), intent(out), optional :: hl_gap - real(c_double), intent(out), optional :: orbital_energies(*) - - if (c_associated(c_wfn)) then - call c_f_pointer(c_wfn, wfn) - if (wfn%n > 0 .and. wfn%nao > 0 .and. wfn%nshell > 0) then - if (present(charges)) charges(1:wfn%n) = wfn%q - if (present(dipoles)) dipoles(1:3, 1:wfn%n) = wfn%dipm - if (present(quadrupoles)) quadrupoles(1:6, 1:wfn%n) = wfn%qp - if (present(bond_orders)) then - bond_orders(1:wfn%n**2) = reshape(wfn%wbo, [wfn%n**2]) - end if - if (present(orbital_energies)) then - orbital_energies(1:wfn%nao) = wfn%emo * evtoau - end if - if (present(hl_gap) .and. wfn%ihomo+1 <= wfn%nao .and. wfn%ihomo >= 1) then - hl_gap = (wfn%emo(wfn%ihomo+1) - wfn%emo(wfn%ihomo)) * evtoau - end if - end if - end if - -end subroutine query_xtb_wavefunction_api - -!> Deconstructor for the tight binding basisset. -subroutine delete_xtb_wavefunction_api(c_wfn) & - & bind(C, name="delete_xTB_wavefunction") - use xtb_type_wavefunction - type(c_ptr), value, intent(in) :: c_wfn - type(TWavefunction), pointer :: wfn - - if (c_associated(c_wfn)) then - call c_f_pointer(c_wfn, wfn) - call wfn%deallocate - deallocate(wfn) - nullify(wfn) - end if - -end subroutine delete_xtb_wavefunction_api - -end module xtb_api_structs diff --git a/src/api/utils.f90 b/src/api/utils.f90 index 32242ee62..168e3ee56 100644 --- a/src/api/utils.f90 +++ b/src/api/utils.f90 @@ -1,6 +1,6 @@ ! This file is part of xtb. ! -! Copyright (C) 2017-2020 Stefan Grimme +! Copyright (C) 2019-2020 Sebastian Ehlert ! ! xtb is free software: you can redistribute it and/or modify it under ! the terms of the GNU Lesser General Public License as published by @@ -15,95 +15,40 @@ ! You should have received a copy of the GNU Lesser General Public License ! along with xtb. If not, see . -!> Utilities to make working with C datatypes as convenient as with Fortran -! datatypes. +!> Utilities to work with data types from C module xtb_api_utils - use iso_c_binding + use, intrinsic :: iso_c_binding use xtb_mctc_accuracy, only : wp - use xtb_chargemodel - use xtb_xtb_data + use xtb_type_molecule, only : TMolecule + use xtb_type_environment, only : init implicit none + private - interface c_return - module procedure :: c_return_int_0d - module procedure :: c_return_double_0d - module procedure :: c_return_double_1d - module procedure :: c_return_double_2d - end interface c_return - - interface c_return_alloc - module procedure :: c_return_double_0d_alloc - module procedure :: c_return_double_1d_alloc - module procedure :: c_return_double_2d_alloc - end interface c_return_alloc - - interface c_get - module procedure :: c_get_double_0d_alloc - module procedure :: c_get_double_1d_alloc - module procedure :: c_get_double_2d_alloc - module procedure :: c_get_double_0d - module procedure :: c_get_double_1d - module procedure :: c_get_double_2d - module procedure :: c_get_int_0d - module procedure :: c_get_int_1d - module procedure :: c_get_int_2d - end interface c_get + public :: c_f_character, verifyMolecule, checkGlobalEnv contains -!> Check if xTB was correctly initialized. -logical(c_bool) function check_xtb_init() & - & result(status) bind(C, name="check_xTB_init") - use xtb_setparam, only : gfn_method - status = gfn_method >= 0 .and. gfn_method <= 2 -end function check_xtb_init -!> Check if xTB was locked on the correct parametrisation. -logical(c_bool) function check_xtb_lock(gfn) & - & result(status) bind(C, name="check_xTB_lock") - use xtb_setparam, only : gfn_method - !> Requested GFN level - integer(c_int), intent(in) :: gfn - status = gfn == gfn_method -end function check_xtb_lock +subroutine c_f_character(rhs, lhs) + character(kind=c_char), intent(in) :: rhs(*) + !> Resulting Fortran string + character(len=:, kind=c_char), allocatable, intent(out) :: lhs -!> Release the lock on the xTB parametrisation. -subroutine reset_xtb_lock() bind(C, name="reset_xTB_lock") - use xtb_setparam, only : gfn_method - !$omp critical(xtb_load_api) - gfn_method = -1 - !$omp end critical(xtb_load_api) -end subroutine + integer :: ii -subroutine eeq_guess_wavefunction(env, mol, wfn, xtbData) - use xtb_setparam, only : gfn_method - use xtb_type_environment - use xtb_type_molecule - use xtb_type_wavefunction - use xtb_type_param - use xtb_disp_ncoord - use xtb_eeq - use xtb_scc_core - type(TMolecule), intent(in) :: mol - type(TWavefunction), intent(inout) :: wfn - type(TEnvironment), intent(inout) :: env - type(TxTBData), intent(in) :: xtbData - type(chrg_parameter) :: chrgeq - real(wp), allocatable :: cn(:) + do ii = 1, huge(ii) - 1 + if (rhs(ii) == c_null_char) then + exit + end if + end do + allocate(character(len=ii-1) :: lhs) + lhs = transfer(rhs(1:ii-1), lhs) - ! do an EEQ guess - allocate( cn(mol%n), source = 0.0_wp ) - call new_charge_model_2019(chrgeq,mol%n,mol%at) - call ncoord_erf(mol%n,mol%at,mol%xyz,cn) - call eeq_chrgeq(mol,env,chrgeq,cn,wfn%q) - deallocate(cn) +end subroutine c_f_character - call iniqshell(xtbData,mol%n,mol%at,mol%z,wfn%nshell,wfn%q,wfn%qsh,gfn_method) -end subroutine eeq_guess_wavefunction !> Cold fusion check -integer function verify_xtb_molecule(mol) result(status) - use xtb_type_molecule +integer function verifyMolecule(mol) result(status) type(TMolecule), intent(in) :: mol integer :: iat, jat status = 0 @@ -112,210 +57,16 @@ integer function verify_xtb_molecule(mol) result(status) if (mol%dist(jat, iat) < 1.0e-9_wp) status = status + 1 end do end do -end function verify_xtb_molecule - -logical function verify_xtb_wavefunction(basis, wfn) result(status) - use xtb_type_basisset - use xtb_type_wavefunction - type(TBasisset), intent(in) :: basis - type(TWavefunction), intent(in) :: wfn - status = wfn%n > 0 .and. wfn%nao > 0 .and. wfn%nshell > 0 & - & .and. wfn%n == basis%n .and. wfn%nshell == basis%nshell & - & .and. wfn%nao == basis%nao -end function verify_xtb_wavefunction - -!> optional return to a c_ptr in case it is not a null pointer and the -! Fortran value has been calculated -subroutine c_return_int_0d(c_array, f_array) - integer(c_int), intent(out), target :: c_array - integer, intent(in) :: f_array - if (c_associated(c_loc(c_array))) then - c_array = f_array - endif -end subroutine c_return_int_0d - -!> optional return to a c_ptr in case it is not a null pointer and the -! Fortran value has been calculated -subroutine c_return_double_0d(c_array, f_array) - real(c_double), intent(out), target :: c_array - real(wp), intent(in) :: f_array - if (c_associated(c_loc(c_array))) then - c_array = f_array - endif -end subroutine c_return_double_0d - -!> optional return to a c_ptr in case it is not a null pointer and the -! Fortran value has been calculated -subroutine c_return_double_0d_alloc(c_array, f_array) - real(c_double), intent(out), target :: c_array - real(wp), allocatable, intent(in) :: f_array - if (c_associated(c_loc(c_array)) .and. allocated(f_array)) then - c_array = f_array - endif -end subroutine c_return_double_0d_alloc - -!> optional return to a c_ptr in case it is not a null pointer and the -! Fortran array has been populated -subroutine c_return_double_1d(c_array, f_array) - real(c_double), intent(out), target :: c_array(:) - real(wp), intent(in) :: f_array(:) - if (c_associated(c_loc(c_array))) then - c_array = f_array - endif -end subroutine c_return_double_1d - -!> optional return to a c_ptr in case it is not a null pointer and the -! Fortran array has been populated -subroutine c_return_double_1d_alloc(c_array, f_array) - real(c_double), intent(out), target :: c_array(:) - real(wp), allocatable, intent(in) :: f_array(:) - if (c_associated(c_loc(c_array)) .and. allocated(f_array)) then - c_array = f_array - endif -end subroutine c_return_double_1d_alloc - -!> optional return to a c_ptr in case it is not a null pointer and the -! Fortran array has been populated (2D version) -subroutine c_return_double_2d(c_array, f_array) - real(c_double), intent(out), target :: c_array(:,:) - real(wp), intent(in) :: f_array(:,:) - if (c_associated(c_loc(c_array))) then - c_array = f_array - endif -end subroutine c_return_double_2d +end function verifyMolecule -!> optional return to a c_ptr in case it is not a null pointer and the -! Fortran array has been populated (2D version) -subroutine c_return_double_2d_alloc(c_array, f_array) - real(c_double), intent(out), target :: c_array(:,:) - real(wp), allocatable, intent(in) :: f_array(:,:) - if (c_associated(c_loc(c_array)) .and. allocated(f_array)) then - c_array = f_array - endif -end subroutine c_return_double_2d_alloc -!> convert vector of chars to deferred size character -subroutine c_string_convert(f_string, c_string) - character(c_char), dimension(*), intent(in) :: c_string - character(len=:), allocatable, intent(out) :: f_string - integer :: i - i = 0 - f_string = '' - do - i = i+1 - if (c_string(i).eq.c_null_char) exit - f_string = f_string//c_string(i) - enddo -end subroutine c_string_convert +subroutine checkGlobalEnv + use xtb_mctc_global, only : persistentEnv + if (.not.allocated(persistentEnv)) then + allocate(persistentEnv) + call init(persistentEnv) + end if +end subroutine checkGlobalEnv -subroutine c_get_double_0d_alloc(f_array, c_array) - real(c_double), intent(in), target :: c_array - real(wp), allocatable, intent(out) :: f_array - if (c_associated(c_loc(c_array))) then - f_array = c_array - endif -end subroutine c_get_double_0d_alloc - -subroutine c_get_double_1d_alloc(f_array, c_array) - real(c_double), intent(in), target :: c_array(:) - real(wp), allocatable, intent(out) :: f_array(:) - if (c_associated(c_loc(c_array))) then - f_array = c_array - endif -end subroutine c_get_double_1d_alloc - -subroutine c_get_double_2d_alloc(f_array, c_array) - real(c_double), intent(in), target :: c_array(:,:) - real(wp), allocatable, intent(out) :: f_array(:,:) - if (c_associated(c_loc(c_array))) then - f_array = c_array - endif -end subroutine c_get_double_2d_alloc - -subroutine c_get_double_0d(f_array, c_array, default) - real(c_double), intent(in), target :: c_array - real(wp), intent(out) :: f_array - real(wp), intent(in) :: default - if (c_associated(c_loc(c_array))) then - f_array = c_array - else - f_array = default - endif -end subroutine c_get_double_0d - -subroutine c_get_double_1d(f_array, c_array, default) - real(c_double), intent(in), target :: c_array(:) - real(wp), intent(out) :: f_array(:) - real(wp), intent(in) :: default - if (c_associated(c_loc(c_array))) then - f_array = c_array - else - f_array = default - endif -end subroutine c_get_double_1d - -subroutine c_get_double_2d(f_array, c_array, default) - real(c_double), intent(in), target :: c_array(:,:) - real(wp), intent(out) :: f_array(:,:) - real(wp), intent(in) :: default - if (c_associated(c_loc(c_array))) then - f_array = c_array - else - f_array = default - endif -end subroutine c_get_double_2d - -subroutine c_get_int_0d(f_array, c_array, default) - integer(c_int), intent(in), target :: c_array - integer, intent(out) :: f_array - integer, intent(in) :: default - if (c_associated(c_loc(c_array))) then - f_array = c_array - else - f_array = default - endif -end subroutine c_get_int_0d - -subroutine c_get_int_1d(f_array, c_array, default) - integer(c_int), intent(in), target :: c_array(:) - integer, intent(out) :: f_array(:) - integer, intent(in) :: default - if (c_associated(c_loc(c_array))) then - f_array = c_array - else - f_array = default - endif -end subroutine c_get_int_1d - -subroutine c_get_int_2d(f_array, c_array, default) - integer(c_int), intent(in), target :: c_array(:,:) - integer, intent(out) :: f_array(:,:) - integer, intent(in) :: default - if (c_associated(c_loc(c_array))) then - f_array = c_array - else - f_array = default - endif -end subroutine c_get_int_2d - -! open a unit for IO -subroutine c_open_file(unit, c_output, prlevel) - use xtb_mctc_io, only : stdout - integer, intent(out) :: unit - character(kind=c_char), intent(in) :: c_output(*) - integer, intent(in) :: prlevel - character(len=:), allocatable :: output - logical :: exist - - call c_string_convert(output, c_output) - - unit = stdout - if (output /= '-' .and. prlevel > 0) then - inquire(file=output, exist=exist) - if (exist) then - open(file=output, newunit=unit) - end if - endif -end subroutine c_open_file end module xtb_api_utils diff --git a/src/calculator.f90 b/src/calculator.f90 deleted file mode 100644 index f29864221..000000000 --- a/src/calculator.f90 +++ /dev/null @@ -1,83 +0,0 @@ -! This file is part of xtb. -! -! Copyright (C) 2017-2020 Stefan Grimme -! -! xtb 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. -! -! xtb 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 xtb. If not, see . -module xtb_calculators - use xtb_mctc_accuracy, only : wp - use xtb_type_environment - implicit none - private - - public :: gfn0_calculation, gfn1_calculation, gfn2_calculation - - interface - module subroutine gfn0_calculation & - (iunit,env,opt,mol,hl_gap,energy,gradient,stress,lattice_gradient) - use xtb_type_options - use xtb_type_molecule - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_param - use xtb_type_data - integer, intent(in) :: iunit - type(TMolecule), intent(inout) :: mol - type(peeq_options), intent(in) :: opt - type(TEnvironment), intent(inout) :: env - real(wp), intent(out) :: energy - real(wp), intent(out) :: hl_gap - real(wp), intent(out) :: gradient(3,mol%n) - real(wp), intent(out) :: stress(3,3) - real(wp), intent(out) :: lattice_gradient(3,3) - end subroutine gfn0_calculation - module subroutine gfn1_calculation & - (iunit,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) - use xtb_type_options - use xtb_type_molecule - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_param - use xtb_type_data - use xtb_type_pcem - integer, intent(in) :: iunit - type(TMolecule), intent(inout) :: mol - type(scc_options), intent(in) :: opt - type(TEnvironment), intent(inout) :: env - type(tb_pcem), intent(inout) :: pcem - type(TWavefunction),intent(inout) :: wfn - real(wp), intent(out) :: energy - real(wp), intent(out) :: hl_gap - real(wp), intent(out) :: gradient(3,mol%n) - end subroutine gfn1_calculation - module subroutine gfn2_calculation & - (iunit,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) - use xtb_type_options - use xtb_type_molecule - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_param - use xtb_type_data - use xtb_type_pcem - integer, intent(in) :: iunit - type(TMolecule), intent(inout) :: mol - type(TWavefunction),intent(inout) :: wfn - type(scc_options), intent(in) :: opt - type(TEnvironment), intent(inout) :: env - type(tb_pcem), intent(inout) :: pcem - real(wp), intent(out) :: energy - real(wp), intent(out) :: hl_gap - real(wp), intent(out) :: gradient(3,mol%n) - end subroutine gfn2_calculation - end interface -end module xtb_calculators diff --git a/src/disp/dftd4.f90 b/src/disp/dftd4.f90 index 4c27baf56..7e9c41772 100644 --- a/src/disp/dftd4.f90 +++ b/src/disp/dftd4.f90 @@ -86,8 +86,6 @@ subroutine d3init(nat,at) enddo ! integrate C6 coefficients - !$omp parallel default(none) private(i,j,ii,jj,ij,alpha,c6) shared(dispm) - !$omp do schedule(runtime) do i = 1, 86 do j = 1, i if (dispm%atoms(i) > 0 .and. dispm%atoms(j) > 0) then @@ -102,8 +100,6 @@ subroutine d3init(nat,at) endif enddo enddo - !$omp end do - !$omp end parallel end subroutine d3init diff --git a/src/filetools.F90 b/src/filetools.F90 index 9e74121b1..8beaa361b 100644 --- a/src/filetools.F90 +++ b/src/filetools.F90 @@ -15,251 +15,86 @@ ! You should have received a copy of the GNU Lesser General Public License ! along with xtb. If not, see . -#ifdef __GNUC__ -#if GCC_VERSION < 70500 -#define HAS_GCC_BUG_84412 -#endif -#endif - subroutine print_filelist(iunit) - use xtb_mctc_filetools + use xtb_mctc_global, only : persistentEnv use xtb_readin implicit none integer,intent(in) :: iunit - character(len=8) :: status - integer :: i - write(iunit,'(1x,"unit",2x,"open",3x,"action",5x,"filename")') - do i = 1, nfiles - select case(filelist(i)%status) - case('r','R'); status = "read " - case('a','A'); status = "appended" - case('w','W'); status = "written " - case('s','S'); status = "replaced" - case('d','D'); status = "deleted " - case('t','T'); status = "touched " - case default; status = " " - end select - write(iunit,'(i5,1x,a5,1x,":",1x,a8,3x,a)') & - abs(filelist(i)%unit), bool2string(filelist(i)%open), status, & - filelist(i)%name - enddo + call persistentEnv%io%list(iunit) + end subroutine print_filelist subroutine delete_file(name) - use xtb_mctc_filetools - use xtb_setparam + use xtb_mctc_global, only : persistentEnv implicit none character(len=*),intent(in) :: name - character(len=:),allocatable :: file - logical :: exist,opened - integer :: iunit,err - file = get_namespace(name) + call persistentEnv%io%deleteFile(name) - !$omp critical (io) - inquire(file=file, exist=exist, opened=opened, number=iunit) - if (opened) then - close(iunit,status='delete') - call pop_file(iunit,'d') - else if (exist) then - inquire(file=file, opened=opened, number=iunit) - open(newunit=iunit, file=file, iostat=err, status='old') - if (err.eq.0) then - close(iunit,status='delete') - call push_file(iunit,file,'d') - endif - endif - !$omp end critical (io) end subroutine delete_file subroutine touch_file(name) - use xtb_mctc_filetools + use xtb_mctc_global, only : persistentEnv use xtb_setparam implicit none character(len=*),intent(in) :: name - character(len=:),allocatable :: file - logical :: exist - integer :: iunit,err - file = get_namespace(name) + call persistentEnv%io%touchFile(name) - !$omp critical (io) - inquire(file=file, exist=exist) - if (.not.exist) then - open(newunit=iunit, file=file, iostat=err, status='new') - if (err.eq.0) then - close(iunit) - call push_file(iunit,file,'t') - endif - endif - !$omp end critical (io) end subroutine touch_file subroutine open_binary(iunit,name,status) - use xtb_mctc_filetools - use xtb_setparam + use xtb_mctc_global, only : persistentEnv implicit none + integer, intent(out) :: iunit character(len=*),intent(in) :: name character(len=1),intent(in) :: status - character(len=:),allocatable :: file - integer,intent(out) :: iunit - character(len=1) :: cstatus - integer :: err - logical :: exist - iunit = -1 ! failed to open file - cstatus = status - - file = get_namespace(name) - - ! acquire mutex for IO - !$omp critical (io) select case(status) case default - err = -1 - case('a','A') - inquire(file=file, exist=exist) - if (exist) then - open(newunit=iunit, file=file, iostat=err, action='write', status='old',& - form='unformatted') - if (err.ne.0) iunit = -1 - else - open(newunit=iunit, file=file, iostat=err, action='write', status='new',& - form='unformatted') - cstatus = 'w' - if (err.ne.0) iunit = -1 - endif + iunit = -1 case('r','R') - inquire(file=file, exist=exist) - if (.not.exist) then - file = name - inquire(file=file, exist=exist) - endif - if (exist) then - open(newunit=iunit, file=file, iostat=err, action='read', status='old',& - form='unformatted') - if (err.ne.0) iunit = -1 - endif + call persistentEnv%io%readBinary(iunit, name) case('w','W') - inquire(file=file, exist=exist) - if (exist) cstatus = 's' - open(newunit=iunit, file=file, iostat=err, action='write',& - form='unformatted') - if (err.ne.0) iunit = -1 + call persistentEnv%io%writeBinary(iunit, name) end select - if (iunit.ne.-1) then - call push_file(iunit,file,cstatus) - endif - !$omp end critical (io) - end subroutine open_binary subroutine open_file(iunit,name,status) - use xtb_mctc_filetools - use xtb_setparam + use xtb_mctc_global, only : persistentEnv implicit none character(len=*),intent(in) :: name character(len=1),intent(in) :: status - character(len=:),allocatable :: file integer,intent(out) :: iunit - character(len=1) :: cstatus - integer :: err - logical :: exist - iunit = -1 ! failed to open file - cstatus = status - - file = get_namespace(name) - - ! acquire mutex for IO - !$omp critical (io) select case(status) case default - err = -1 - case('a','A') - inquire(file=file, exist=exist) - if (exist) then - open(newunit=iunit, file=file, iostat=err, action='write', status='old') - if (err.ne.0) iunit = -1 - else - open(newunit=iunit, file=file, iostat=err, action='write', status='new') - cstatus = 'w' - if (err.ne.0) iunit = -1 - endif + iunit = -1 case('r','R') - inquire(file=file, exist=exist) - if (.not.exist) then - file = name - inquire(file=file, exist=exist) - endif - if (exist) then - open(newunit=iunit, file=file, iostat=err, action='read', status='old') - if (err.ne.0) iunit = -1 - endif + call persistentEnv%io%readFile(iunit, name) case('w','W') - inquire(file=file, exist=exist) - if (exist) cstatus = 's' - open(newunit=iunit, file=file, iostat=err, action='write') - if (err.ne.0) iunit = -1 + call persistentEnv%io%writeFile(iunit, name) end select - if (iunit.ne.-1) then - call push_file(iunit,file,cstatus) - endif - !$omp end critical (io) - end subroutine open_file subroutine close_file(unit) - use xtb_mctc_filetools + use xtb_mctc_global, only : persistentEnv implicit none integer,intent(in) :: unit logical :: opened - integer :: i - !$omp critical (io) - ! GCC7 cannot inquire on files with negative unit, fixed in GCC7.5.0 - ! see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=84412 -#ifdef HAS_GCC_BUG_84412 - opened = .false. - do i = 1, nfiles - if (unit == filelist(i)%unit) then - opened = opened .or. filelist(i)%open - end if - end do -#else - inquire(unit=unit,opened=opened) -#endif - if (opened) then - close(unit) - call pop_file(unit) - endif - !$omp end critical (io) + + call persistentEnv%io%closeFile(unit) + end subroutine close_file subroutine remove_file(unit) - use xtb_mctc_filetools + use xtb_mctc_global, only : persistentEnv implicit none integer,intent(in) :: unit - logical :: opened - integer :: i - !$omp critical (io) - ! GCC7 cannot inquire on files with negative unit, fixed in GCC7.5.0 - ! see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=84412 -#ifdef HAS_GCC_BUG_84412 - opened = .false. - do i = 1, nfiles - if (unit == filelist(i)%unit) then - opened = opened .or. filelist(i)%open - end if - end do -#else - inquire(unit=unit,opened=opened) -#endif - if (opened) then - close(unit,status='delete') - call pop_file(unit,'d') - endif - !$omp end critical (io) + + call persistentEnv%io%closeFile(unit, remove=.true.) + end subroutine remove_file diff --git a/src/gfn0_calculator.f90 b/src/gfn0_calculator.f90 deleted file mode 100644 index a24d13569..000000000 --- a/src/gfn0_calculator.f90 +++ /dev/null @@ -1,173 +0,0 @@ -! This file is part of xtb. -! -! Copyright (C) 2017-2019 Stefan Grimme -! -! xtb 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. -! -! xtb 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 xtb. If not, see . -submodule(xtb_calculators) gfn0_calc_implementation - implicit none -contains -! ======================================================================== -!> periodic GFN0-xTB (PEEQ) calculation -module subroutine gfn0_calculation & - (iunit,env,opt,mol,hl_gap,energy,gradient,stress,lattice_gradient) - use xtb_mctc_accuracy, only : wp - - use xtb_mctc_systools - - use xtb_type_options - use xtb_type_molecule - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_param - use xtb_type_data - - use xtb_setparam, only : gfn_method, ngrida - - use xtb_pbc_tools - use xtb_basis - use xtb_peeq - use xtb_solv_gbobc - use xtb_readparam - use xtb_paramset - - use xtb_main_setup - use xtb_xtb_calculator - use xtb_xtb_data - use xtb_xtb_gfn0 - - implicit none - - character(len=*), parameter :: source = 'calculator_gfn0' - - integer, intent(in) :: iunit - - type(TMolecule), intent(inout) :: mol - type(peeq_options), intent(in) :: opt - type(TEnvironment), intent(inout) :: env - - real(wp), intent(out) :: energy - real(wp), intent(out) :: hl_gap - real(wp), intent(out) :: gradient(3,mol%n) - real(wp), intent(out) :: stress(3,3) - real(wp), intent(out) :: lattice_gradient(3,3) - real(wp) :: sigma(3,3) - real(wp) :: inv_lat(3,3) - - integer, parameter :: wsc_rep(3) = [1,1,1] ! FIXME - - type(TWavefunction) :: wfn - type(TxTBCalculator) :: calc - type(scc_results) :: res - - character(len=*),parameter :: outfmt = & - '(9x,"::",1x,a,f24.12,1x,a,1x,"::")' - character(len=*), parameter :: p_fnv_gfn0 = 'param_gfn0-xtb.txt' - character(len=:), allocatable :: fnv - type(TxTBParameter) :: globpar - integer :: ipar,i - logical :: exist - logical :: exitRun - - logical :: okbas,diff - - gfn_method = 0 - sigma = 0.0_wp - - ! ==================================================================== - ! STEP 1: prepare geometry input - ! ==================================================================== - ! we assume that the user provides a resonable molecule input - ! -> all atoms are inside the unit cell, all data is set and consistent - - wfn%nel = nint(sum(mol%z) - mol%chrg) - wfn%nopen = mol%uhf - ! at this point, don't complain about odd multiplicities for even electron - ! systems and just fix it silently, the API is supposed catch this - if (mod(wfn%nopen,2) == 0.and.mod(wfn%nel,2) /= 0) wfn%nopen = 1 - if (mod(wfn%nopen,2) /= 0.and.mod(wfn%nel,2) == 0) wfn%nopen = 0 - - if (mol%npbc > 0) then - ! get Wigner-Seitz cell if necessary - ! -> this modifies the molecule, since the WSC is bound to a molecule - call generate_wsc(mol,mol%wsc,wsc_rep) - endif - - ! give an optional summary on the geometry used - if (opt%prlevel > 2) then - call main_geometry(iunit,mol) - endif - - ! ==================================================================== - ! STEP 2: get the parametrisation - ! ==================================================================== - ! we could require our user to perform this step, but if we want - ! to be sure about getting the correct parameters, we should do it here - ! let's check if we can find the parameter file - call rdpath(env%xtbpath, p_fnv_gfn0, fnv, exist) - if (exist) then - call newXTBCalculator(env, mol, calc, fnv) - else - call newXTBCalculator(env, mol, calc, p_fnv_gfn0) - end if - - if (opt%prlevel > 1) then - call gfn0_header(iunit) - endif - - lgbsa = len_trim(opt%solvent).gt.0 .and. opt%solvent.ne."none" & - & .and. mol%npbc == 0 ! GBSA is not yet periodic - if (lgbsa) then - call init_gbsa(iunit,trim(opt%solvent),0,opt%etemp,gfn_method,ngrida,& - & opt%prlevel > 0) - endif - - ! ==================================================================== - ! STEP 4: setup the initial wavefunction - ! ==================================================================== - - call wfn%allocate(mol%n,calc%basis%nshell,calc%basis%nao) - ! do a SAD guess since we are not need any of this information later - wfn%q = mol%chrg / real(mol%n,wp) - - ! ==================================================================== - ! STEP 5: do the calculation - ! ==================================================================== - - call peeq(env,mol,wfn,calc%basis,calc%xtbData,hl_gap,opt%etemp,opt%prlevel, & - & opt%grad,opt%ccm,opt%acc,energy,gradient,sigma,res) - - call env%check(exitRun) - if (exitRun) then - call env%error("Single point calculation terminated", source) - return - end if - - if (mol%npbc > 0) then - inv_lat = mat_inv_3x3(mol%lattice) - call sigma_to_latgrad(sigma,inv_lat,lattice_gradient) - stress = sigma/mol%volume - endif - - if (opt%prlevel > 0) then - write(iunit,'(9x,53(":"))') - write(iunit,outfmt) "total energy ", res%e_total,"Eh " - write(iunit,outfmt) "gradient norm ", res%gnorm, "Eh/α" - if (mol%npbc > 0) & - write(iunit,outfmt) "gradlatt norm ", norm2(lattice_gradient), "Eh/α" - write(iunit,outfmt) "HOMO-LUMO gap ", res%hl_gap, "eV " - write(iunit,'(9x,53(":"))') - endif - -end subroutine gfn0_calculation -end submodule gfn0_calc_implementation diff --git a/src/gfn1_calculator.f90 b/src/gfn1_calculator.f90 deleted file mode 100644 index 3e3b2bd2e..000000000 --- a/src/gfn1_calculator.f90 +++ /dev/null @@ -1,184 +0,0 @@ -! This file is part of xtb. -! -! Copyright (C) 2017-2019 Stefan Grimme -! -! xtb 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. -! -! xtb 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 xtb. If not, see . -submodule(xtb_calculators) gfn1_calc_implementation - implicit none -contains -! ======================================================================== -!> GFN1-xTB calculation -module subroutine gfn1_calculation & - (iunit,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) - use xtb_mctc_accuracy, only : wp - - use xtb_mctc_systools - - use xtb_type_options - use xtb_type_molecule - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_param - use xtb_type_data - use xtb_type_pcem - - use xtb_setparam, only : gfn_method, ngrida - - use xtb_basis - use xtb_eeq - use xtb_chargemodel - use xtb_disp_ncoord - use xtb_scc_core - use xtb_scf - use xtb_solv_gbobc - use xtb_embedding - use xtb_restart - use xtb_readparam - use xtb_paramset - - use xtb_main_setup - use xtb_xtb_calculator - use xtb_xtb_data - use xtb_xtb_gfn1 - - implicit none - - character(len=*), parameter :: source = 'calculator_gfn1' - - integer, intent(in) :: iunit - - type(TMolecule), intent(inout) :: mol - type(scc_options), intent(in) :: opt - type(TEnvironment), intent(inout) :: env - type(tb_pcem), intent(inout) :: pcem - type(TWavefunction),intent(inout) :: wfn - - real(wp), intent(out) :: energy - real(wp), intent(out) :: hl_gap - real(wp), intent(out) :: gradient(3,mol%n) - - integer, parameter :: wsc_rep(3) = [1,1,1] ! FIXME - - type(TxTBCalculator) :: calc - type(scc_results) :: res - type(chrg_parameter) :: chrgeq - - real(wp), allocatable :: cn(:) - - character(len=*),parameter :: outfmt = & - '(9x,"::",1x,a,f24.12,1x,a,1x,"::")' - character(len=*), parameter :: p_fnv_gfn1 = 'param_gfn1-xtb.txt' - character(len=:), allocatable :: fnv - type(TxTBParameter) :: globpar - integer :: ipar - logical :: exist - logical :: exitRun - - logical :: okbas,diff - - gfn_method = 1 - call init_pcem - - ! ==================================================================== - ! STEP 1: prepare geometry input - ! ==================================================================== - ! we assume that the user provides a resonable molecule input - ! -> all atoms are inside the unit cell, all data is set and consistent - - wfn%nel = nint(sum(mol%z) - mol%chrg) - wfn%nopen = mol%uhf - ! at this point, don't complain about odd multiplicities for even electron - ! systems and just fix it silently, the API is supposed catch this - if (mod(wfn%nopen,2) == 0.and.mod(wfn%nel,2) /= 0) wfn%nopen = 1 - if (mod(wfn%nopen,2) /= 0.and.mod(wfn%nel,2) == 0) wfn%nopen = 0 - - ! give an optional summary on the geometry used - if (opt%prlevel > 2) then - call main_geometry(iunit,mol) - endif - - ! ==================================================================== - ! STEP 2: get the parametrisation - ! ==================================================================== - ! we could require our user to perform this step, but if we want - ! to be sure about getting the correct parameters, we should do it here - call rdpath(env%xtbpath, p_fnv_gfn1, fnv, exist) - if (exist) then - call newXTBCalculator(env, mol, calc, fnv) - else - call newXTBCalculator(env, mol, calc, p_fnv_gfn1) - end if - - if (opt%prlevel > 1) then - call gfn1_header(iunit) - endif - - lgbsa = len_trim(opt%solvent).gt.0 .and. opt%solvent.ne."none" - if (lgbsa) then - call init_gbsa(iunit,trim(opt%solvent),0,opt%etemp,gfn_method,ngrida, & - & opt%prlevel > 0) - endif - - ! ==================================================================== - ! STEP 4: setup the initial wavefunction - ! ==================================================================== - - call wfn%allocate(mol%n,calc%basis%nshell,calc%basis%nao) - - if (mol%npbc == 0) then - ! do a EEQ guess - allocate( cn(mol%n), source = 0.0_wp ) - call new_charge_model_2019(chrgeq,mol%n,mol%at) - call ncoord_erf(mol%n,mol%at,mol%xyz,cn) - call eeq_chrgeq(mol,env,chrgeq,cn,wfn%q) - deallocate(cn) - call env%check(exitRun) - if (exitRun) then - call env%error("EEQ quess failed", source) - end if - else - wfn%q(:) = 0.0_wp - end if - - call iniqshell(calc%xtbData,mol%n,mol%at,mol%z,calc%basis%nshell,wfn%q,wfn%qsh,gfn_method) - - if (opt%restart) & - call readRestart(env,wfn,'xtbrestart',mol%n,mol%at,gfn_method,exist,.false.) - - ! ==================================================================== - ! STEP 5: do the calculation - ! ==================================================================== - call scf(env,mol,wfn,calc%basis,pcem,calc%xtbData,hl_gap, & - & opt%etemp,opt%maxiter,opt%prlevel,.false.,opt%grad,opt%acc, & - & energy,gradient,res) - - call env%check(exitRun) - if (exitRun) then - call env%error("SCF calculation terminated", source) - end if - - if (opt%restart) then - call writeRestart(env,wfn,'xtbrestart',gfn_method) - endif - - if (opt%prlevel > 0) then - write(iunit,'(9x,53(":"))') - write(iunit,outfmt) "total energy ", res%e_total,"Eh " - write(iunit,outfmt) "gradient norm ", res%gnorm, "Eh/α" - write(iunit,outfmt) "HOMO-LUMO gap ", res%hl_gap, "eV " - write(iunit,'(9x,53(":"))') - endif - -end subroutine gfn1_calculation -end submodule gfn1_calc_implementation diff --git a/src/gfn2_calculator.f90 b/src/gfn2_calculator.f90 deleted file mode 100644 index 2d5a7aa9e..000000000 --- a/src/gfn2_calculator.f90 +++ /dev/null @@ -1,180 +0,0 @@ -! This file is part of xtb. -! -! Copyright (C) 2017-2019 Stefan Grimme -! -! xtb 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. -! -! xtb 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 xtb. If not, see . -submodule(xtb_calculators) gfn2_calc_implementation - implicit none -contains -! ======================================================================== -!> GFN2-xTB calculation -module subroutine gfn2_calculation & - (iunit,env,opt,mol,pcem,wfn,hl_gap,energy,gradient) - use xtb_mctc_accuracy, only : wp - - use xtb_mctc_systools - - use xtb_type_options - use xtb_type_molecule - use xtb_type_wavefunction - use xtb_type_basisset - use xtb_type_param - use xtb_type_data - use xtb_type_pcem - - use xtb_setparam, only : gfn_method, ngrida - - use xtb_basis - use xtb_eeq - use xtb_chargemodel - use xtb_disp_ncoord - use xtb_scc_core - use xtb_scf - use xtb_solv_gbobc - use xtb_embedding - use xtb_restart - use xtb_readparam - use xtb_paramset - - use xtb_main_setup - use xtb_xtb_calculator - use xtb_xtb_data - use xtb_xtb_gfn2 - - implicit none - - character(len=*), parameter :: source = 'calculator_gfn2' - - integer, intent(in) :: iunit - - type(TMolecule), intent(inout) :: mol - type(TWavefunction),intent(inout) :: wfn - type(scc_options), intent(in) :: opt - type(TEnvironment), intent(inout) :: env - type(tb_pcem), intent(inout) :: pcem - - real(wp), intent(out) :: energy - real(wp), intent(out) :: hl_gap - real(wp), intent(out) :: gradient(3,mol%n) - - integer, parameter :: wsc_rep(3) = [1,1,1] ! FIXME - - type(TxTBCalculator) :: calc - type(scc_results) :: res - type(chrg_parameter) :: chrgeq - - real(wp), allocatable :: cn(:) - - character(len=*),parameter :: outfmt = & - '(9x,"::",1x,a,f24.12,1x,a,1x,"::")' - character(len=*), parameter :: p_fnv_gfn2 = 'param_gfn2-xtb.txt' - character(len=:), allocatable :: fnv - type(TxTBParameter) :: globpar - integer :: ipar - logical :: exist,diff - logical :: exitRun - - logical :: okbas - - gfn_method = 2 - call init_pcem - - ! ==================================================================== - ! STEP 1: prepare geometry input - ! ==================================================================== - ! we assume that the user provides a resonable molecule input - ! -> all atoms are inside the unit cell, all data is set and consistent - - wfn%nel = nint(sum(mol%z) - mol%chrg) - wfn%nopen = mol%uhf - ! at this point, don't complain about odd multiplicities for even electron - ! systems and just fix it silently, the API is supposed catch this - if (mod(wfn%nopen,2) == 0.and.mod(wfn%nel,2) /= 0) wfn%nopen = 1 - if (mod(wfn%nopen,2) /= 0.and.mod(wfn%nel,2) == 0) wfn%nopen = 0 - - ! give an optional summary on the geometry used - if (opt%prlevel > 2) then - call main_geometry(iunit,mol) - endif - - ! ==================================================================== - ! STEP 2: get the parametrisation - ! ==================================================================== - ! we could require our user to perform this step, but if we want - ! to be sure about getting the correct parameters, we should do it here - call rdpath(env%xtbpath, p_fnv_gfn2, fnv, exist) - if (exist) then - call newXTBCalculator(env, mol, calc, fnv) - else - call newXTBCalculator(env, mol, calc, p_fnv_gfn2) - end if - - if (opt%prlevel > 1) then - call gfn2_header(iunit) - endif - - lgbsa = len_trim(opt%solvent).gt.0 .and. opt%solvent.ne."none" - if (lgbsa) then - call init_gbsa(iunit,trim(opt%solvent),0,opt%etemp,gfn_method,ngrida, & - & opt%prlevel > 0) - endif - - ! ==================================================================== - ! STEP 4: setup the initial wavefunction - ! ==================================================================== - - call wfn%allocate(mol%n,calc%basis%nshell,calc%basis%nao) - - ! do an EEQ guess - allocate( cn(mol%n), source = 0.0_wp ) - call new_charge_model_2019(chrgeq,mol%n,mol%at) - call ncoord_erf(mol%n,mol%at,mol%xyz,cn) - call eeq_chrgeq(mol,env,chrgeq,cn,wfn%q) - deallocate(cn) - call env%check(exitRun) - if (exitRun) then - call env%error("EEQ quess failed", source) - end if - - call iniqshell(calc%xtbData,mol%n,mol%at,mol%z,calc%basis%nshell,wfn%q,wfn%qsh,gfn_method) - - if (opt%restart) & - call readRestart(env,wfn,'xtbrestart',mol%n,mol%at,gfn_method,exist,.false.) - - ! ==================================================================== - ! STEP 5: do the calculation - ! ==================================================================== - call scf(env,mol,wfn,calc%basis,pcem,calc%xtbData,hl_gap, & - & opt%etemp,opt%maxiter,opt%prlevel,.false.,opt%grad,opt%acc, & - & energy,gradient,res) - - call env%check(exitRun) - if (exitRun) then - call env%error("SCF calculation terminated", source) - end if - - if (opt%restart) then - call writeRestart(env,wfn,'xtbrestart',gfn_method) - endif - - if (opt%prlevel > 0) then - write(iunit,'(9x,53(":"))') - write(iunit,outfmt) "total energy ", res%e_total,"Eh " - write(iunit,outfmt) "gradient norm ", res%gnorm, "Eh/α" - write(iunit,outfmt) "HOMO-LUMO gap ", res%hl_gap, "eV " - write(iunit,'(9x,53(":"))') - endif - -end subroutine gfn2_calculation -end submodule gfn2_calc_implementation diff --git a/src/main/property.f90 b/src/main/property.f90 index ba0389bf8..82305eb54 100644 --- a/src/main/property.f90 +++ b/src/main/property.f90 @@ -72,6 +72,7 @@ end subroutine write_energy_gff subroutine main_property & (iunit,mol,wfx,basis,xtbData,res,acc) + use xtb_mctc_global, only : persistentEnv use xtb_mctc_convert !! ======================================================================== diff --git a/src/main/setup.f90 b/src/main/setup.f90 index 61db38f82..376152426 100644 --- a/src/main/setup.f90 +++ b/src/main/setup.f90 @@ -18,26 +18,30 @@ !> TODO module xtb_main_setup use xtb_mctc_accuracy, only : wp + use xtb_mctc_systools, only : rdpath use xtb_type_calculator, only : TCalculator use xtb_type_dummycalc, only : TDummyCalculator use xtb_type_environment, only : TEnvironment use xtb_type_molecule, only : TMolecule - use xtb_type_param, only : TxTBParameter + use xtb_type_param, only : TxTBParameter, chrg_parameter use xtb_type_wavefunction, only : TWavefunction use xtb_readparam, only : readParam use xtb_paramset, only : use_parameterset use xtb_basis, only : newBasisset use xtb_xtb_calculator, only : TxTBCalculator use xtb_gfnff_calculator, only : TGFFCalculator - use xtb_eeq, only : goedecker_chrgeq + use xtb_eeq, only : eeq_chrgeq use xtb_iniq, only : iniqcn use xtb_scc_core, only : iniqshell + use xtb_embedding, only : read_pcem use xtb_setparam use xtb_disp_ncoord + use xtb_chargemodel + use xtb_solv_gbobc, only : lgbsa, init_gbsa implicit none private - public :: newCalculator + public :: newCalculator, newWavefunction, addSolvationModel public :: newXTBCalculator, newGFFCalculator @@ -71,7 +75,7 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy) case(p_ext_eht, p_ext_xtb) allocate(xtb) - call newXTBCalculator(env, mol, xtb, fname) + call newXTBCalculator(env, mol, xtb, fname, gfn_method, accuracy) call env%check(exitRun) if (exitRun) then @@ -80,7 +84,6 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy) end if call move_alloc(xtb, calc) - calc%accuracy = accuracy case(p_ext_gfnff) allocate(gfnff) @@ -101,7 +104,7 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy) end subroutine newCalculator -subroutine newXTBCalculator(env, mol, calc, fname) +subroutine newXTBCalculator(env, mol, calc, fname, method, accuracy) character(len=*), parameter :: source = 'main_setup_newXTBCalculator' @@ -111,28 +114,71 @@ subroutine newXTBCalculator(env, mol, calc, fname) type(TxTBCalculator), intent(out) :: calc - character(len=*), intent(in) :: fname + integer, intent(in), optional :: method + + real(wp), intent(in), optional :: accuracy + + character(len=*), intent(in), optional :: fname + character(len=:), allocatable :: filename type(TxTBParameter) :: globpar integer :: ich logical :: exist, okbas logical :: exitRun + if (present(fname)) then + filename = fname + else + if (present(method)) then + select case(method) + case(0) + call rdpath(env%xtbpath, 'param_gfn0-xtb.txt', filename, exist) + if (.not.exist) filename = 'param_gfn0-xtb.txt' + case(1) + call rdpath(env%xtbpath, 'param_gfn1-xtb.txt', filename, exist) + if (.not.exist) filename = 'param_gfn1-xtb.txt' + case(2) + call rdpath(env%xtbpath, 'param_gfn2-xtb.txt', filename, exist) + if (.not.exist) filename = 'param_gfn2-xtb.txt' + end select + end if + end if + if (.not.allocated(filename)) then + call env%error("No parameter file or parametrisation info provided", source) + return + end if + + if (present(accuracy)) then + calc%accuracy = accuracy + else + calc%accuracy = 1.0_wp + end if + + calc%etemp = 300.0_wp + calc%maxiter = 250 + !> Obtain the parameter file allocate(calc%xtbData) - call open_file(ich, fname, 'r') + call open_file(ich, filename, 'r') exist = ich /= -1 if (exist) then call readParam(env, ich, globpar, calc%xtbData, .true.) call close_file(ich) else ! no parameter file, check if we have one compiled into the code - call use_parameterset(fname, globpar, calc%xtbData, exist) + call use_parameterset(filename, globpar, calc%xtbData, exist) if (.not.exist) then - call env%error('Parameter file '//fname//' not found!', source) + call env%error('Parameter file '//filename//' not found!', source) return end if endif + if (present(method)) then + if (method /= calc%xtbData%level) then + call env%error("Requested method does not match loaded method", source) + return + end if + end if + call env%check(exitRun) if (exitRun) then call env%error("Could not load parameters", source) @@ -147,6 +193,15 @@ subroutine newXTBCalculator(env, mol, calc, fname) return end if + !> check for external point charge field + if (allocated(pcem_file)) then + call open_file(ich, pcem_file, 'r') + if (ich /= -1) then + call read_pcem(ich, env, calc%pcem, calc%xtbData%coulomb) + call close_file(ich) + end if + end if + end subroutine newXTBCalculator subroutine newGFFCalculator(env, mol, calc, fname, restart, version) @@ -214,4 +269,85 @@ subroutine newGFFCalculator(env, mol, calc, fname, restart, version) end subroutine newGFFCalculator + +subroutine newWavefunction(env, mol, calc, wfn) + type(TEnvironment), intent(inout) :: env + type(TWavefunction), intent(inout) :: wfn + type(TxTBCalculator), intent(in) :: calc + type(TMolecule), intent(in) :: mol + real(wp), allocatable :: cn(:) + type(chrg_parameter) :: chrgeq + + allocate(cn(mol%n)) + call wfn%allocate(mol%n,calc%basis%nshell,calc%basis%nao) + + wfn%nel = idint(sum(mol%z)) - nint(mol%chrg) + wfn%nopen = mol%uhf + if(wfn%nopen == 0 .and. mod(wfn%nel,2) /= 0) wfn%nopen=1 + + if (mol%npbc > 0) then + wfn%q = mol%chrg/real(mol%n,wp) + else + if (guess_charges.eq.p_guess_gasteiger) then + call iniqcn(mol%n,wfn%nel,mol%at,mol%z,mol%xyz,nint(mol%chrg),1.0_wp, & + & wfn%q,cn,calc%xtbData%level,.true.) + else if (guess_charges.eq.p_guess_goedecker) then + call new_charge_model_2019(chrgeq,mol%n,mol%at) + call ncoord_erf(mol%n,mol%at,mol%xyz,cn) + call eeq_chrgeq(mol,env,chrgeq,cn,wfn%q) + else + wfn%q = mol%chrg/real(mol%n,wp) + end if + end if + + call iniqshell(calc%xtbData,mol%n,mol%at,mol%z,calc%basis%nshell, & + & wfn%q,wfn%qsh,calc%xtbData%level) +end subroutine newWavefunction + + +subroutine addSolvationModel(env, calc, solvent, state, temp, nang) + type(TEnvironment), intent(inout) :: env + class(TCalculator), intent(in) :: calc + character(len=*), intent(in) :: solvent + integer, intent(in), optional :: state + real(wp), intent(in), optional :: temp + integer, intent(in), optional :: nang + integer :: solvState + real(wp) :: temperature + integer :: gridSize + integer :: level + + if (present(state)) then + solvState = state + else + solvState = 0 + end if + + if (present(temp)) then + temperature = temp + else + temperature = 298.15_wp + end if + + if (present(nang)) then + gridSize = nang + else + gridSize = ngrida + end if + + level = 0 + select type(calc) + type is(TxTBCalculator) + level = calc%xtbData%level + end select + + + lgbsa = len_trim(solvent).gt.0 .and. solvent.ne."none" + if (lgbsa) then + call init_gbsa(env%unit, trim(solvent), solvState, temperature, level, & + & gridSize, .false.) + endif + +end subroutine addSolvationModel + end module xtb_main_setup diff --git a/src/mctc/mctc_global.f90 b/src/mctc/mctc_global.f90 index 64af7d171..313b1469c 100644 --- a/src/mctc/mctc_global.f90 +++ b/src/mctc/mctc_global.f90 @@ -29,7 +29,8 @@ module xtb_mctc_global ! sanity status of the mctc environment, beware if it goes insane logical :: good = .true. - type(TEnvironment) :: persistentEnv + !> Global environment + type(TEnvironment), allocatable :: persistentEnv type :: errormsg character(len=:),allocatable :: msg diff --git a/src/mctc/mctc_init.f90 b/src/mctc/mctc_init.f90 index 11527400d..17a185880 100644 --- a/src/mctc/mctc_init.f90 +++ b/src/mctc/mctc_init.f90 @@ -53,6 +53,7 @@ subroutine mctc_init(progname,ntimer,verbose) ! set this for xtb_mctc_global name = progname + allocate(persistentEnv) call init(persistentEnv) end subroutine mctc_init diff --git a/src/meson.build b/src/meson.build index 26d9fd0b0..d8c96b905 100644 --- a/src/meson.build +++ b/src/meson.build @@ -41,7 +41,6 @@ srcs += files( 'bias_path.f90', 'blowsy.f90', 'broyden.f90', - 'calculator.f90', 'charge_model.f90', 'coffee.f90', 'constrain_param.f90', @@ -73,10 +72,7 @@ srcs += files( 'geosum.f90', 'getname.f90', 'getsymnum.f90', - 'gfn0_calculator.f90', 'gfn0param.f90', - 'gfn1_calculator.f90', - 'gfn2_calculator.f90', 'gfn_paramset.f90', 'gfn_prparam.f90', 'grad_core.f90', diff --git a/src/printout.f90 b/src/printout.f90 index 278e1089a..51a1c0e4e 100644 --- a/src/printout.f90 +++ b/src/printout.f90 @@ -58,6 +58,7 @@ end subroutine writecosmofile subroutine setup_summary(iunit,n,fname,xcontrol,wfx,xrc,exist) use xtb_mctc_accuracy, only : wp + use xtb_mctc_global, only : persistentEnv use xtb_mctc_systools use xtb_type_wavefunction use xtb_setparam @@ -84,8 +85,8 @@ subroutine setup_summary(iunit,n,fname,xcontrol,wfx,xrc,exist) call rdvar('HOSTNAME',cdum,err) if (err.eq.0) & write(iunit,'(10x,a,":",1x,a)') 'hostname ',cdum - if (allocated(xenv%namespace)) & - write(iunit,'(10x,a,":",1x,a)') 'calculation namespace ',xenv%namespace + if (allocated(persistentEnv%io%namespace)) & + write(iunit,'(10x,a,":",1x,a)') 'calculation namespace ',persistentEnv%io%namespace ! ---------------------------------------------------------------------- ! print the home and path to check if there are set correctly write(iunit,'(10x,a,":",1x,a)') 'coordinate file ',fname diff --git a/src/prog/main.F90 b/src/prog/main.F90 index 5d5907b25..6382c7d6a 100644 --- a/src/prog/main.F90 +++ b/src/prog/main.F90 @@ -366,6 +366,7 @@ subroutine xtbMain(env, argParser) endif mol%chrg = real(chrg, wp) + mol%uhf = nalphabeta wfn%nel = idint(sum(mol%z)) - chrg wfn%nopen = nalphabeta if(wfn%nopen == 0 .and. mod(wfn%nel,2) /= 0) wfn%nopen=1 @@ -1030,6 +1031,7 @@ end subroutine xtbMain !> Parse command line arguments and forward them to settings subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, & & restart, gsolvstate, strict, copycontrol, coffee) + use xtb_mctc_global, only : persistentEnv !> Name of error producer character(len=*), parameter :: source = "prog_main_parseArguments" @@ -1157,8 +1159,8 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, & end if case('--namespace') - call args%nextArg(xenv%namespace) - if (.not.allocated(xenv%namespace)) then + call args%nextArg(persistentEnv%io%namespace) + if (.not.allocated(persistentEnv%io%namespace)) then call env%error("Namespace argument is missing", source) end if diff --git a/src/prog/topology.f90 b/src/prog/topology.f90 index 53622aa02..4f33221dc 100644 --- a/src/prog/topology.f90 +++ b/src/prog/topology.f90 @@ -132,6 +132,7 @@ end subroutine xtbTopology !> Parse command line arguments subroutine parseArguments(env, args, param) + use xtb_mctc_global, only : persistentEnv !> Name of error producer character(len=*), parameter :: source = "prog_topology_parseArguments" @@ -177,8 +178,8 @@ subroutine parseArguments(env, args, param) end if case('--namespace') - call args%nextArg(xenv%namespace) - if (.not.allocated(xenv%namespace)) then + call args%nextArg(persistentEnv%io%namespace) + if (.not.allocated(persistentEnv%io%namespace)) then call env%error("Namespace argument is missing", source) end if diff --git a/src/scf_module.f90 b/src/scf_module.f90 index 7d57f10f8..bd90e3b5f 100644 --- a/src/scf_module.f90 +++ b/src/scf_module.f90 @@ -311,7 +311,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & ! initialize the fgb matrix (dielectric screening of the Coulomb potential) call compute_fgb(gbsa,fgb,fhb) ! initialize the CM5 charges computation - if(gfn_method > 1) then !GFN2 does not use CM5 charges + if(xtbData%level > 1) then !GFN2 does not use CM5 charges cm5=wfn%q cm5a=0.d0 dcm5a=0.d0 @@ -331,7 +331,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & allocate(selfEnergy(maxval(xtbData%nshell), mol%n)) allocate(dSEdcn(maxval(xtbData%nshell), mol%n)) - call setzshell(xtbData,mol%n,mol%at,basis%nshell,mol%z,zsh,eatoms,gfn_method) + call setzshell(xtbData,mol%n,mol%at,basis%nshell,mol%z,zsh,eatoms,xtbData%level) ! fill levels if(wfn%nel.ne.0) then @@ -414,7 +414,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & idnum(ii) = mol%at(jat) end do - if(gfn_method == 1)then + if(xtbData%level == 1)then gav = gamAverage%harmonic else !GFN2 gav = gamAverage%arithmetic @@ -434,7 +434,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & allocate(Vpc(basis%nshell)) vpc(:) = 0.0_wp if(lpcem)then - if (gfn_method == 1)then + if (xtbData%level == 1)then call jpot_pcem_gfn1(xtbData%coulomb,mol%n,pcem,xtbData%nshell,mol%at, & & mol%xyz,xtbData%coulomb%gExp,Vpc) else ! GFN2 @@ -454,9 +454,9 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & if (allocated(xtbData%halogen)) & write(env%unit,intfmt) "# halogen bonds ",nxb write(env%unit,intfmt) "max. iterations ",maxiter - if (gfn_method == 2) & + if (xtbData%level == 2) & write(env%unit,chrfmt) "Hamiltonian ","GFN2-xTB" - if (gfn_method == 1) & + if (xtbData%level == 1) & write(env%unit,chrfmt) "Hamiltonian ","GFN1-xTB" write(env%unit,chrfmt) "restarted? ",bool2string(restart) write(env%unit,chrfmt) "GBSA solvation ",bool2string(lgbsa) @@ -506,7 +506,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & ! ------------------------------------------------------------------------ ! Coordination number - if (gfn_method == 1) then + if (xtbData%level == 1) then ! D3 part first because we need CN call getCoordinationNumber(mol, trans, 40.0_wp, cnType%exp, cn, dcndr, dcndL) else @@ -517,7 +517,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & ! ------------------------------------------------------------------------ ! dispersion (DFT-D type correction) call latp%getLatticePoints(trans, 60.0_wp) - if (gfn_method == 1) then + if (xtbData%level == 1) then call d3_gradient & & (mol, trans, xtbData%dispersion%dpar, 4.0_wp, 60.0_wp, & & cn, dcndr, dcndL, ed, gradient, sigma) @@ -723,7 +723,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & ! ------------------------------------------------------------------------ ! Solvation contributions from GBSA if (lgbsa) then - if (gfn_method > 1) then + if (xtbData%level > 1) then call compute_gb_egrad(gbsa, wfn%q, gborn, ghb, gradient, minpr) else cm5(:)=wfn%q+cm5a @@ -749,7 +749,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, & ! ------------------------------------------------------------------------ ! ES point charge embedding if (lpcem) then - if (gfn_method == 1) then + if (xtbData%level == 1) then call pcem_grad_gfn1(xtbData%coulomb,gradient,pcem%grd,mol%n,pcem,mol%at, & & xtbData%nshell,mol%xyz,xtbData%coulomb%gExp,wfn%qsh) else diff --git a/src/setparam.f90 b/src/setparam.f90 index 0641673d8..d980f28ef 100644 --- a/src/setparam.f90 +++ b/src/setparam.f90 @@ -430,6 +430,7 @@ subroutine initrand end subroutine initrand function get_namespace(string) result(name) + use xtb_mctc_global, only : persistentEnv implicit none character(len=*),intent(in) :: string character(len=:),allocatable :: name @@ -437,11 +438,11 @@ function get_namespace(string) result(name) name = string return endif - if (allocated(xenv%namespace)) then + if (allocated(persistentEnv%io%namespace)) then if (string(1:1).eq.dot) then - name = dot//xenv%namespace//string + name = dot//persistentEnv%io%namespace//string else - name = xenv%namespace//dot//string + name = persistentEnv%io%namespace//dot//string endif else name = string diff --git a/src/type/CMakeLists.txt b/src/type/CMakeLists.txt index f0a600f26..7b5204a17 100644 --- a/src/type/CMakeLists.txt +++ b/src/type/CMakeLists.txt @@ -30,6 +30,7 @@ list(APPEND srcs "${dir}/environment.f90" "${dir}/fragments.f90" "${dir}/identitymap.f90" + "${dir}/iohandler.f90" "${dir}/molecule.f90" "${dir}/neighbourlist.f90" "${dir}/latticepoint.f90" diff --git a/src/type/environment.f90 b/src/type/environment.f90 index bfc5c509c..dad42ca0b 100644 --- a/src/type/environment.f90 +++ b/src/type/environment.f90 @@ -20,6 +20,7 @@ module xtb_type_environment use xtb_mctc_accuracy, only : wp use xtb_mctc_io, only : stdout use xtb_mctc_systools, only : rdvar, rdarg + use xtb_type_iohandler, only : TIOHandler, init_ => init implicit none private @@ -66,6 +67,9 @@ module xtb_type_environment !> Calculations home directory character(len=:),allocatable :: xtbhome + !> File input/output handling + type(TIOHandler) :: io + contains !> Add a warning to the message log @@ -142,6 +146,8 @@ subroutine initEnvironment(self, strict) self%strict = .false. end if + call init_(self%io) + end subroutine initEnvironment diff --git a/src/type/iohandler.f90 b/src/type/iohandler.f90 new file mode 100644 index 000000000..2a9694f07 --- /dev/null +++ b/src/type/iohandler.f90 @@ -0,0 +1,470 @@ +! This file is part of xtb. +! +! Copyright (C) 2019-2020 Sebastian Ehlert +! +! xtb 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. +! +! xtb 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 xtb. If not, see . + +!> +module xtb_type_iohandler + use xtb_mctc_filetypes, only : generateFileMetaInfo + implicit none + private + + public :: TIOHandler, init + + + type :: TFileStatusEnum + integer :: readin = 1 + integer :: append = 2 + integer :: written = 3 + integer :: replaced = 4 + integer :: deleted = 5 + integer :: created = 6 + end type + + type(TFileStatusEnum), parameter :: fileStatus = TFileStatusEnum() + + + type :: TFileHandle + character(len=:),allocatable :: name + integer :: status + integer :: unit + logical :: open + end type TFileHandle + + interface TFileHandle + module procedure :: initFileHandle + end interface TFileHandle + + + type :: TIOHandler + + character(len=:), allocatable :: namespace + + integer :: count + type(TFileHandle), allocatable :: log(:) + + contains + + procedure :: readFile + procedure :: writeFile + procedure :: readBinary + procedure :: writeBinary + procedure :: touchFile + procedure :: closeFile + procedure :: deleteFile + procedure :: list + + procedure, private :: getName + procedure, private :: pushBack + procedure, private :: findUnit + + end type TIOHandler + + + interface init + module procedure :: initIOHandler + end interface init + + +contains + + +subroutine initIOHandler(self, namespace) + type(TIOHandler), intent(out) :: self + character(len=*), intent(in), optional :: namespace + + if (present(namespace)) then + self%namespace = namespace + end if + self%count = 0 + allocate(self%log(20)) + +end subroutine initIOHandler + + +elemental function initFileHandle(name, status, unit, open) result(self) + character(len=*), intent(in) :: name + integer, intent(in) :: status + integer, intent(in) :: unit + logical, intent(in) :: open + type(TFileHandle) :: self + + self%name = name + self%status = status + self%unit = unit + self%open = open + +end function initFileHandle + + +subroutine getName(self, file, filename) + class(TIOHandler), intent(inout) :: self + character(len=*), intent(in) :: file + character(len=:), allocatable, intent(out) :: filename + + if (allocated(self%namespace)) then + if (index(file, '/') == 0) then + if (file(1:1) == '.') then + filename = '.'//self%namespace//file + else + filename = self%namespace//'.'//file + end if + else + filename = file + end if + else + filename = file + end if + +end subroutine getName + + +subroutine pushBack(self, fileHandle) + class(TIOHandler), intent(inout) :: self + type(TFileHandle), intent(in) :: fileHandle + type(TFileHandle), allocatable :: tmp(:) + integer :: n + + self%count = self%count + 1 + if (self%count > size(self%log)) then + n = size(self%log) + call move_alloc(self%log, tmp) + allocate(self%log(n + n/2 + 1)) + self%log(1:n) = tmp + deallocate(tmp) + end if + self%log(self%count) = fileHandle + +end subroutine pushBack + + +subroutine findUnit(self, unit, pos) + class(TIOHandler), intent(in) :: self + integer, intent(in) :: unit + integer, intent(out) :: pos + integer :: i + + pos = 0 + do i = 1, self%count + if (self%log(i)%open .and. self%log(i)%unit == unit) then + pos = unit + exit + end if + end do + +end subroutine findUnit + + +subroutine list(self, unit) + class(TIOHandler), intent(in) :: self + integer, intent(in) :: unit + character(len=8) :: status + integer :: i + + if (self%count > 0) then + write(unit,'(1x,"unit",2x,"open",3x,"action",5x,"filename")') + end if + do i = 1, self%count + select case(self%log(i)%status) + case(fileStatus%readin); status = "read " + case(fileStatus%written); status = "written " + case(fileStatus%deleted); status = "deleted " + case(fileStatus%replaced); status = "replaced" + case(fileStatus%created); status = "created " + case default; status = "unknown " + end select + if (self%log(i)%open) then + write(unit,'(i5,1x,a5,1x,":",1x,a8,3x,a)') & + & abs(self%log(i)%unit), "true", status, self%log(i)%name + else + write(unit,'(i5,1x,a5,1x,":",1x,a8,3x,a)') & + & abs(self%log(i)%unit), "false", status, self%log(i)%name + end if + end do + +end subroutine list + + +subroutine readFile(self, unit, file, iostat) + class(TIOHandler), intent(inout) :: self + integer, intent(out) :: unit + character(len=*), intent(in) :: file + integer, intent(out), optional :: iostat + character(len=:), allocatable :: name + logical :: exist + integer :: error + + unit = -1 + error = 0 + + call self%getName(file, name) + !$omp critical(io) + inquire(file=name, exist=exist) + + if (exist) then + open(file=name, newunit=unit, status='old', action='read', iostat=error) + if (error == 0) then + call self%pushBack(TFileHandle(name, fileStatus%readin, unit, .true.)) + else + unit = -1 + end if + else + inquire(file=file, exist=exist) + if (exist) then + open(file=file, newunit=unit, status='old', action='read', iostat=error) + if (error == 0) then + call self%pushBack(TFileHandle(file, fileStatus%readin, unit, .true.)) + else + unit = -1 + end if + else + error = 1 + end if + end if + !$omp end critical(io) + + if (present(iostat)) then + iostat = error + end if + +end subroutine readFile + + +subroutine writeFile(self, unit, file, iostat) + class(TIOHandler), intent(inout) :: self + integer, intent(out) :: unit + character(len=*), intent(in) :: file + integer, intent(out), optional :: iostat + character(len=:), allocatable :: name + logical :: exist + integer :: error + + unit = -1 + error = 0 + + call self%getName(file, name) + !$omp critical(io) + inquire(file=name, exist=exist) + + open(file=name, newunit=unit, action='write', iostat=error) + if (error == 0) then + if (exist) then + call self%pushBack(TFileHandle(name, fileStatus%replaced, unit, .true.)) + else + call self%pushBack(TFileHandle(name, fileStatus%written, unit, .true.)) + end if + else + unit = -1 + end if + !$omp end critical(io) + + if (present(iostat)) then + iostat = error + end if + +end subroutine writeFile + + +subroutine readBinary(self, unit, file, iostat) + class(TIOHandler), intent(inout) :: self + integer, intent(out) :: unit + character(len=*), intent(in) :: file + integer, intent(out), optional :: iostat + character(len=:), allocatable :: name + logical :: exist + integer :: error + + unit = -1 + error = 0 + + call self%getName(file, name) + !$omp critical(io) + inquire(file=name, exist=exist) + + if (exist) then + open(file=name, newunit=unit, status='old', action='read', & + & form='unformatted', iostat=error) + if (error == 0) then + call self%pushBack(TFileHandle(name, fileStatus%readin, unit, .true.)) + else + unit = -1 + end if + else + inquire(file=file, exist=exist) + if (exist) then + open(file=file, newunit=unit, status='old', action='read', & + & form='unformatted', iostat=error) + if (error == 0) then + call self%pushBack(TFileHandle(file, fileStatus%readin, unit, .true.)) + else + unit = -1 + end if + else + error = 1 + end if + end if + !$omp end critical(io) + + if (present(iostat)) then + iostat = error + end if + +end subroutine readBinary + + +subroutine writeBinary(self, unit, file, iostat) + class(TIOHandler), intent(inout) :: self + integer, intent(out) :: unit + character(len=*), intent(in) :: file + integer, intent(out), optional :: iostat + character(len=:), allocatable :: name + logical :: exist + integer :: error + + unit = -1 + error = 0 + + call self%getName(file, name) + !$omp critical(io) + inquire(file=name, exist=exist) + + open(file=name, newunit=unit, action='write', form='unformatted', iostat=error) + if (error == 0) then + if (exist) then + call self%pushBack(TFileHandle(name, fileStatus%replaced, unit, .true.)) + else + call self%pushBack(TFileHandle(name, fileStatus%written, unit, .true.)) + end if + else + unit = -1 + end if + !$omp end critical(io) + + if (present(iostat)) then + iostat = error + end if + +end subroutine writeBinary + + +subroutine touchFile(self, file, iostat) + class(TIOHandler), intent(inout) :: self + character(len=*), intent(in) :: file + integer, intent(out), optional :: iostat + character(len=:), allocatable :: name + integer :: unit, error + logical :: exist + + error = 0 + + call self%getName(file, name) + !$omp critical(io) + inquire(file=name, exist=exist) + if (exist) then + error = 1 + else + open(file=name, newunit=unit, status='new', iostat=error) + if (error == 0) then + call self%pushBack(TFileHandle(name, fileStatus%created, unit, .false.)) + close(unit, iostat=error) + end if + end if + !$omp end critical(io) + + if (present(iostat)) then + iostat = error + end if + +end subroutine touchFile + + +subroutine closeFile(self, unit, iostat, remove) + class(TIOHandler), intent(inout) :: self + integer, intent(in) :: unit + integer, intent(out), optional :: iostat + logical, intent(in), optional :: remove + integer :: pos, error + logical :: delete + + if (present(remove)) then + delete = remove + else + delete = .false. + end if + + error = 0 + + call self%findUnit(unit, pos) + + !$omp critical(io) + if (pos > 0) then + if (delete) then + close(unit, iostat=error, status='delete') + else + close(unit, iostat=error) + end if + if (error == 0) then + self%log(pos)%open = .false. + if (delete) then + self%log(pos)%status = fileStatus%deleted + end if + end if + else + error = 1 + end if + !$omp end critical(io) + + if (present(iostat)) then + iostat = error + end if + +end subroutine closeFile + + +subroutine deleteFile(self, file, iostat) + class(TIOHandler), intent(inout) :: self + character(len=*), intent(in) :: file + integer, intent(out), optional :: iostat + character(len=:), allocatable :: name + integer :: unit + logical :: exist + integer :: error + + unit = -1 + error = 0 + + call self%getName(file, name) + !$omp critical(io) + inquire(file=name, exist=exist) + + if (exist) then + open(file=name, newunit=unit, status='old', iostat=error) + if (error == 0) then + call self%pushBack(TFileHandle(name, fileStatus%deleted, unit, .false.)) + close(unit, status='delete', iostat=error) + end if + else + error = 1 + end if + !$omp end critical(io) + + if (present(iostat)) then + iostat = error + end if + +end subroutine deleteFile + + +end module xtb_type_iohandler diff --git a/src/type/meson.build b/src/type/meson.build index 4780cc8cf..533a36bde 100644 --- a/src/type/meson.build +++ b/src/type/meson.build @@ -28,6 +28,7 @@ srcs += files( 'environment.f90', 'fragments.f90', 'identitymap.f90', + 'iohandler.f90', 'molecule.f90', 'neighbourlist.f90', 'latticepoint.f90', diff --git a/src/xtb/calculator.f90 b/src/xtb/calculator.f90 index c9501cfbf..01b2c6a74 100644 --- a/src/xtb/calculator.f90 +++ b/src/xtb/calculator.f90 @@ -59,6 +59,9 @@ module xtb_xtb_calculator !> Maximum number of cycles for SCC convergence integer :: maxiter + !> External potential + type(tb_pcem) :: pcem + contains !> Perform xTB single point calculation @@ -115,7 +118,6 @@ subroutine singlepoint(self, env, mol, wfn, printlevel, restart, & !> Detailed results type(scc_results), intent(out) :: results - type(tb_pcem) :: pcem integer :: i,ich integer :: mode_sp_run = 1 real(wp) :: efix @@ -132,21 +134,11 @@ subroutine singlepoint(self, env, mol, wfn, printlevel, restart, & hlgap = 0.0_wp efix = 0.0_wp - ! ------------------------------------------------------------------------ - !> external constrains which can be applied beforehand - - !> check for external point charge field - call open_file(ich, pcem_file, 'r') - if (ich /= -1) then - call read_pcem(ich, env, pcem, self%xtbData%coulomb) - call close_file(ich) - end if - ! ------------------------------------------------------------------------ ! actual calculation select case(self%xtbData%level) case(1, 2) - call scf(env,mol,wfn,self%basis,pcem,self%xtbData, & + call scf(env,mol,wfn,self%basis,self%pcem,self%xtbData, & & hlgap,self%etemp,self%maxiter,printlevel,restart,.true., & & self%accuracy,energy,gradient,results) @@ -166,10 +158,10 @@ subroutine singlepoint(self, env, mol, wfn, printlevel, restart, & ! post processing of gradient and energy ! point charge embedding gradient file - if (pcem%n > 0) then + if (allocated(pcem_grad) .and. self%pcem%n > 0) then call open_file(ich,pcem_grad,'w') - do i=1,pcem%n - write(ich,'(3f12.8)')pcem%grd(1:3,i) + do i=1,self%pcem%n + write(ich,'(3f12.8)')self%pcem%grd(1:3,i) enddo call close_file(ich) endif @@ -215,9 +207,9 @@ subroutine singlepoint(self, env, mol, wfn, printlevel, restart, & write(env%unit,outfmt) "LUMO orbital eigv.", wfn%emo(wfn%ihomo+1),"eV " endif write(env%unit,'(9x,"::",49("."),"::")') - if (gfn_method.eq.2) call print_gfn2_results(env%unit,results,verbose,lgbsa) - if (gfn_method.eq.1) call print_gfn1_results(env%unit,results,verbose,lgbsa) - if (gfn_method.eq.0) call print_gfn0_results(env%unit,results,verbose,lgbsa) + if (self%xtbData%level.eq.2) call print_gfn2_results(env%unit,results,verbose,lgbsa) + if (self%xtbData%level.eq.1) call print_gfn1_results(env%unit,results,verbose,lgbsa) + if (self%xtbData%level.eq.0) call print_gfn0_results(env%unit,results,verbose,lgbsa) write(env%unit,outfmt) "add. restraining ", efix, "Eh " if (verbose) then write(env%unit,'(9x,"::",49("."),"::")') From ccff38226f0e56ea7803e0612185d279b6c03c37 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Fri, 22 May 2020 16:59:04 +0200 Subject: [PATCH 2/4] Also update the README to reflect the change in the API --- python/README.md | 135 +++++++++++++++++++++++------------------------ 1 file changed, 66 insertions(+), 69 deletions(-) diff --git a/python/README.md b/python/README.md index a44f201af..5fe48a26d 100644 --- a/python/README.md +++ b/python/README.md @@ -1,91 +1,88 @@ # Python Wrapper for the Extended Tight Binding Program -This is the Python-side of the wrapper around the `xtb` library. +This is a basic `ctypes` Python API around the `xtb` library. This wrapper provides access to the three Hamiltonians (GFN0-xTB, GFN1-xTB -and GFN2-xTB) implemented in the `xtb` program, making accessible -energies, gradients and related properties like dipole moments, partial charges -and bond orders. +and GFN2-xTB and experimental access to the GFN-FF) implemented in the `xtb` +program, making accessible energies, gradients and related properties like +dipole moments, partial charges and bond orders. It is not yet a full replacement for the `xtb` binary since it is mainly intended to supplement its functionality. ## Usage -For the everyday use we recommend to use `xtb` together with the -Atomic Simulation Model (ASE) like +To use the C-API import the wrapper classes around the `xtb` objects: ```python -from ase.io import read -from xtb import GFN2 -mol = read('coord') -mol.set_calculator(GFN2()) -energy = mol.get_potential_energy() -forces = mol.get_forces() -charges = mol.get_charges() +from xtb.interface import Calculator, Results, Param +import numpy as np +# atomic numbers range from 1 to 86 +numbers = np.array( + [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], +) +# cartesian coordinates in atomic units (Bohr) +positions = np.array( + [ + [ 2.02799738646442, 0.09231312124713, -0.14310895950963], + [ 4.75011007621000, 0.02373496014051, -0.14324124033844], + [ 6.33434307654413, 2.07098865582721, -0.14235306905930], + [ 8.72860718071825, 1.38002919517619, -0.14265542523943], + [ 8.65318821103610, -1.19324866489847, -0.14231527453678], + [ 6.23857175648671, -2.08353643730276, -0.14218299370797], + [ 5.63266886875962, -4.69950321056008, -0.13940509630299], + [ 3.44931709749015, -5.48092386085491, -0.14318454855466], + [ 7.77508917214346, -6.24427872938674, -0.13107140408805], + [10.30229550927022, -5.39739796609292, -0.13672168520430], + [12.07410272485492, -6.91573621641911, -0.13666499342053], + [10.70038521493902, -2.79078533715849, -0.14148379504141], + [13.24597858727017, -1.76969072232377, -0.14218299370797], + [ 7.40891694074004, -8.95905928176407, -0.11636933482904], + [ 1.38702118184179, 2.05575746325296, -0.14178615122154], + [ 1.34622199478497, -0.86356704498496, 1.55590600570783], + [ 1.34624089204623, -0.86133716815647, -1.84340893849267], + [ 5.65596919189118, 4.00172183859480, -0.14131371969009], + [14.67430918222276, -3.26230980007732, -0.14344911021228], + [13.50897177220290, -0.60815166181684, 1.54898960808727], + [13.50780014200488, -0.60614855212345, -1.83214617078268], + [ 5.41408424778406, -9.49239668625902, -0.11022772492007], + [ 8.31919801555568, -9.74947502841788, 1.56539243085954], + [ 8.31511620712388, -9.76854236502758, -1.79108242206824], + ], +) + +calc = Calculator(Param.GFN2xTB, numbers, positions) +results = calc.singlepoint() + +print(results.get_energy()) ``` -The three calculators implement the corresponding Hamiltonians with the -following properties - -| |`GFN0`|`GFN1`|`GFN2`| -|--------------------|------|------|------| -| energy | x | x | x | -| forces | x | x | x | -| stress tensor | x | | | -| dipole moment | | x | x | -| partial charges | | x | x | -| atomic dipoles | | | x | -| atomic quadrupoles | | | x | -| bond order | | x | x | - -The atomic dipoles and quadrupoles as well as the bond orders are currently -not accessible from the atoms object and must be taken directly from the -calculator object. - -Additionally there is a solvation model calculator based on a -generalized Born model with a solvent accessible surface area, `GBSA`, -which can calculate Born radii and the surface area by itself -and together with an extended tight-binding calculator also a -solvation free energy. Use it as +The Python classes provide pretty much a one to one translation from the +C-API to Python, with the exception that the calculation environment and +molecular structure data are absorbed into the single point calculator +for convenience. + +To update the `xtb` objects with new coordinates or lattice parameters use ```python -from ase.io import read -from xtb import GFN2 -from xtb.solvation import GBSA -mol = read('coord') -mol.set_calculator(GBSA(solvent='ch2cl2', calc=GFN2())) -gsolv = mol.get_potential_energy() +calc.update(positions) +results = calc.singlepoint(results) ``` -which evaluates the GBSA model and performs two GFN2-xTB calculations. - -## Technical Details - -We provide a multilayered API for `xtb` starting from the Fortran -side with a set of standalone calculator subroutines in `xtb/calculator.f90`, -these calculators are exposed by using the `iso_c_binding` module as C-API -(see `xtb/c_api.f90` and `include/xtb.h`). -All non-Fortran compatible languages (that makes all languages except for Fortran), -should use this C-API to interface with `xtb`. +Previously calculated results objects can be used to restart single point +calculations. +To change the molecular structure data the calculator has to be reconstructed, +which means the number of atoms and the atomic types are immutable for an +already existing calculator. -To use the C-API in Python we decided to define the interface on the Python -side rather than on the Fortran side using the `ctypes` module. -This interface (or Python-API) hides some of the implementation dependent details -from the user and provides access to all main calculators implemented in `xtb`. +## Troubleshooting -Using the Python-API I wrapped everything up for the Atomic Simulation Environment -(ASE) by implementing a Calculator that handles the conversion between the -Atoms objects and the bundle of ndarrays required for Python-API. +Importing the shared library via `ctypes` is somewhat errorprone, in case the +default setup of the library breaks use -As a Python user I recommend using the ASE Calculators gives access to a -feature-rich environment for computational chemistry. -For more information on ASE I refer to their [detailed documentation](https://wiki.fysik.dtu.dk/ase/). - -I you want to interface your Python program with `xtb` and are afraid to add -a such heavy module as ASE as your dependency, you can interface directly -to the Python-API which lives in `xtb.interface` and does only require -`numpy` and `ctypes`. +```python +from xtb.interface import _libxtb, XTBLibrary, load_library +_libxtb = XTBLibrary(libarary=load_library('path/to/libxtb') +``` -Coming from every other C-compatible language I would recommend to wrap -the C-API in a similar way like I did for Python. +to point the interface to the correct library. From 0ca975b79f680590910538b7293c82e90a2f418e Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Fri, 22 May 2020 17:03:28 +0200 Subject: [PATCH 3/4] Version bump --- CMakeLists.txt | 2 +- meson.build | 2 +- python/setup.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 53e0f22f6..25f5a2823 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.9) # Setup the XTB Project project(xtb - VERSION 6.2.2 + VERSION 6.3.0 ) enable_language(Fortran) enable_testing() diff --git a/meson.build b/meson.build index 0f49a8f44..66a54dc1c 100644 --- a/meson.build +++ b/meson.build @@ -19,7 +19,7 @@ project( 'xtb', 'fortran', 'c', - version: '6.2.3', + version: '6.3.0', license: 'LGPL3', meson_version: '>=0.51', default_options: [ diff --git a/python/setup.py b/python/setup.py index 5d33dbabb..682936e4d 100644 --- a/python/setup.py +++ b/python/setup.py @@ -21,7 +21,7 @@ long_description = fh.read() setup(name="xtb", - version="6.2.3", + version="6.3.0", author="Sebastian Ehlert", author_email="ehlert@thch.uni-bonn.de", description="Wrapper for the extended tight binding program", From 14291f81fea08ba2b9340bba0feaab5cc35675c9 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Fri, 22 May 2020 17:09:20 +0200 Subject: [PATCH 4/4] Remove global import --- src/main/property.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/main/property.f90 b/src/main/property.f90 index 82305eb54..ba0389bf8 100644 --- a/src/main/property.f90 +++ b/src/main/property.f90 @@ -72,7 +72,6 @@ end subroutine write_energy_gff subroutine main_property & (iunit,mol,wfx,basis,xtbData,res,acc) - use xtb_mctc_global, only : persistentEnv use xtb_mctc_convert !! ========================================================================