From aa9f2e224ff84313515dff5e4d312b30acd5739e Mon Sep 17 00:00:00 2001 From: awvwgk <28669218+awvwgk@users.noreply.github.com> Date: Thu, 3 Oct 2019 12:55:23 +0200 Subject: [PATCH] start reworking Python API - add interface class wrapping the CDLL - sanity check for ndarrays (type, size) - rework ASE Calculator to be independent of the ctypes interface - add unit test for the three Hamiltonians - add mypy, pytest and pylint to meson.build - ignore __pycache__ in the working directory - moved old xtb.py to scripts - flake8 test in meson.build - additional test for the library interface - added GBSA Python API - GFN1 return now charges, dipole moments and bond orders - unified interface for molecular and periodic GFN0 - updated setup.py and corresponding README - allows passing number of unpaired electrons to the API - fix test with unpaired electrons - added stress tensor as property to GFN0 Calculator - updated description for Python package - added a sanity check for multiplicities and electron number - adjusted optimization script for new setup --- .gitignore | 6 + TESTSUITE/cpp_api_example.cpp | 4 +- TESTSUITE/gfn0.f90 | 9 +- TESTSUITE/gfn1.f90 | 31 +- TESTSUITE/gfn2.f90 | 16 +- TESTSUITE/peeq.f90 | 8 +- include/xtb.h | 140 +++++-- meson.build | 28 ++ python/LICENSE | 165 ++++++++ python/README.md | 91 +++++ python/setup.py | 43 ++ python/xtb.py | 589 ---------------------------- python/xtb/__init__.py | 23 ++ python/xtb/calculators.py | 260 ++++++++++++ python/xtb/interface.py | 433 ++++++++++++++++++++ python/xtb/solvation.py | 131 +++++++ python/xtb/test/test_calculators.py | 117 ++++++ python/xtb/test/test_interface.py | 174 ++++++++ python/xtb/test/test_solvation.py | 125 ++++++ scripts/xtb_opt.py | 280 +++++++++++++ xtb/c_api.f90 | 219 ++++++++--- xtb/calculator.f90 | 217 ++-------- xtb/scf_module.f90 | 38 +- xtb/tbdef_molecule.f90 | 3 + xtb/tbdef_options.f90 | 19 + 25 files changed, 2272 insertions(+), 897 deletions(-) create mode 100644 python/LICENSE create mode 100644 python/README.md create mode 100644 python/setup.py delete mode 100755 python/xtb.py create mode 100644 python/xtb/__init__.py create mode 100755 python/xtb/calculators.py create mode 100644 python/xtb/interface.py create mode 100644 python/xtb/solvation.py create mode 100644 python/xtb/test/test_calculators.py create mode 100644 python/xtb/test/test_interface.py create mode 100644 python/xtb/test/test_solvation.py create mode 100755 scripts/xtb_opt.py diff --git a/.gitignore b/.gitignore index d39dd76e1..6c3e9e51b 100644 --- a/.gitignore +++ b/.gitignore @@ -59,3 +59,9 @@ xtbopt.xyz xtbopt.sdf xtbopt.coord xtbopt.log + +# python build artefacts +__pycache__/ +python/dist +python/build +*.egg-info diff --git a/TESTSUITE/cpp_api_example.cpp b/TESTSUITE/cpp_api_example.cpp index 01ae625d3..0ddeca080 100644 --- a/TESTSUITE/cpp_api_example.cpp +++ b/TESTSUITE/cpp_api_example.cpp @@ -22,7 +22,7 @@ main (int argc, char **argv) 1.92825631079613, 0.00000000000000,-2.53624948351102, 0.00000000000000, 0.00000000000000, 5.23010455462158}; - const xtb::SCC_options opt {2, 0, 1.0, 300.0, true, false, 30, "none"}; + const xtb::SCC_options opt {2, 0, 1.0, 300.0, true, false, true, 30, "none"}; double energy {0.0}; double dipole[3] {0.0}; @@ -30,7 +30,7 @@ main (int argc, char **argv) double qp[6*natoms] {0.0}; double wbo[natoms*natoms] {0.0}; - int stat = xtb::GFN2_calculation(&natoms, attyp, &charge, coord, &opt, "-", + int stat = xtb::GFN2_calculation(&natoms, attyp, &charge, nullptr, coord, &opt, "-", &energy, nullptr, dipole, q, nullptr, qp, wbo); assert(stat == 0); diff --git a/TESTSUITE/gfn0.f90 b/TESTSUITE/gfn0.f90 index 1cf8f625f..9436680b2 100644 --- a/TESTSUITE/gfn0.f90 +++ b/TESTSUITE/gfn0.f90 @@ -64,7 +64,7 @@ subroutine test_gfn0_sp mol%xyz = xyz mol%chrg = 0.0_wp call mol%set_nuclear_charge - call mol%calculate_distances + call mol%update wfn%nel = idint(sum(mol%z)) wfn%nopen = 0 @@ -160,7 +160,7 @@ subroutine test_gfn0_api real(wp) :: energy real(wp) :: hl_gap - real(wp) :: gradlatt(3,3) + real(wp) :: dum(3,3) real(wp),allocatable :: gradient(:,:) ! setup the environment variables @@ -169,14 +169,15 @@ subroutine test_gfn0_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp gradient = 0.0_wp call gfn0_calculation & - (istdout,env,opt,mol,gfn,hl_gap,energy,gradient) + (istdout,env,opt,mol,gfn,hl_gap,energy,gradient,dum,dum) call assert_close(hl_gap, 5.5384029314207_wp,thr) call assert_close(energy,-8.6908532561691_wp,thr) diff --git a/TESTSUITE/gfn1.f90 b/TESTSUITE/gfn1.f90 index ae92d4278..6f722b8e7 100644 --- a/TESTSUITE/gfn1.f90 +++ b/TESTSUITE/gfn1.f90 @@ -54,12 +54,13 @@ subroutine test_gfn1_scc mol%xyz = xyz mol%chrg = 0.0_wp call mol%set_nuclear_charge + call mol%update wfn%nel = idint(sum(mol%z)) wfn%nopen = 0 allocate( g(3,mol%n), source = 0.0_wp ) - + call use_parameterset('.param_gfn.xtb',globpar,okpar) call assert(okpar) @@ -124,6 +125,7 @@ subroutine test_gfn1_api use tbdef_molecule use tbdef_param use tbdef_pcem + use tbdef_wavefunction use tb_calculators @@ -147,6 +149,7 @@ subroutine test_gfn1_api type(tb_environment) :: env type(gfn_parameter) :: gfn type(tb_pcem) :: pcem + type(tb_wavefunction):: wfn real(wp) :: energy real(wp) :: hl_gap @@ -158,14 +161,15 @@ subroutine test_gfn1_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp gradient = 0.0_wp call gfn1_calculation & - (istdout,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (istdout,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) call assert_close(hl_gap, 5.6067613075402_wp,thr) call assert_close(energy,-8.4156335932985_wp,thr) @@ -189,6 +193,7 @@ subroutine test_gfn1gbsa_api use tbdef_molecule use tbdef_param use tbdef_pcem + use tbdef_wavefunction use tb_calculators @@ -214,6 +219,7 @@ subroutine test_gfn1gbsa_api type(tb_environment) :: env type(gfn_parameter) :: gfn type(tb_pcem) :: pcem + type(tb_wavefunction):: wfn real(wp) :: energy real(wp) :: hl_gap @@ -225,14 +231,15 @@ subroutine test_gfn1gbsa_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp gradient = 0.0_wp call gfn1_calculation & - (istdout,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (istdout,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) call assert_close(hl_gap, 6.641641300724_wp,1e-4_wp) call assert_close(energy,-14.215790820910_wp,thr) @@ -255,6 +262,7 @@ subroutine test_gfn1_pcem_api use tbdef_molecule use tbdef_param use tbdef_pcem + use tbdef_wavefunction use aoparam @@ -288,6 +296,7 @@ subroutine test_gfn1_pcem_api type(tb_environment) :: env type(gfn_parameter) :: gfn type(tb_pcem) :: pcem + type(tb_wavefunction):: wfn real(wp) :: energy real(wp) :: hl_gap @@ -299,14 +308,15 @@ subroutine test_gfn1_pcem_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp gradient = 0.0_wp call gfn1_calculation & - (istdout,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (istdout,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) call assert_close(hl_gap, 9.0155275960407_wp,thr*10) call assert_close(energy,-23.113490916186_wp,thr) @@ -325,7 +335,8 @@ subroutine test_gfn1_pcem_api call mol%allocate(nat2) mol%at = at(:nat2) mol%xyz = xyz(:,:nat2) - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update call pcem%allocate(nat2) pcem%xyz = xyz(:,nat2+1:) @@ -335,7 +346,7 @@ subroutine test_gfn1_pcem_api pcem%grd = 0.0_wp call gfn1_calculation & - (istdout,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (istdout,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) call assert_close(hl_gap, 8.7253450666347_wp,thr) call assert_close(energy,-11.559896105984_wp,thr) @@ -359,7 +370,7 @@ subroutine test_gfn1_pcem_api pcem%gam = 999.0_wp ! point charges call gfn1_calculation & - (istdout,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (istdout,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) call assert_close(hl_gap, 8.9183046297437_wp,thr) call assert_close(energy,-11.565012263827_wp,thr) diff --git a/TESTSUITE/gfn2.f90 b/TESTSUITE/gfn2.f90 index 755b88f31..ff934bdc2 100644 --- a/TESTSUITE/gfn2.f90 +++ b/TESTSUITE/gfn2.f90 @@ -54,6 +54,7 @@ subroutine test_gfn2_scc mol%xyz = xyz mol%chrg = 0.0_wp call mol%set_nuclear_charge + call mol%update wfn%nel = idint(sum(mol%z)) wfn%nopen = 0 @@ -167,7 +168,8 @@ subroutine test_gfn2_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp @@ -239,7 +241,8 @@ subroutine test_gfn2gbsa_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp @@ -309,7 +312,8 @@ subroutine test_gfn2salt_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp @@ -390,7 +394,8 @@ subroutine test_gfn2_pcem_api call mol%allocate(nat) mol%at = at mol%xyz = xyz - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp @@ -416,7 +421,8 @@ subroutine test_gfn2_pcem_api call mol%allocate(nat2) mol%at = at(:nat2) mol%xyz = xyz(:,:nat2) - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update call pcem%allocate(nat2) pcem%xyz = xyz(:,nat2+1:) diff --git a/TESTSUITE/peeq.f90 b/TESTSUITE/peeq.f90 index 639dcb333..61bee8151 100644 --- a/TESTSUITE/peeq.f90 +++ b/TESTSUITE/peeq.f90 @@ -202,6 +202,7 @@ subroutine test_peeq_api real(wp) :: energy real(wp) :: hl_gap real(wp) :: gradlatt(3,3) + real(wp) :: stress(3,3) real(wp),allocatable :: gradient(:,:) ! setup the environment variables @@ -217,7 +218,8 @@ subroutine test_peeq_api call dlat_to_cell(lattice,mol%cellpar) call dlat_to_rlat(lattice,mol%rec_lat) call coord_trafo(nat,lattice,abc,mol%xyz) - call mol%calculate_distances + call mol%set_nuclear_charge + call mol%update allocate(gradient(3,mol%n)) energy = 0.0_wp @@ -226,8 +228,8 @@ subroutine test_peeq_api call mctc_mute - call peeq_calculation & - (istdout,env,opt,mol,gfn,hl_gap,energy,gradient,gradlatt) + call gfn0_calculation & + (istdout,env,opt,mol,gfn,hl_gap,energy,gradient,stress,gradlatt) call assert_close(hl_gap, 4.8620892163953_wp,thr) call assert_close(energy,-8.4930019025474_wp,thr) diff --git a/include/xtb.h b/include/xtb.h index c5346118d..77f851729 100644 --- a/include/xtb.h +++ b/include/xtb.h @@ -29,6 +29,7 @@ typedef struct _PEEQ_options { double etemp; bool grad; bool ccm; + char solvent[20]; } PEEQ_options; typedef struct _SCC_options { @@ -38,6 +39,7 @@ typedef struct _SCC_options { double etemp; bool grad; bool restart; + bool ccm; int maxiter; char solvent[20]; } SCC_options; @@ -47,50 +49,128 @@ extern "C" { #endif extern int -GFN0_PBC_calculation(const int* natoms, const int* attyp, const double* charge, - const double* coord, const double* lattice, const bool* pbc, - const PEEQ_options* opt, const char* output, - double* energy, double* grad, double* glat); +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 double* coord, const SCC_options* opt, const char* output, - double* energy, double* grad, double* dipole, double* q, - double* dipm, double* qp, double* wbo); +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 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); +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 double* coord, const SCC_options* opt, const char* output, - double* energy, double* grad); +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 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); +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 double* coord, const PEEQ_options* opt, const char* output, - double* energy, double* grad); +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); +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); +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); #ifdef __cplusplus } diff --git a/meson.build b/meson.build index e030fde5b..449d01579 100644 --- a/meson.build +++ b/meson.build @@ -466,3 +466,31 @@ test('GFN0-xTB: SP', xtb_test,args: ['gfn0','sp']) test('GFN0-xTB: API',xtb_test,args: ['gfn0','api']) test('GFN0-xTB: SP (PBC)', xtb_test,args: ['peeq','sp']) test('GFN0-xTB: API (PBC)', xtb_test,args: ['peeq','api']) + +py_xtb = meson.current_source_dir() / 'python' / 'xtb' +pypath = environment() +pypath.prepend('LD_LIBRARY_PATH', meson.current_build_dir()) +pypath.prepend('PYTHONPATH', meson.current_source_dir() / 'python') + +pytest = find_program('pytest', required: false) +if pytest.found() + test('pytest: xtb.py', pytest, args: ['--pyargs', 'xtb'], env: pypath) +endif + +mypy = find_program('mypy', required: false) +if mypy.found() + test('mypy: xtb.py', mypy, args: [py_xtb, '--ignore-missing-imports'], + env: pypath) +endif + +pylint = find_program('pylint', required: false) +if pylint.found() + test('pylint: xtb.py', pylint, args: [py_xtb, '--disable=invalid-names,duplicate-code'], + env: pypath) +endif + +flake8 = find_program('flake8', required: false) +if pylint.found() + test('flake8: xtb.py', flake8, args: [py_xtb, '--max-line-length=83'], + env: pypath) +endif diff --git a/python/LICENSE b/python/LICENSE new file mode 100644 index 000000000..0a041280b --- /dev/null +++ b/python/LICENSE @@ -0,0 +1,165 @@ + GNU LESSER GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + + This version of the GNU Lesser General Public License incorporates +the terms and conditions of version 3 of the GNU General Public +License, supplemented by the additional permissions listed below. + + 0. Additional Definitions. + + As used herein, "this License" refers to version 3 of the GNU Lesser +General Public License, and the "GNU GPL" refers to version 3 of the GNU +General Public License. + + "The Library" refers to a covered work governed by this License, +other than an Application or a Combined Work as defined below. + + An "Application" is any work that makes use of an interface provided +by the Library, but which is not otherwise based on the Library. +Defining a subclass of a class defined by the Library is deemed a mode +of using an interface provided by the Library. + + A "Combined Work" is a work produced by combining or linking an +Application with the Library. The particular version of the Library +with which the Combined Work was made is also called the "Linked +Version". + + The "Minimal Corresponding Source" for a Combined Work means the +Corresponding Source for the Combined Work, excluding any source code +for portions of the Combined Work that, considered in isolation, are +based on the Application, and not on the Linked Version. + + The "Corresponding Application Code" for a Combined Work means the +object code and/or source code for the Application, including any data +and utility programs needed for reproducing the Combined Work from the +Application, but excluding the System Libraries of the Combined Work. + + 1. Exception to Section 3 of the GNU GPL. + + You may convey a covered work under sections 3 and 4 of this License +without being bound by section 3 of the GNU GPL. + + 2. Conveying Modified Versions. + + If you modify a copy of the Library, and, in your modifications, a +facility refers to a function or data to be supplied by an Application +that uses the facility (other than as an argument passed when the +facility is invoked), then you may convey a copy of the modified +version: + + a) under this License, provided that you make a good faith effort to + ensure that, in the event an Application does not supply the + function or data, the facility still operates, and performs + whatever part of its purpose remains meaningful, or + + b) under the GNU GPL, with none of the additional permissions of + this License applicable to that copy. + + 3. Object Code Incorporating Material from Library Header Files. + + The object code form of an Application may incorporate material from +a header file that is part of the Library. You may convey such object +code under terms of your choice, provided that, if the incorporated +material is not limited to numerical parameters, data structure +layouts and accessors, or small macros, inline functions and templates +(ten or fewer lines in length), you do both of the following: + + a) Give prominent notice with each copy of the object code that the + Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the object code with a copy of the GNU GPL and this license + document. + + 4. Combined Works. + + You may convey a Combined Work under terms of your choice that, +taken together, effectively do not restrict modification of the +portions of the Library contained in the Combined Work and reverse +engineering for debugging such modifications, if you also do each of +the following: + + a) Give prominent notice with each copy of the Combined Work that + the Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the Combined Work with a copy of the GNU GPL and this license + document. + + c) For a Combined Work that displays copyright notices during + execution, include the copyright notice for the Library among + these notices, as well as a reference directing the user to the + copies of the GNU GPL and this license document. + + d) Do one of the following: + + 0) Convey the Minimal Corresponding Source under the terms of this + License, and the Corresponding Application Code in a form + suitable for, and under terms that permit, the user to + recombine or relink the Application with a modified version of + the Linked Version to produce a modified Combined Work, in the + manner specified by section 6 of the GNU GPL for conveying + Corresponding Source. + + 1) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (a) uses at run time + a copy of the Library already present on the user's computer + system, and (b) will operate properly with a modified version + of the Library that is interface-compatible with the Linked + Version. + + e) Provide Installation Information, but only if you would otherwise + be required to provide such information under section 6 of the + GNU GPL, and only to the extent that such information is + necessary to install and execute a modified version of the + Combined Work produced by recombining or relinking the + Application with a modified version of the Linked Version. (If + you use option 4d0, the Installation Information must accompany + the Minimal Corresponding Source and Corresponding Application + Code. If you use option 4d1, you must provide the Installation + Information in the manner specified by section 6 of the GNU GPL + for conveying Corresponding Source.) + + 5. Combined Libraries. + + You may place library facilities that are a work based on the +Library side by side in a single library together with other library +facilities that are not Applications and are not covered by this +License, and convey such a combined library under terms of your +choice, if you do both of the following: + + a) Accompany the combined library with a copy of the same work based + on the Library, uncombined with any other library facilities, + conveyed under the terms of this License. + + b) Give prominent notice with the combined library that part of it + is a work based on the Library, and explaining where to find the + accompanying uncombined form of the same work. + + 6. Revised Versions of the GNU Lesser General Public License. + + The Free Software Foundation may publish revised and/or new versions +of the GNU Lesser General Public License from time to time. Such new +versions will be similar in spirit to the present version, but may +differ in detail to address new problems or concerns. + + Each version is given a distinguishing version number. If the +Library as you received it specifies that a certain numbered version +of the GNU Lesser General Public License "or any later version" +applies to it, you have the option of following the terms and +conditions either of that published version or of any later version +published by the Free Software Foundation. If the Library as you +received it does not specify a version number of the GNU Lesser +General Public License, you may choose any version of the GNU Lesser +General Public License ever published by the Free Software Foundation. + + If the Library as you received it specifies that a proxy can decide +whether future versions of the GNU Lesser General Public License shall +apply, that proxy's public statement of acceptance of any version is +permanent authorization for you to choose that version for the +Library. diff --git a/python/README.md b/python/README.md new file mode 100644 index 000000000..23c1a3921 --- /dev/null +++ b/python/README.md @@ -0,0 +1,91 @@ +# Python Wrapper for the Extended Tight Binding Program + +This is the Python-side of the wrapper 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 accessable +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 usa we recommend to use `xtb` together with the +Atomic Simulation Model (ASE) like + +```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() +``` + +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 accessable from the atoms object and must be taken directely from the +calculator object. + +Additionally there is a solvation model calculator based on a +generalized Born model with a solvent accessable 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 + +```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() +``` + +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`. + +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`. + +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. + +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 directely +to the Python-API which lives in `xtb.interface` and does only require +`numpy` and `ctypes`. + +Coming from every other C-compatible language I would recommend to wrap +the C-API in a similar way like I did for Python. diff --git a/python/setup.py b/python/setup.py new file mode 100644 index 000000000..09249cef5 --- /dev/null +++ b/python/setup.py @@ -0,0 +1,43 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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 . + +from setuptools import setup, find_packages + +with open("README.md", "r") as fh: + long_description = fh.read() + +setup(name="xtb", + version="6.2.1", + author="Sebastian Ehlert", + author_email="ehlert@thch.uni-bonn.de", + description="Wrapper for the extended tight binding program", + long_description=long_description, + long_description_content_type="text/markdown", + keywords="tight-binding", + url="https://github.com/grimme-lab/xtb", + project_urls={ + "Documentation": "https://xtb-docs.readthedocs.io/en/latest/contents.html", + }, + packages=find_packages(), + install_requires=['ase'], + classifiers=[ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: LGPL3", + "Operating System :: Linux", + ], + python_requires=">=3.5", +) diff --git a/python/xtb.py b/python/xtb.py deleted file mode 100755 index 8e5153b28..000000000 --- a/python/xtb.py +++ /dev/null @@ -1,589 +0,0 @@ -#!/usr/bin/env python - -from __future__ import print_function - -import numpy as np -import ase -import ctypes - -from ctypes import cdll, Structure, c_double, c_int, c_bool, POINTER, c_char_p - -from ase import Atoms -from ase.calculators.calculator import Calculator, all_changes -from ase.units import Hartree, Bohr - -class SCC_options(Structure): - _fields_ = [('prlevel', c_int), - ('parallel',c_int), - ('acc', c_double), - ('etemp', c_double), - ('grad', c_bool), - ('restart', c_bool), - ('maxiter', c_int), - ('solvent', ctypes.c_char*20)] - -class PEEQ_options(Structure): - _fields_ = [('prlevel', c_int), - ('parallel',c_int), - ('acc', c_double), - ('etemp', c_double), - ('grad', c_bool), - ('ccm', c_bool)] - -class XTB(Calculator): - implemented_properties = ['energy', 'free_energy', 'forces'] - command = None - - def __init__(self, charge = None, accuracy = None, temperature = None, - print_level = None, - restart=None, ignore_bad_restart_file=False, - label="saw", atoms=None, command=None, debug=False, **kwargs): - """Construct the xtb-calculator object.""" - - self._debug = debug - self._lib = None - self.label = None - self.parameters = None - self.results = None - self.atoms = None - - if charge is not None: - self.charge = charge - else: - self.charge = 0.0 - - if accuracy is not None: - self.accuracy = accuracy - else: - self.accuracy = 1.0 - - if temperature is not None: - self.temperature = temperature - else: - self.temperature = 300.0 - - if print_level is not None: - self.print_level = print_level - else: - self.print_level = 2 - - if command is not None: - self.command = command - else: - self.command = 'xtb' # default - - Calculator.__init__(self, restart, ignore_bad_restart_file, - label, atoms, **kwargs) - - if restart is not None: - try: - self.read(restart) - except: - if ignore_bad_restart_file: - self.reset() - else: - raise - - def set(self, **kwargs): - """Set parameters like set(key1=value1, key2=value2, ...).""" - changed_parameters = Calculator.set(self, **kwargs) - if changed_parameters: - self.reset() - - def write(self, label): - """Write atoms, parameters and calculated results into restart files.""" - if self._debug: - print("Writting restart to: ", label) - self.atoms.write(label + '_restart.traj') - self.parameters.write(label + '_params.ase') - open(label + '_results.ase', 'w').write(repr(self.results)) - - def read(self, label): - """Read atoms, parameters and calculated results from restart files.""" - self.atoms = ase.io.read(label + '_restart.traj') - self.parameters = Parameters.read(label + '_params.ase') - results_txt = open(label + '_results.ase').read() - self.results = eval(results_txt, {'array': np.array}) - - def _open_lib(self): - # load xtb-library - self._lib = cdll.LoadLibrary('libxtb.so') - - c_int_p = POINTER(c_int) - c_bool_p = POINTER(c_bool) - c_double_p = POINTER(c_double) - - # define periodic GFN0-xTB interface - self._lib.GFN0_PBC_calculation.argtypes = [c_int_p, c_int_p, - c_double_p, c_double_p, c_double_p, - c_bool_p, POINTER(PEEQ_options), c_char_p, - c_double_p, c_double_p, c_double_p] - - # define GFN0-xTB interface - self._lib.GFN0_calculation.argtypes = [c_int_p, c_int_p, - c_double_p, c_double_p, POINTER(PEEQ_options), c_char_p, - c_double_p, c_double_p] - - # define GFN1-xTB interface - self._lib.GFN1_calculation.argtypes = [c_int_p, c_int_p, - c_double_p, c_double_p, POINTER(SCC_options), c_char_p, - c_double_p, c_double_p] - - # define GFN2-xTB interface - self._lib.GFN2_calculation.argtypes = [c_int_p, c_int_p, - c_double_p, c_double_p, POINTER(SCC_options), c_char_p, - c_double_p, c_double_p, c_double_p, c_double_p, c_double_p, - c_double_p, c_double_p] - - def calculate(self, atoms=None, properties=None, system_changes=all_changes): - """Do the calculation.""" - - if not properties: - properties = ['energy'] - Calculator.calculate(self, atoms, properties, system_changes) - -class GFN0(XTB): - implemented_properties = ['energy', 'free_energy', 'forces', 'stress'] - command = None - - def __init__(self, charge = None, accuracy = None, temperature = None, - print_level = None, - restart=None, ignore_bad_restart_file=False, - label="saw", atoms=None, command=None, debug=False, **kwargs): - """Construct the xtb-calculator object.""" - - XTB.__init__(self, charge, accuracy, temperature, print_level, - restart, ignore_bad_restart_file, label, atoms, **kwargs) - - - def calculate(self, atoms=None, properties=None, system_changes=all_changes): - """Do the calculation.""" - - if not properties: - properties = ['energy'] - XTB.calculate(self, atoms, properties, system_changes) - - if self._debug: - print("system_changes:", system_changes) - - if self._lib is None: - self._open_lib() - - natoms = c_int(atoms.get_number_of_atoms()) - attyp = np.array(atoms.get_atomic_numbers(), dtype=np.int32) - attyp_p = attyp.ctypes.data_as(POINTER(c_int)) - coord = atoms.get_positions()/Bohr # from Angstrom to Bohr - coord_p = coord.ctypes.data_as(POINTER(c_double)) - charge = c_double(self.charge) - - energy = c_double(0.0) - gradient = np.zeros((3,natoms.value),order='F') - grad_p = gradient.ctypes.data_as(POINTER(c_double)) - - opt = PEEQ_options(c_int(self.print_level), - c_int(0), # parallel -> automatic - c_double(self.accuracy), - c_double(self.temperature), - c_bool(True), # gradient -> always - c_bool(True)) # CCM -> yes - - outfile = "gfn0.out".encode('utf-8') - - stat = self._lib.GFN0_calculation(natoms,attyp_p,charge, - coord_p,opt,outfile,energy,grad_p) - # in case Fortran behaves we find a useful return value - if (stat != 0): - raise Exception("xtb terminated in error.") - - self.results['energy'] = energy.value * Hartree - self.results['free_energy'] = energy.value * Hartree - self.results['forces'] = - gradient.T * Hartree / Bohr - -class GFN1(XTB): - implemented_properties = ['energy', 'free_energy', 'forces'] - command = None - - def __init__(self, charge = None, accuracy = None, temperature = None, - max_iterations = None, print_level = None, solvent = None, - restart=None, ignore_bad_restart_file=False, - label="saw", atoms=None, command=None, debug=False, **kwargs): - """Construct the GFN2-calculator object.""" - - if max_iterations is not None: - self.max_iterations = max_iterations - else: - self.max_iterations = 250 - - if solvent is not None: - self.solvent = solvent - else: - self.solvent = "none" - - XTB.__init__(self, charge, accuracy, temperature, print_level, - restart, ignore_bad_restart_file, label, atoms, **kwargs) - - def calculate(self, atoms=None, properties=None, system_changes=all_changes): - """Do the calculation.""" - - if not properties: - properties = ['energy'] - XTB.calculate(self, atoms, properties, system_changes) - - if self._debug: - print("system_changes:", system_changes) - - if self._lib is None: - self._open_lib() - - natoms = c_int(atoms.get_number_of_atoms()) - attyp = np.array(atoms.get_atomic_numbers(), dtype=np.int32) - attyp_p = attyp.ctypes.data_as(POINTER(c_int)) - coord = atoms.get_positions()/Bohr # from Angstrom to Bohr - coord_p = coord.ctypes.data_as(POINTER(c_double)) - charge = c_double(self.charge) - - energy = c_double(0.0) - gradient = np.zeros((3,natoms.value),order='F') - grad_p = gradient.ctypes.data_as(POINTER(c_double)) - - opt = SCC_options(c_int(self.print_level), - c_int(0), # parallel -> automatic - c_double(self.accuracy), - c_double(self.temperature), - c_bool(True), # gradient -> always - c_bool(False), # restart -> no - c_int(self.max_iterations), - self.solvent.encode('utf-8')) - - outfile = "gfn1.out".encode('utf-8') - - stat = self._lib.GFN1_calculation(natoms,attyp_p,charge, - coord_p,opt,outfile,energy,grad_p) - # in case Fortran behaves we find a useful return value - if (stat != 0): - raise Exception("xtb terminated in error.") - - - self.results['energy'] = energy.value * Hartree - self.results['free_energy'] = energy.value * Hartree - self.results['forces'] = - gradient.T * Hartree / Bohr - -class GFN2(XTB): - implemented_properties = ['energy', 'free_energy', 'forces', 'dipole', - 'charges', 'dipoles', 'quadrupoles', 'wiberg'] - command = None - - def __init__(self, charge = None, accuracy = None, temperature = None, - max_iterations = None, print_level = None, solvent = None, - restart=None, ignore_bad_restart_file=False, - label="saw", atoms=None, command=None, debug=False, **kwargs): - """Construct the GFN2-calculator object.""" - - if max_iterations is not None: - self.max_iterations = max_iterations - else: - self.max_iterations = 250 - - if solvent is not None: - self.solvent = solvent - else: - self.solvent = "none" - - XTB.__init__(self, charge, accuracy, temperature, print_level, - restart, ignore_bad_restart_file, label, atoms, **kwargs) - - def calculate(self, atoms=None, properties=None, system_changes=all_changes): - """Do the calculation.""" - - if not properties: - properties = ['energy'] - XTB.calculate(self, atoms, properties, system_changes) - - if self._debug: - print("system_changes:", system_changes) - - if self._lib is None: - self._open_lib() - - natoms = c_int(atoms.get_number_of_atoms()) - attyp = np.array(atoms.get_atomic_numbers(), dtype=np.int32) - attyp_p = attyp.ctypes.data_as(POINTER(c_int)) - coord = atoms.get_positions()/Bohr # from Angstrom to Bohr - coord_p = coord.ctypes.data_as(POINTER(c_double)) - charge = c_double(self.charge) - self.results['charges'] = np.zeros(natoms.value) - self.results['dipoles'] = np.zeros((3,natoms.value),order='F') - self.results['quadrupoles'] = np.zeros((6,natoms.value),order='F') - self.results['dipole'] = np.zeros(3) - self.results['wiberg'] = np.zeros((natoms.value,natoms.value),order='F') - q_p = self.results['charges'].ctypes.data_as(POINTER(c_double)) - dipm_p = self.results['dipoles'].ctypes.data_as(POINTER(c_double)) - qp_p = self.results['quadrupoles'].ctypes.data_as(POINTER(c_double)) - dipole_p = self.results['dipole'].ctypes.data_as(POINTER(c_double)) - wbo_p = self.results['wiberg'].ctypes.data_as(POINTER(c_double)) - - energy = c_double(0.0) - gradient = np.zeros((3,natoms.value),order='F') - grad_p = gradient.ctypes.data_as(POINTER(c_double)) - - opt = SCC_options(c_int(self.print_level), - c_int(0), # parallel -> automatic - c_double(self.accuracy), - c_double(self.temperature), - c_bool(True), # gradient -> always - c_bool(True), # restart -> if possible - c_int(self.max_iterations), - self.solvent.encode('utf-8')) - - outfile = "gfn2.out".encode('utf-8') - - stat = self._lib.GFN2_calculation(natoms,attyp_p,charge, - coord_p,opt,outfile,energy,grad_p,dipole_p,q_p,dipm_p,qp_p,wbo_p) - # in case Fortran behaves we find a useful return value - if (stat != 0): - raise Exception("xtb terminated in error.") - - self.results['energy'] = energy.value * Hartree - #print("-> GFN2:", energy.value * Hartree) - self.results['free_energy'] = energy.value * Hartree - self.results['forces'] = - gradient.T * Hartree / Bohr - - # some unit conversion for the sake of consistency - self.results['dipole'] = self.results['dipole'] * Bohr - self.results['dipoles'] = self.results['dipoles'] * Bohr - self.results['quadrupoles'] = self.results['quadrupoles'] * Bohr*Bohr - -class GFN0_PBC(XTB): - implemented_properties = ['energy', 'free_energy', 'forces', 'stress'] - command = None - - def __init__(self, charge = None, accuracy = None, temperature = None, - print_level = None, - restart=None, ignore_bad_restart_file=False, - label="saw", atoms=None, command=None, debug=False, **kwargs): - """Construct the xtb-calculator object.""" - - XTB.__init__(self, charge, accuracy, temperature, print_level, - restart, ignore_bad_restart_file, label, atoms, **kwargs) - - - def calculate(self, atoms=None, properties=None, system_changes=all_changes): - """Do the calculation.""" - - if not properties: - properties = ['energy'] - XTB.calculate(self, atoms, properties, system_changes) - - if self._debug: - print("system_changes:", system_changes) - - if self._lib is None: - self._open_lib() - - natoms = c_int(atoms.get_number_of_atoms()) - attyp = np.array(atoms.get_atomic_numbers(), dtype=np.int32) - attyp_p = attyp.ctypes.data_as(POINTER(c_int)) - coord = atoms.get_positions()/Bohr # from Angstrom to Bohr - coord_p = coord.ctypes.data_as(POINTER(c_double)) - lattice = atoms.get_cell()/Bohr # from Angstrom to Bohr - dlat_p = lattice.ctypes.data_as(POINTER(c_double)) - charge = c_double(self.charge) - - energy = c_double(0.0) - gradient = np.zeros((3,natoms.value),order='F') - grad_p = gradient.ctypes.data_as(POINTER(c_double)) - gradlatt = np.zeros((3,3),order='F') - glat_p = gradlatt.ctypes.data_as(POINTER(c_double)) - - opt = PEEQ_options(c_int(self.print_level), - c_int(0), # parallel -> automatic - c_double(self.accuracy), - c_double(self.temperature), - c_bool(True), # gradient -> always - c_bool(True)) # CCM -> yes - - #pbc = np.array([True,True,True], dtype=np.bool) - pbc = atoms.get_pbc() - pbc_p = pbc.ctypes.data_as(POINTER(c_bool)) - - outfile = "gfn0.out".encode('utf-8') - - # try/except is not working on Fortran - stat = self._lib.GFN0_PBC_calculation(natoms,attyp_p,charge, - coord_p,dlat_p,pbc_p,opt,outfile,energy,grad_p,glat_p) - # in case Fortran behaves we find a useful return value - if (stat != 0): - raise Exception("xtb terminated in error.") - - self.results['energy'] = energy.value * Hartree - self.results['free_energy'] = energy.value * Hartree - self.results['forces'] = - gradient.T * Hartree / Bohr - # xtb returns the cell gradient, so we have to backtransform - # taking into account that the stress tensor should be symmetric, - # we interface F and C ordered arrays and have to transpose we - # can simply drop all reshape and transpose steps by writing - stress = np.dot(self.atoms.cell, - gradlatt * Hartree / Bohr / self.atoms.get_volume()) - self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]] - -if __name__ == "__main__": - - import argparse - - def get_cmd_args(): - parser = argparse.ArgumentParser(prog='xtb', - description='Wrapper for xTB calculation.') - - parser.add_argument("-f","--format",dest="filetype",action="store", - default=None,help="Format of input geometry (automatic)") - - parser.add_argument(dest="filename",action="store", - help="Input geometry") - - parser.add_argument("-x","--method",dest="method",action="store", - default='gfn2',help="SQM method for calculation (gfn2)") - - parser.add_argument("-c","--charge",dest="charge",action="store", - type=int,default=0,help="Total charge (0)") - - parser.add_argument("--etemp",dest="etemp",action="store", - type=float,default=300.0,help="Electronic temperature (300K)") - - parser.add_argument("--accuracy",dest="acc",action="store", - type=float,default=1.0,help="Calculation accuracy (1.0)") - - parser.add_argument("--maxiter",dest="maxiter",action="store", - type=int,default=250,help="Maximum number of SCC iterations (250)") - - parser.add_argument("--lcovb",dest="rnn",action="store", - type=float,default=5.0,help="Longest covalent bond for preconditioning (5.0)") - - parser.add_argument("--econv",dest="econv",action="store", - type=float,default=5.4423e-4,help="Convergence threshold for energy") - - parser.add_argument("--gconv",dest="fconv",action="store", - type=float,default=3.88938e-05,help="Convergence threshold for forces") - - parser.add_argument("--linesearch",dest="armijo",action="store_true", - help="Use line search") - - parser.add_argument("--precon",dest="precon",action="store_true", - help="Use preconditioning") - - parser.add_argument("--optcell",dest="optcell",action="store_true", - help="Relax cellparameters") - - parser.add_argument("--logfile",dest="logfile",action="store", - default='-',help="File for logging information (STDOUT)") - - parser.add_argument("--trajectory",dest="trajectory",action="store", - default='xtbopt.traj',help="File for trajectory output (xtbopt.traj)") - - parser.add_argument("--solvent",dest="solvent",action="store", - default="none",help="Solvent for GBSA (none)") - - return parser.parse_args() - - import numpy - import ase - from ase.io import read, write - from ase.units import * - from ase.optimize.precon import Exp, PreconFIRE - from ase.constraints import ExpCellFilter - - # overwrite convergence thresholds of PreconFIRE optimizer - class PatchedOptimizer(PreconFIRE): - - def initialize(self): - PreconFIRE.initialize(self) - self.elast = None - - def run(self, steps=100000000, econv = 0.00054423,fconv = 3.88938e-05): - self.econv = econv - self.fconv = fconv - PreconFIRE.run(self, 0.05, steps) - - def converged(self, forces=None): - - # get current energy and check for convergence in energy change - ecurr = self.atoms.get_potential_energy() - if self.elast is not None: - ediff = ecurr - self.elast - econverged = ediff < 0.0 and abs(ediff) < self.econv - else: - econverged = False - # save potential energy - self.elast = ecurr - - # check for convergence in forces - if forces is None: - forces = self.atoms.get_forces() - fnorm = sqrt((forces**2).sum()) - fconverged = fnorm < self.fconv - - #print("--> norm(F):", fnorm, "converged?", fconverged) - #print("--> energy :", ecurr, "converged?", econverged) - - return econverged and fconverged - - args = get_cmd_args() - - if args.filetype is not None: - mol = ase.io.read(args.filename, format = str(args.filetype)) - else: - mol = ase.io.read(args.filename) - - if args.method == 'gfn0': - if mol.pbc.any(): - calc = GFN0_PBC(charge=args.charge,accuracy=args.acc, - temperature=args.etemp,print_level=2) - else: - calc = GFN0(charge=args.charge,accuracy=args.acc,temperature=args.etemp, - print_level=2) - elif args.method == 'gfn1': - if mol.pbc.any(): - raise Exception("GFN1-xTB is not available with PBC") - else: - calc = GFN1(charge=args.charge,accuracy=args.acc,temperature=args.etemp, - max_iterations=args.maxiter,print_level=2,solvent=args.solvent) - elif args.method == 'gfn2': - if mol.pbc.any(): - raise Exception("GFN2-xTB is not available with PBC") - else: - calc = GFN2(charge=args.charge,accuracy=args.acc,temperature=args.etemp, - max_iterations=args.maxiter,print_level=2,solvent=args.solvent) - else: - raise Exception("Method not implemented.") - mol.set_calculator(calc) - - e = mol.get_potential_energy() - - print('Initial energy: eV, Eh',e, e/Hartree) - - if args.precon: - precon = Exp(A=3, r_NN=args.rnn + 0.1, r_cut=args.rnn + 0.6) - else: - precon = None - - if args.optcell: - sf = ExpCellFilter(mol) - relax = PatchedOptimizer(sf, precon=precon, - trajectory=args.trajectory, logfile=args.logfile, use_armijo=args.armijo) - else: - relax = PatchedOptimizer(mol, precon=precon, - trajectory=args.trajectory, logfile=args.logfile, use_armijo=args.armijo) - - try: - relax.run(econv=args.econv,fconv=args.fconv) - except KeyboardInterrupt: - print('User got impatient') - except RuntimeError: - print('Optimization terminated due to internal error') - - e = mol.get_potential_energy() - print('Final energy: eV, Eh',e, e/Hartree) - - if (mol.pbc.any()): - ase.io.write("xtbopt.POSCAR",mol,format='vasp') - else: - ase.io.write("xtbopt.xyz",mol,format='xyz') diff --git a/python/xtb/__init__.py b/python/xtb/__init__.py new file mode 100644 index 000000000..9c4e0c12a --- /dev/null +++ b/python/xtb/__init__.py @@ -0,0 +1,23 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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 . + +"""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"] diff --git a/python/xtb/calculators.py b/python/xtb/calculators.py new file mode 100755 index 000000000..2be47fc5e --- /dev/null +++ b/python/xtb/calculators.py @@ -0,0 +1,260 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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 + 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.""" + 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.""" + 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 new file mode 100644 index 000000000..1f4e1e593 --- /dev/null +++ b/python/xtb/interface.py @@ -0,0 +1,433 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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 . +"""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.util import find_library + +import numpy as np + +__all__ = ['SCCOptions', 'PEEQOptions', 'XTBLibrary'] +__XTB_LIBRARY_NAME__ = find_library('xtb') or 'libxtb.so' + + +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 np.ctypeslib.as_ctypes_type(array.dtype) != ctype: + raise ValueError("{} must be of {} compatible type" + .format(name, ctype.__class__.__name__)) + + +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] + + def __init__(self, library: Optional[CDLL] = None): + """construct library from CDLL object.""" + if library is None: + library = cdll.LoadLibrary(__XTB_LIBRARY_NAME__) + + self.library = library + 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") + if periodic: + 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) + 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, + } + + # 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, + } + + # 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)), + ) + stat = self.library.GBSA_calculation(*args) + + if stat != 0: + raise RuntimeError("GBSA calculation failed in xtb.") + return { + 'born': born, + 'sasa': sasa, + } + + 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)), + ) + + stat = self.library.GBSA_model_preload(*args) + + if stat != 0: + raise RuntimeError("GBSA parameters could not be loaded by xtb.") diff --git a/python/xtb/solvation.py b/python/xtb/solvation.py new file mode 100644 index 000000000..daa5c3606 --- /dev/null +++ b/python/xtb/solvation.py @@ -0,0 +1,131 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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.""" + 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 new file mode 100644 index 000000000..f61c5cc6b --- /dev/null +++ b/python/xtb/test/test_calculators.py @@ -0,0 +1,117 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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.09848149585026728), + "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 new file mode 100644 index 000000000..2e3d03a2f --- /dev/null +++ b/python/xtb/test/test_interface.py @@ -0,0 +1,174 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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 ctypes interface to xtb.""" + + +def test_library(): + """check if we can find the library and it looks okay""" + from ctypes import cdll + from ctypes.util import find_library + + name = find_library("xtb") + assert name is not None + # check if we can load this one + lib = cdll.LoadLibrary(name) + # 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 + + 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") + + +class Testing: + """test all defined interfaces""" + from ctypes import c_int, c_double + from numpy import array + from xtb.interface import XTBLibrary + + natoms = 24 + numbers = 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], + dtype=c_int, + ) + positions = 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], + ], + 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 + assert approx(gnorm, thr) == 0.00510542946913145 + + def test_gfn1_interface(self): + """check if the GFN1-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, + "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 + assert approx(gnorm, thr) == 0.004842110219908473 + + def test_gfn2_gbsa_interface(self): + """check if the GFN2-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, + "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, + ) + assert results + gnorm = abs(results["gradient"].flatten()).mean() + assert approx(results["energy"], thr) == -42.170584528028506 + assert approx(gnorm, thr) == 0.004685246019376688 + assert approx(results["charges"].sum(), thr) == 0.0 diff --git a/python/xtb/test/test_solvation.py b/python/xtb/test/test_solvation.py new file mode 100644 index 000000000..de096b50e --- /dev/null +++ b/python/xtb/test/test_solvation.py @@ -0,0 +1,125 @@ +# This file is part of xtb. +# +# Copyright (C) 2019 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/scripts/xtb_opt.py b/scripts/xtb_opt.py new file mode 100755 index 000000000..9689276fd --- /dev/null +++ b/scripts/xtb_opt.py @@ -0,0 +1,280 @@ +#!/usr/bin/env python +# This file is part of xtb. +# +# Copyright (C) 2019 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 . + +"""Wrapper script to perform geometry optimizations with xtb Python wrapper.""" + +if __name__ == "__main__": + import argparse + + def get_cmd_args(): + """Commandline arguments for xtb wrapper in script mode""" + parser = argparse.ArgumentParser( + prog="xtb", description="Wrapper for xTB calculation." + ) + + parser.add_argument( + "-f", + "--format", + dest="filetype", + action="store", + default=None, + help="Format of input geometry (automatic)", + ) + + parser.add_argument(dest="filename", action="store", help="Input geometry") + + parser.add_argument( + "-x", + "--method", + dest="method", + action="store", + default="gfn2", + help="SQM method for calculation (gfn2)", + ) + + parser.add_argument( + "-c", + "--charge", + dest="charge", + action="store", + type=int, + default=0, + help="Total charge (0)", + ) + + parser.add_argument( + "--etemp", + dest="etemp", + action="store", + type=float, + default=300.0, + help="Electronic temperature (300K)", + ) + + parser.add_argument( + "--accuracy", + dest="acc", + action="store", + type=float, + default=1.0, + help="Calculation accuracy (1.0)", + ) + + parser.add_argument( + "--maxiter", + dest="maxiter", + action="store", + type=int, + default=250, + help="Maximum number of SCC iterations (250)", + ) + + parser.add_argument( + "--lcovb", + dest="rnn", + action="store", + type=float, + default=5.0, + help="Longest covalent bond for preconditioning (5.0)", + ) + + parser.add_argument( + "--econv", + dest="econv", + action="store", + type=float, + default=5.4423e-4, + help="Convergence threshold for energy", + ) + + parser.add_argument( + "--gconv", + dest="fconv", + action="store", + type=float, + default=3.88938e-05, + help="Convergence threshold for forces", + ) + + parser.add_argument( + "--linesearch", + dest="armijo", + action="store_true", + help="Use line search" + ) + + parser.add_argument( + "--precon", + dest="precon", + action="store_true", + help="Use preconditioning" + ) + + parser.add_argument( + "--optcell", + dest="optcell", + action="store_true", + help="Relax cellparameters", + ) + + parser.add_argument( + "--logfile", + dest="logfile", + action="store", + default="-", + help="File for logging information (STDOUT)", + ) + + parser.add_argument( + "--trajectory", + dest="trajectory", + action="store", + default="xtbopt.traj", + help="File for trajectory output (xtbopt.traj)", + ) + + parser.add_argument( + "--solvent", + dest="solvent", + action="store", + default="none", + help="Solvent for GBSA (none)", + ) + + return parser.parse_args() + + from math import sqrt + from ase.io import read, write + from ase.optimize.precon import Exp, PreconFIRE + from ase.constraints import ExpCellFilter + from ase.units import Hartree + from xtb import GFN1, GFN2, GFN0 + + # overwrite convergence thresholds of PreconFIRE optimizer + class PatchedOptimizer(PreconFIRE): + """Patches the ASE optimizer to use a different convergence threshold""" + + econv = None + fconv = None + elast = None + + def initialize(self): + PreconFIRE.initialize(self) + self.elast = None + + def run( + self, steps=100000000, econv=0.00054423, fconv=3.88938e-05 + ): # pylint: disable=arguments-differ + self.econv = econv + self.fconv = fconv + PreconFIRE.run(self, 0.05, steps) + + def converged(self, forces=None): + + # get current energy and check for convergence in energy change + ecurr = self.atoms.get_potential_energy() + if self.elast is not None: + ediff = ecurr - self.elast + econverged = ediff < 0.0 and abs(ediff) < self.econv + else: + econverged = False + # save potential energy + self.elast = ecurr + + # check for convergence in forces + if forces is None: + forces = self.atoms.get_forces() + fnorm = sqrt((forces ** 2).sum()) + fconverged = fnorm < self.fconv + + # print("--> norm(F):", fnorm, "converged?", fconverged) + # print("--> energy :", ecurr, "converged?", econverged) + + return econverged and fconverged + + # pylint: disable=invalid-name + args = get_cmd_args() + + if args.filetype is not None: + mol = read(args.filename, format=str(args.filetype)) + else: + mol = read(args.filename) + + parameters = { + "accuracy": args.acc, + "electronic_temperature": args.etemp, + "max_iterations": args.maxiter, + "solvent": args.solvent, + "print_level": 2, + } + + if args.method == "gfn0": + calc = GFN0(**parameters) + elif args.method == "gfn1": + if mol.pbc.any(): + raise Exception("GFN1-xTB is not available with PBC") + calc = GFN1(**parameters) + elif args.method == "gfn2": + if mol.pbc.any(): + raise Exception("GFN2-xTB is not available with PBC") + calc = GFN2(**parameters) + else: + raise Exception("Method not implemented.") + mol.set_calculator(calc) + + e = mol.get_potential_energy() + + print("Initial energy: eV, Eh", e, e / Hartree) + + if args.precon: + precon = Exp(A=3, r_NN=args.rnn + 0.1, r_cut=args.rnn + 0.6) + else: + precon = None + + if args.optcell: + sf = ExpCellFilter(mol) + relax = PatchedOptimizer( + sf, + precon=precon, + trajectory=args.trajectory, + logfile=args.logfile, + use_armijo=args.armijo, + ) + else: + relax = PatchedOptimizer( + mol, + precon=precon, + trajectory=args.trajectory, + logfile=args.logfile, + use_armijo=args.armijo, + ) + + try: + relax.run(econv=args.econv, fconv=args.fconv) + except KeyboardInterrupt: + print("User got impatient") + except RuntimeError: + print("Optimization terminated due to internal error") + + e = mol.get_potential_energy() + print("Final energy: eV, Eh", e, e / Hartree) + + if mol.pbc.any(): + write("xtbopt.POSCAR", mol, format="vasp") + else: + write("xtbopt.xyz", mol, format="xyz") diff --git a/xtb/c_api.f90 b/xtb/c_api.f90 index f5e37e857..c67fbb47e 100644 --- a/xtb/c_api.f90 +++ b/xtb/c_api.f90 @@ -21,7 +21,7 @@ module c_api implicit none - !> some overloading for convience + !> some overloading for convenience interface c_return module procedure :: c_return_double_0d module procedure :: c_return_double_1d @@ -35,15 +35,22 @@ module c_api 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 contains function peeq_api & - & (natoms,attyp,charge,coord,lattice,pbc,opt_in,file_in,etot,grad,glat) & + & (natoms,attyp,charge,uhf,coord,lattice,pbc,opt_in,file_in, & + & etot,grad,stress,glat) & & result(status) bind(C,name="GFN0_PBC_calculation") use tbdef_molecule @@ -57,6 +64,7 @@ function peeq_api & 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) @@ -67,6 +75,7 @@ function peeq_api & 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(tb_molecule) :: mol @@ -82,6 +91,7 @@ function peeq_api & real(wp) :: energy real(wp) :: hl_gap real(wp) :: lattice_gradient(3,3) + real(wp) :: stress_tensor(3,3) real(wp),allocatable :: gradient(:,:) ! ==================================================================== @@ -94,7 +104,7 @@ function peeq_api & ! STEP 2: convert the options from C struct to actual Fortran type ! ==================================================================== opt = opt_in - + call c_string_convert(outfile, file_in) if (outfile.ne.'-'.and.opt%prlevel > 0) then @@ -126,9 +136,12 @@ function peeq_api & ! get atomtypes, coordinates and total charge mol%at = attyp mol%xyz = coord - mol%chrg = charge + call c_get(mol%chrg, charge, 0.0_wp) + call c_get(mol%uhf, uhf, 0) ! update everything from xyz and lattice + call mol%set_nuclear_charge + call mol%set_atomic_masses call mol%update ! ==================================================================== @@ -144,16 +157,13 @@ function peeq_api & ! shut down fatal errors from the MCTC library, so it will not kill the caller call mctc_mute - call peeq_calculation & - (iunit,env,opt,mol,gfn,hl_gap,energy,gradient,lattice_gradient) + call gfn0_calculation & + (iunit,env,opt,mol,gfn,hl_gap,energy,gradient,stress_tensor,lattice_gradient) ! check if the MCTC environment is still sane, if not tell the caller call mctc_sanity(sane) - if (.not.sane) then !! it's stark raving mad and on fire !! - - ! at least try to destroy the molecule structure data + if (.not.sane) then call finalize - status = 1 return endif @@ -164,6 +174,7 @@ function peeq_api & call c_return(etot, energy) call c_return(grad, gradient) call c_return(glat, lattice_gradient) + call c_return(stress, stress_tensor) call finalize @@ -178,7 +189,7 @@ end subroutine finalize end function peeq_api function gfn2_api & - & (natoms,attyp,charge,coord,opt_in,file_in, & + & (natoms,attyp,charge,uhf,coord,opt_in,file_in, & & etot,grad,dipole,q,dipm,qp,wbo) & & result(status) bind(C,name="GFN2_calculation") @@ -195,6 +206,7 @@ function gfn2_api & 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(*) @@ -235,7 +247,7 @@ function gfn2_api & ! STEP 2: convert the options from C struct to actual Fortran type ! ==================================================================== opt = opt_in - + call c_string_convert(outfile, file_in) if (outfile.ne.'-'.and.opt%prlevel > 0) then @@ -263,8 +275,13 @@ function gfn2_api & ! get atomtypes, coordinates and total charge mol%at = attyp mol%xyz = coord - mol%chrg = charge - call mol%calculate_distances + call c_get(mol%chrg, charge, 0.0_wp) + call c_get(mol%uhf, uhf, 0) + + ! update everything from xyz and lattice + call mol%set_nuclear_charge + call mol%set_atomic_masses + call mol%update ! ==================================================================== ! STEP 4: reserve some memory @@ -284,11 +301,8 @@ function gfn2_api & ! check if the MCTC environment is still sane, if not tell the caller call mctc_sanity(sane) - if (.not.sane) then !! it's stark raving mad and on fire !! - - ! at least try to destroy the molecule structure data + if (.not.sane) then call finalize - status = 1 return endif @@ -318,13 +332,14 @@ end subroutine finalize end function gfn2_api function gfn1_api & - & (natoms,attyp,charge,coord,opt_in,file_in,etot,grad) & + & (natoms,attyp,charge,uhf,coord,opt_in,file_in,etot,grad,dipole,q,wbo) & & result(status) bind(C,name="GFN1_calculation") use tbdef_molecule use tbdef_param use tbdef_options use tbdef_pcem + use tbdef_wavefunction use tb_calculators @@ -333,6 +348,7 @@ function gfn1_api & 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(*) @@ -341,12 +357,16 @@ function gfn1_api & 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) type(tb_molecule) :: mol type(gfn_parameter) :: gfn type(scc_options) :: opt type(tb_environment) :: env type(tb_pcem) :: pcem + type(tb_wavefunction):: wfn character(len=:),allocatable :: outfile @@ -395,8 +415,13 @@ function gfn1_api & ! get atomtypes, coordinates and total charge mol%at = attyp mol%xyz = coord - mol%chrg = charge - call mol%calculate_distances + call c_get(mol%chrg, charge, 0.0_wp) + call c_get(mol%uhf, uhf, 0) + + ! update everything from xyz and lattice + call mol%set_nuclear_charge + call mol%set_atomic_masses + call mol%update ! ==================================================================== ! STEP 4: reserve some memory @@ -412,15 +437,12 @@ function gfn1_api & call mctc_mute call gfn1_calculation & - (iunit,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (iunit,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) ! check if the MCTC environment is still sane, if not tell the caller call mctc_sanity(sane) - if (.not.sane) then !! it's stark raving mad and on fire !! - - ! at least try to destroy the molecule structure data + if (.not.sane) then call finalize - status = 1 return endif @@ -430,6 +452,9 @@ function gfn1_api & ! ==================================================================== 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)) call finalize @@ -438,13 +463,14 @@ function gfn1_api & contains subroutine finalize call mol%deallocate + call wfn%deallocate deallocate(gradient) if (iunit.ne.istdout) call close_file(iunit) end subroutine finalize end function gfn1_api function gfn0_api & - & (natoms,attyp,charge,coord,opt_in,file_in,etot,grad) & + & (natoms,attyp,charge,uhf,coord,opt_in,file_in,etot,grad) & & result(status) bind(C,name="GFN0_calculation") use tbdef_molecule @@ -458,6 +484,7 @@ function gfn0_api & 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(*) @@ -479,6 +506,7 @@ function gfn0_api & integer :: i real(wp) :: energy real(wp) :: hl_gap + real(wp) :: dum(3,3) real(wp),allocatable :: gradient(:,:) ! ==================================================================== @@ -491,7 +519,7 @@ function gfn0_api & ! STEP 2: convert the options from C struct to actual Fortran type ! ==================================================================== opt = opt_in - + call c_string_convert(outfile, file_in) if (outfile.ne.'-'.and.opt%prlevel > 0) then @@ -519,8 +547,13 @@ function gfn0_api & ! get atomtypes, coordinates and total charge mol%at = attyp mol%xyz = coord - mol%chrg = charge - call mol%calculate_distances + call c_get(mol%chrg, charge, 0.0_wp) + call c_get(mol%uhf, uhf, 0) + + ! update everything from xyz and lattice + call mol%set_nuclear_charge + call mol%set_atomic_masses + call mol%update ! ==================================================================== ! STEP 4: reserve some memory @@ -536,15 +569,12 @@ function gfn0_api & call mctc_mute call gfn0_calculation & - (iunit,env,opt,mol,gfn,hl_gap,energy,gradient) + (iunit,env,opt,mol,gfn,hl_gap,energy,gradient,dum,dum) ! check if the MCTC environment is still sane, if not tell the caller call mctc_sanity(sane) - if (.not.sane) then !! it's stark raving mad and on fire !! - - ! at least try to destroy the molecule structure data + if (.not.sane) then call finalize - status = 1 return endif @@ -568,7 +598,7 @@ end subroutine finalize end function gfn0_api function gfn2_pcem_api & - & (natoms,attyp,charge,coord,opt_in,file_in, & + & (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") @@ -587,6 +617,7 @@ function gfn2_pcem_api & 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(*) @@ -657,8 +688,13 @@ function gfn2_pcem_api & ! get atomtypes, coordinates and total charge mol%at = attyp mol%xyz = coord - mol%chrg = charge - call mol%calculate_distances + call c_get(mol%chrg, charge, 0.0_wp) + call c_get(mol%uhf, uhf, 0) + + ! update everything from xyz and lattice + call mol%set_nuclear_charge + call mol%set_atomic_masses + call mol%update call pcem%allocate(npc) pcem%q = pc_q @@ -687,11 +723,8 @@ function gfn2_pcem_api & ! check if the MCTC environment is still sane, if not tell the caller call mctc_sanity(sane) - if (.not.sane) then !! it's stark raving mad and on fire !! - - ! at least try to destroy the molecule structure data + if (.not.sane) then call finalize - status = 1 return endif @@ -711,13 +744,14 @@ function gfn2_pcem_api & subroutine finalize call mol%deallocate call pcem%deallocate + call wfn%deallocate deallocate(gradient) if (iunit.ne.istdout) call close_file(iunit) end subroutine finalize end function gfn2_pcem_api function gfn1_pcem_api & - & (natoms,attyp,charge,coord,opt_in,file_in, & + & (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") @@ -725,6 +759,7 @@ function gfn1_pcem_api & use tbdef_param use tbdef_options use tbdef_pcem + use tbdef_wavefunction use aoparam @@ -735,6 +770,7 @@ function gfn1_pcem_api & 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(*) @@ -756,6 +792,7 @@ function gfn1_pcem_api & type(scc_options) :: opt type(tb_environment) :: env type(tb_pcem) :: pcem + type(tb_wavefunction):: wfn character(len=:),allocatable :: outfile @@ -804,8 +841,13 @@ function gfn1_pcem_api & ! get atomtypes, coordinates and total charge mol%at = attyp mol%xyz = coord - mol%chrg = charge - call mol%calculate_distances + call c_get(mol%chrg, charge, 0.0_wp) + call c_get(mol%uhf, uhf, 0) + + ! update everything from xyz and lattice + call mol%set_nuclear_charge + call mol%set_atomic_masses + call mol%update call pcem%allocate(npc) pcem%q = pc_q @@ -830,15 +872,12 @@ function gfn1_pcem_api & call mctc_mute call gfn1_calculation & - (iunit,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (iunit,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) ! check if the MCTC environment is still sane, if not tell the caller call mctc_sanity(sane) - if (.not.sane) then !! it's stark raving mad and on fire !! - - ! at least try to destroy the molecule structure data + if (.not.sane) then call finalize - status = 1 return endif @@ -858,6 +897,7 @@ function gfn1_pcem_api & subroutine finalize call mol%deallocate call pcem%deallocate + call wfn%deallocate deallocate(gradient) if (iunit.ne.istdout) call close_file(iunit) end subroutine finalize @@ -874,9 +914,10 @@ function gbsa_model_preload_api & !> Dielectric data real(c_double), intent(in) :: epsv_in real(wp), allocatable :: epsv - !> Solvent density (g/cm^3) and molar mass (g/mol) + !> 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 @@ -1102,28 +1143,94 @@ subroutine c_string_convert(f_string, c_string) enddo end subroutine c_string_convert -subroutine c_get_double_0d(f_array, c_array) +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 +end subroutine c_get_double_0d_alloc -subroutine c_get_double_1d(f_array, c_array) +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 +end subroutine c_get_double_1d_alloc -subroutine c_get_double_2d(f_array, c_array) +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 + end module c_api diff --git a/xtb/calculator.f90 b/xtb/calculator.f90 index 465db5b44..8e5889165 100644 --- a/xtb/calculator.f90 +++ b/xtb/calculator.f90 @@ -20,8 +20,8 @@ module tb_calculators ! ======================================================================== !> periodic GFN0-xTB (PEEQ) calculation -subroutine peeq_calculation & - (iunit,env,opt,mol,gfn,hl_gap,energy,gradient,lattice_gradient) +subroutine gfn0_calculation & + (iunit,env,opt,mol,gfn,hl_gap,energy,gradient,stress,lattice_gradient) use iso_fortran_env, wp => real64 use mctc_systools @@ -33,13 +33,14 @@ subroutine peeq_calculation & use tbdef_param use tbdef_data - use setparam, only : gfn_method + use setparam, only : gfn_method, ngrida use aoparam, only : use_parameterset use ehtparam, only : import_basisset use pbc_tools use xbasis use peeq_module + use gbobc implicit none @@ -53,6 +54,7 @@ subroutine peeq_calculation & 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) @@ -83,12 +85,12 @@ subroutine peeq_calculation & ! we assume that the user provides a resonable molecule input ! -> all atoms are inside the unit cell, all data is set and consistent - ! get nuclear charges (without core) - call mol%set_nuclear_charge - call mol%set_atomic_masses - - wfn%nel = idint(sum(mol%z)) - mol%chrg - wfn%nopen = 0 + 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 @@ -132,6 +134,12 @@ subroutine peeq_calculation & call gfn0_prparam(iunit,mol%n,mol%at,param) 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) + endif + ! ==================================================================== ! STEP 3: expand our Slater basis set in contracted Gaussians ! ==================================================================== @@ -157,22 +165,23 @@ subroutine peeq_calculation & call peeq(iunit,mol,wfn,basis,param,hl_gap,opt%etemp,opt%prlevel,opt%grad, & & opt%ccm,opt%acc,energy,gradient,sigma,res) - ! some post processing for the sigma tensor, since there is a lot of - ! confusion in the literature regarding cell gradients, we have to decide - ! for one, here I pick the cell gradients over the stress tensor - inv_lat = mat_inv_3x3(mol%lattice) - call sigma_to_latgrad(sigma,inv_lat,lattice_gradient) + 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 peeq_calculation +end subroutine gfn0_calculation ! ======================================================================== !> GFN2-xTB calculation @@ -245,12 +254,12 @@ subroutine gfn2_calculation & ! we assume that the user provides a resonable molecule input ! -> all atoms are inside the unit cell, all data is set and consistent - ! get nuclear charges (without core) - call mol%set_nuclear_charge - call mol%set_atomic_masses - - wfn%nel = idint(sum(mol%z)) - mol%chrg - wfn%nopen = 0 + 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 @@ -288,8 +297,8 @@ subroutine gfn2_calculation & call gfn2_prparam(iunit,mol%n,mol%at,param) endif - if (len_trim(opt%solvent).gt.0 .and. opt%solvent.ne."none") then - lgbsa = .true. + 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) endif @@ -344,7 +353,7 @@ end subroutine gfn2_calculation ! ======================================================================== !> GFN1-xTB calculation subroutine gfn1_calculation & - (iunit,env,opt,mol,gfn,pcem,hl_gap,energy,gradient) + (iunit,env,opt,mol,gfn,pcem,wfn,hl_gap,energy,gradient) use iso_fortran_env, wp => real64 use mctc_systools @@ -378,6 +387,7 @@ subroutine gfn1_calculation & type(scc_options), intent(in) :: opt type(tb_environment), intent(in) :: env type(tb_pcem), intent(inout) :: pcem + type(tb_wavefunction),intent(inout) :: wfn real(wp), intent(out) :: energy real(wp), intent(out) :: hl_gap @@ -385,7 +395,6 @@ subroutine gfn1_calculation & integer, parameter :: wsc_rep(3) = [1,1,1] ! FIXME - type(tb_wavefunction) :: wfn type(tb_basisset) :: basis type(scc_parameter) :: param type(scc_results) :: res @@ -412,12 +421,12 @@ subroutine gfn1_calculation & ! we assume that the user provides a resonable molecule input ! -> all atoms are inside the unit cell, all data is set and consistent - ! get nuclear charges (without core) - call mol%set_nuclear_charge - call mol%set_atomic_masses - - wfn%nel = idint(sum(mol%z)) - mol%chrg - wfn%nopen = 0 + 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 @@ -455,8 +464,8 @@ subroutine gfn1_calculation & call gfn1_prparam(iunit,mol%n,mol%at,param) endif - if (len_trim(opt%solvent).gt.0 .and. opt%solvent.ne."none") then - lgbsa = .true. + 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) endif @@ -501,148 +510,6 @@ subroutine gfn1_calculation & end subroutine gfn1_calculation -! ======================================================================== -!> GFN0-xTB calculation -subroutine gfn0_calculation & - (iunit,env,opt,mol,gfn,hl_gap,energy,gradient) - use iso_fortran_env, wp => real64 - - use mctc_systools - - use tbdef_options - use tbdef_molecule - use tbdef_wavefunction - use tbdef_basisset - use tbdef_param - use tbdef_data - - use setparam, only : gfn_method - use aoparam, only : use_parameterset - use ehtparam, only : import_basisset - - use pbc_tools - use xbasis - use peeq_module - - implicit none - - integer, intent(in) :: iunit - - type(tb_molecule), intent(inout) :: mol - type(gfn_parameter), intent(in) :: gfn - type(peeq_options), intent(in) :: opt - type(tb_environment), intent(in) :: env - - real(wp), intent(out) :: energy - real(wp), intent(out) :: hl_gap - real(wp), intent(out) :: gradient(3,mol%n) - real(wp) :: 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(tb_wavefunction) :: wfn - type(tb_basisset) :: basis - type(scc_parameter) :: param - 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' - character(len=:), allocatable :: fnv - real(wp) :: globpar(25) - integer :: ipar - logical :: exist - - logical :: okbas,diff - - gfn_method = 0 - - ! ==================================================================== - ! 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 - - ! get nuclear charges (without core) - call mol%set_nuclear_charge - call mol%set_atomic_masses - - wfn%nel = idint(sum(mol%z)) - mol%chrg - wfn%nopen = 0 - - ! give an optional summary on the geometry used - if (opt%prlevel > 2) then - call print_pbcsum(iunit,mol) - call print_geosum(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 - - ! we will try an internal parameter file first to avoid IO - call use_parameterset(p_fnv_gfn0,globpar,exist) - ! no luck, we have to fire up some IO to get our parameters - if (.not.exist) then - ! let's check if we can find the parameter file - call rdpath(env%xtbpath,p_fnv_gfn0,fnv,exist) - ! 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_gfn0 - call open_file(ipar,fnv,'r') - if (ipar.eq.-1) then - ! at this point there is no chance to recover from this error - ! THEREFORE, we have to kill the program - call raise('E',"Parameter file '"//fnv//"' not found!",1) - return - endif - call read_gfn_param(ipar,globpar,.true.) - call close_file(ipar) - endif - call set_gfn0_parameter(param,globpar,mol%n,mol%at) - if (opt%prlevel > 1) then - call gfn0_header(iunit) - call gfn0_prparam(iunit,mol%n,mol%at,param) - endif - - ! ==================================================================== - ! STEP 3: expand our Slater basis set in contracted Gaussians - ! ==================================================================== - - call xbasis0(mol%n,mol%at,basis) - call xbasis_gfn0(mol%n,mol%at,basis,okbas,diff) - call xbasis_cao2sao(mol%n,mol%at,basis) - ! make the basis set available globally, this should be removed at some point - call import_basisset(basis) - - ! ==================================================================== - ! STEP 4: setup the initial wavefunction - ! ==================================================================== - - call wfn%allocate(mol%n,basis%nshell,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(iunit,mol,wfn,basis,param,hl_gap,opt%etemp,opt%prlevel,opt%grad, & - & .false.,opt%acc,energy,gradient,sigma,res) - - 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 gfn0_calculation - !> interface to the DFT-D4 module subroutine d4_calculation(iunit,opt,mol,dparam,energy,gradient) use iso_fortran_env, wp => real64 diff --git a/xtb/scf_module.f90 b/xtb/scf_module.f90 index b309653a9..6aaa1126d 100644 --- a/xtb/scf_module.f90 +++ b/xtb/scf_module.f90 @@ -437,18 +437,18 @@ subroutine scf(iunit,mol,wfn,basis,param,pcem, & if (profile) call timer%measure(2) if (profile) call timer%measure(3,"integral evaluation") + allocate(dpint(3,basis%nao*(basis%nao+1)/2), & + & qpint(6,basis%nao*(basis%nao+1)/2), & + & source = 0.0_wp) + ! compute integrals and prescreen to set up list arrays + call sdqint(mol%n,mol%at,basis%nbf,basis%nao,mol%xyz,neglect,ndp,nqp,intcut, & + & basis%caoshell,basis%saoshell,basis%nprim,basis%primcount, & + & basis%alp,basis%cont,S,dpint,qpint) + ! prepare aes stuff if(gfn_method.gt.1) then ! CN/dCN replaced by special smoother and faster decaying function call dncoord_gfn(mol%n,mol%at,mol%xyz,cn,dcn) - ii=basis%nao*(basis%nao+1)/2 - allocate(dpint(3,ii),qpint(6,ii)) - dpint=0.0_wp - qpint=0.0_wp -! compute integrals and prescreen to set up list arrays - call sdqint(mol%n,mol%at,basis%nbf,basis%nao,mol%xyz,neglect,ndp,nqp,intcut, & - & basis%caoshell,basis%saoshell,basis%nprim,basis%primcount, & - & basis%alp,basis%cont,S,dpint,qpint) ! allocate arrays for lists and fill (to exploit sparsity) allocate(mdlst(2,ndp),mqlst(2,nqp)) @@ -459,16 +459,6 @@ subroutine scf(iunit,mol,wfn,basis,param,pcem, & call get_radcn(mol%n,mol%at,cn,param%cn_shift,param%cn_expo,param%cn_rmax,radcn) call mmomgabzero(mol%n,mol%at,mol%xyz,param%xbrad,param%xbdamp,radcn,gab3,gab5) ! zero damping, xbrad=kdmp3,xbdamp=kdmp5 ! allocate CAMM arrays - else -! in GFN1 only AO overlap matrix needed - ii=basis%nao*(basis%nao+1)/2 - allocate(dpint(3,ii),qpint(6,ii),source = 0.0_wp) -! compute integrals and prescreen to set up list arrays - call sdqint(mol%n,mol%at,basis%nbf,basis%nao,mol%xyz,neglect,ndp,nqp,intcut, & - & basis%caoshell,basis%saoshell,basis%nprim,basis%primcount, & - & basis%alp,basis%cont,S,dpint,qpint) - deallocate(qpint) ! keep dipole integrals for later - endif if (profile) call timer%measure(3) @@ -477,7 +467,8 @@ subroutine scf(iunit,mol,wfn,basis,param,pcem, & if(gfn_method.gt.1) then ! if no CAMMs were read, get them from P (e.g., in geometry opts., MD runs) if(.not.restart) & - & call mmompop(mol%n,basis%nao,basis%aoat2,mol%xyz,wfn%p,s,dpint,qpint,wfn%dipm,wfn%qp) + & call mmompop(mol%n,basis%nao,basis%aoat2,mol%xyz,wfn%p,s,dpint,qpint, & + & wfn%dipm,wfn%qp) ! scale CAMMs before setting up the potential ! call scalecamm(mol%n,mol%at,dipm,qp) @@ -629,9 +620,6 @@ subroutine scf(iunit,mol,wfn,basis,param,pcem, & ! free some memory (this stuff is not needed for gradients) deallocate(ves) if (allocated(vpc)) deallocate(vpc) - if(gfn_method.gt.1) then - deallocate(qpint) ! potential stuff not needed any more - endif if(allocated(wdispmat)) deallocate( wdispmat ) 9999 continue @@ -727,9 +715,13 @@ subroutine scf(iunit,mol,wfn,basis,param,pcem, & endif printing - ! ------------------------------------------------------------------------ ! dipole calculation (always done because its free) + if (gfn_method.lt.2) then + call mmompop(mol%n,basis%nao,basis%aoat2,mol%xyz,wfn%p,s,dpint,qpint, & + & wfn%dipm,wfn%qp) + endif + call calc_dipole(mol%n,mol%at,mol%xyz,mol%z,basis%nao,wfn%P,dpint,dip,dipol) if (profile) call timer%measure(7) diff --git a/xtb/tbdef_molecule.f90 b/xtb/tbdef_molecule.f90 index c589e675a..d59cf02bf 100644 --- a/xtb/tbdef_molecule.f90 +++ b/xtb/tbdef_molecule.f90 @@ -34,6 +34,7 @@ module tbdef_molecule type :: tb_molecule integer :: n = 0 !< number of atoms real(wp) :: chrg = 0.0_wp !< total charge + integer :: uhf = 0 !< number of unpaired electrons logical :: pbc(3) = .false. !< periodic dimensions integer :: npbc = 0 !< periodicity of system character(len=2),allocatable :: sym(:) !< element symbols @@ -90,6 +91,7 @@ subroutine deallocate_molecule(self) self%n = 0 self%pbc = .false. self%chrg = 0.0_wp + self%uhf = 0 self%lattice = 0.0_wp if (allocated(self%at)) deallocate(self%at) if (allocated(self%sym)) deallocate(self%sym) @@ -116,6 +118,7 @@ subroutine write_molecule(self,iunit,comment) write(iunit,'(1x,"*",1x,a)') "status of the fields" write(iunit,dfmt) "integer :: n ",self%n write(iunit,dfmt) "real :: chrg ",self%chrg + write(iunit,dfmt) "integer :: uhf ",self%uhf write(iunit,dfmt) "integer :: npbc ",self%npbc write(iunit,dfmt) "logical :: pbc(1) ",self%pbc(1) write(iunit,dfmt) " & pbc(2) ",self%pbc(2) diff --git a/xtb/tbdef_options.f90 b/xtb/tbdef_options.f90 index 79f84fa12..f995b8c7b 100644 --- a/xtb/tbdef_options.f90 +++ b/xtb/tbdef_options.f90 @@ -81,6 +81,7 @@ module tbdef_options real(wp) :: etemp = 0.0_wp logical :: grad = .false. logical :: restart = .false. + logical :: ccm = .false. integer :: maxiter = 0 character(len=20) :: solvent = "none" end type scc_options @@ -92,6 +93,7 @@ module tbdef_options real(c_double) :: etemp = 0.0_c_double logical(c_bool) :: grad = .false._c_bool logical(c_bool) :: restart = .false._c_bool + logical(c_bool) :: ccm = .false._c_bool integer(c_int) :: maxiter = 0_c_int character(c_char) :: solvent(20) = c_null_char end type c_scc_options @@ -104,6 +106,7 @@ module tbdef_options real(wp) :: etemp = 0.0_wp logical :: grad = .false. logical :: ccm = .false. + character(len=20) :: solvent = "none" end type peeq_options type,bind(C) :: c_peeq_options @@ -113,6 +116,7 @@ module tbdef_options real(c_double) :: etemp = 0.0_c_double logical(c_bool) :: grad = .false._c_bool logical(c_bool) :: ccm = .false._c_bool + character(c_char) :: solvent(20) = c_null_char end type c_peeq_options type :: eeq_options @@ -292,6 +296,7 @@ pure elemental subroutine convert_scc_options_c_to_f & f_opt%grad = c_opt%grad f_opt%restart = c_opt%restart f_opt%parallel= c_opt%parallel + f_opt%ccm = c_opt%ccm f_opt%solvent = '' do i = 1, 20 @@ -315,6 +320,7 @@ pure elemental subroutine convert_scc_options_f_to_c & c_opt%grad = f_opt%grad c_opt%restart = f_opt%restart c_opt%parallel= f_opt%parallel + c_opt%ccm = f_opt%ccm do i = 1, len_trim(f_opt%solvent) c_opt%solvent(i) = f_opt%solvent(i:i) @@ -338,12 +344,19 @@ pure elemental subroutine convert_peeq_options_c_to_f & implicit none type(c_peeq_options), intent(in) :: c_opt type(peeq_options), intent(out) :: f_opt + integer :: i f_opt%prlevel = c_opt%prlevel f_opt%ccm = c_opt%ccm f_opt%acc = c_opt%acc f_opt%etemp = c_opt%etemp f_opt%grad = c_opt%grad f_opt%parallel= c_opt%parallel + + f_opt%solvent = '' + do i = 1, 20 + if (c_opt%solvent(i).eq.c_null_char) exit + f_opt%solvent(i:i) = c_opt%solvent(i) + enddo end subroutine convert_peeq_options_c_to_f pure elemental subroutine convert_peeq_options_f_to_c & @@ -351,12 +364,18 @@ pure elemental subroutine convert_peeq_options_f_to_c & implicit none type(peeq_options), intent(in) :: f_opt type(c_peeq_options),intent(out) :: c_opt + integer :: i c_opt%prlevel = f_opt%prlevel c_opt%ccm = f_opt%ccm c_opt%acc = f_opt%acc c_opt%etemp = f_opt%etemp c_opt%grad = f_opt%grad c_opt%parallel= f_opt%parallel + + do i = 1, len_trim(f_opt%solvent) + c_opt%solvent(i) = f_opt%solvent(i:i) + enddo + c_opt%solvent(len_trim(f_opt%solvent)+1) = c_null_char end subroutine convert_peeq_options_f_to_c pure elemental subroutine convert_eeq_options_c_to_f &