From 75b99213cea2afb1523401bde72a3446fbdf17ea Mon Sep 17 00:00:00 2001 From: Christian Blau Date: Fri, 13 Mar 2020 18:45:52 +0100 Subject: [PATCH] Remove dysfunctional QMMM interface pt2 Uncoupling from forcerec, deleting unused files. Change-Id: Id72281a26b8b39b6d2a37f5f8a43330f399ce5aa --- docs/user-guide/environment-variables.rst | 33 - src/gromacs/gmxpreprocess/grompp.cpp | 21 - src/gromacs/mdlib/force.cpp | 7 - src/gromacs/mdlib/forcerec.cpp | 27 +- src/gromacs/mdlib/mdatoms.cpp | 1 - src/gromacs/mdlib/qm_gamess.cpp | 254 ------ src/gromacs/mdlib/qm_gamess.h | 60 -- src/gromacs/mdlib/qm_gaussian.cpp | 894 ---------------------- src/gromacs/mdlib/qm_gaussian.h | 72 -- src/gromacs/mdlib/qm_mopac.cpp | 260 ------- src/gromacs/mdlib/qm_mopac.h | 67 -- src/gromacs/mdlib/qm_orca.cpp | 369 --------- src/gromacs/mdlib/qm_orca.h | 44 -- src/gromacs/mdlib/qmmm.cpp | 889 --------------------- src/gromacs/mdlib/qmmm.h | 152 ---- src/gromacs/mdlib/sim_util.cpp | 7 - src/gromacs/mdrun/runner.cpp | 2 +- src/gromacs/mdtypes/forcerec.h | 7 - 18 files changed, 4 insertions(+), 3162 deletions(-) delete mode 100644 src/gromacs/mdlib/qm_gamess.cpp delete mode 100644 src/gromacs/mdlib/qm_gamess.h delete mode 100644 src/gromacs/mdlib/qm_gaussian.cpp delete mode 100644 src/gromacs/mdlib/qm_gaussian.h delete mode 100644 src/gromacs/mdlib/qm_mopac.cpp delete mode 100644 src/gromacs/mdlib/qm_mopac.h delete mode 100644 src/gromacs/mdlib/qm_orca.cpp delete mode 100644 src/gromacs/mdlib/qm_orca.h delete mode 100644 src/gromacs/mdlib/qmmm.cpp delete mode 100644 src/gromacs/mdlib/qmmm.h diff --git a/docs/user-guide/environment-variables.rst b/docs/user-guide/environment-variables.rst index 8325547968..4a19a2e208 100644 --- a/docs/user-guide/environment-variables.rst +++ b/docs/user-guide/environment-variables.rst @@ -4,8 +4,6 @@ .. Another useful one-liner to find undocumentedvariables: .. ( export INPUT_FILE=docs/user-guide/environment-variables.rst; GIT_PAGER="cat "; for ss in `for s in $(git grep getenv | sed 's/.*getenv("\(.*\)".*/\1/' | sort -u | grep '^[A-Z]'); do [ $(grep $s $INPUT_FILE -c) -eq 0 ] && echo $s; done `; do git grep $ss ; done ) -.. todo:: still undocumented GMX_QM_GAUSSIAN_NCPUS - Environment Variables ===================== @@ -499,30 +497,11 @@ compilation of OpenCL kernels, but they are also used in device selection. Analysis and Core Functions --------------------------- -``GMX_QM_ACCURACY`` - accuracy in Gaussian L510 (MC-SCF) component program. - -``GMX_QM_ORCA_BASENAME`` - prefix of :ref:`tpr` files, used in Orca calculations - for input and output file names. - -``GMX_QM_CPMCSCF`` - when set to a nonzero value, Gaussian QM calculations will - iteratively solve the CP-MCSCF equations. - -``GMX_QM_MODIFIED_LINKS_DIR`` - location of modified links in Gaussian. ``DSSP`` used by :ref:`gmx do_dssp` to point to the ``dssp`` executable (not just its path). -``GMX_QM_GAUSS_DIR`` - directory where Gaussian is installed. - -``GMX_QM_GAUSS_EXE`` - name of the Gaussian executable. - ``GMX_DIPOLE_SPACING`` spacing used by :ref:`gmx dipoles`. @@ -545,9 +524,6 @@ Analysis and Core Functions the time unit used in output files, can be anything in fs, ps, ns, us, ms, s, m or h. -``GMX_QM_GAUSSIAN_MEMORY`` - memory used for Gaussian QM calculation. - ``MULTIPROT`` name of the ``multiprot`` executable, used by the contributed program ``do_multiprot``. @@ -555,15 +531,6 @@ Analysis and Core Functions ``NCPUS`` number of CPUs to be used for Gaussian QM calculation -``GMX_ORCA_PATH`` - directory where Orca is installed. - -``GMX_QM_SA_STEP`` - simulated annealing step size for Gaussian QM calculation. - -``GMX_QM_GROUND_STATE`` - defines state for Gaussian surface hopping calculation. - ``GMX_TOTAL`` name of the ``total`` executable used by the contributed ``do_shift`` program. diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 0ba7a763c5..20ad419a74 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -80,7 +80,6 @@ #include "gromacs/mdlib/compute_io.h" #include "gromacs/mdlib/constr.h" #include "gromacs/mdlib/perf_est.h" -#include "gromacs/mdlib/qmmm.h" #include "gromacs/mdlib/vsite.h" #include "gromacs/mdrun/mdmodules.h" #include "gromacs/mdtypes/inputrec.h" @@ -1771,7 +1770,6 @@ int gmx_grompp(int argc, char* argv[]) double reppow; const char* mdparin; bool bNeedVel, bGenVel; - gmx_bool have_atomnumber; gmx_output_env_t* oenv; gmx_bool bVerbose = FALSE; warninp* wi; @@ -1959,25 +1957,6 @@ int gmx_grompp(int argc, char* argv[]) warning_note(wi, "Temperature coupling is ignored with SD integrators."); } - /* If we are doing QM/MM, check that we got the atom numbers */ - have_atomnumber = TRUE; - for (gmx::index i = 0; i < gmx::ssize(atypes); i++) - { - have_atomnumber = have_atomnumber && (atypes.atomNumberFromAtomType(i) >= 0); - } - if (!have_atomnumber && ir->bQMMM) - { - warning_error( - wi, - "\n" - "It appears as if you are trying to run a QM/MM calculation, but the force\n" - "field you are using does not contain atom numbers fields. This is an\n" - "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n" - "for QM/MM. The good news is that it is easy to add - put the atom number as\n" - "an integer just before the mass column in ffXXXnb.itp.\n" - "NB: United atoms have the same atom numbers as normal ones.\n\n"); - } - /* Check for errors in the input now, since they might cause problems * during processing further down. */ diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index 6a6a124895..671343c09b 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -55,7 +55,6 @@ #include "gromacs/math/vec.h" #include "gromacs/math/vecdump.h" #include "gromacs/mdlib/forcerec_threading.h" -#include "gromacs/mdlib/qmmm.h" #include "gromacs/mdlib/rf_util.h" #include "gromacs/mdlib/wall.h" #include "gromacs/mdtypes/commrec.h" @@ -125,12 +124,6 @@ void do_force_lowlevel(t_forcerec* fr, auto& forceWithVirial = forceOutputs->forceWithVirial(); - /* do QMMM first if requested */ - if (fr->bQMMM) - { - enerd->term[F_EQM] = calculate_QMMM(cr, &forceOutputs->forceWithShiftForces(), fr->qr); - } - /* Call the short range functions all in one go. */ if (ir->nwall) diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index e2eeddd818..aa4c7b540b 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -69,7 +69,6 @@ #include "gromacs/mdlib/forcerec_threading.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdlib/md_support.h" -#include "gromacs/mdlib/qmmm.h" #include "gromacs/mdlib/rf_util.h" #include "gromacs/mdlib/wall.h" #include "gromacs/mdtypes/commrec.h" @@ -1327,30 +1326,10 @@ void init_forcerec(FILE* fp, } // QM/MM initialization if requested - fr->bQMMM = ir->bQMMM; - if (fr->bQMMM) + if (ir->bQMMM) { - // Initialize QM/MM if supported - if (GMX_QMMM) - { - GMX_LOG(mdlog.info) - .asParagraph() - .appendText( - "Large parts of the QM/MM support is deprecated, and may be removed in " - "a future " - "version. Please get in touch with the developers if you find the " - "support useful, " - "as help is needed if the functionality is to continue to be " - "available."); - fr->qr = mk_QMMMrec(); - init_QMMMrec(cr, mtop, ir, fr); - } - else - { - gmx_incons( - "QM/MM was requested, but is only available when GROMACS " - "is configured with QM/MM support"); - } + + gmx_incons("QM/MM was requested, but is no longer available in GROMACS"); } /* Set all the static charge group info */ diff --git a/src/gromacs/mdlib/mdatoms.cpp b/src/gromacs/mdlib/mdatoms.cpp index 3dc2cf61a3..1426046998 100644 --- a/src/gromacs/mdlib/mdatoms.cpp +++ b/src/gromacs/mdlib/mdatoms.cpp @@ -47,7 +47,6 @@ #include "gromacs/gpu_utils/hostallocator.h" #include "gromacs/math/functions.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" -#include "gromacs/mdlib/qmmm.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" #include "gromacs/mdtypes/mdatom.h" diff --git a/src/gromacs/mdlib/qm_gamess.cpp b/src/gromacs/mdlib/qm_gamess.cpp deleted file mode 100644 index 4b86adeeeb..0000000000 --- a/src/gromacs/mdlib/qm_gamess.cpp +++ /dev/null @@ -1,254 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gmxpre.h" - -#include "qm_gamess.h" - -#include "config.h" - -#include -#include -#include - -#include - -#include "gromacs/fileio/confio.h" -#include "gromacs/gmxlib/network.h" -#include "gromacs/gmxlib/nrnb.h" -#include "gromacs/math/units.h" -#include "gromacs/math/vec.h" -#include "gromacs/mdlib/qmmm.h" -#include "gromacs/mdtypes/commrec.h" -#include "gromacs/mdtypes/forcerec.h" -#include "gromacs/mdtypes/md_enums.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/smalloc.h" - -// When not built in a configuration with QMMM support, much of this -// code is unreachable by design. Tell clang not to warn about it. -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wmissing-noreturn" - -/* QMMM sub routines */ -/* mopac interface routines */ - - -static void F77_FUNC(inigms, IMIGMS)(); - -static void F77_FUNC(grads, GRADS)(const int* nrqmat, - real* qmcrd, - const int* nrmmat, - const real* mmchrg, - real* mmcrd, - real* qmgrad, - real* mmgrad, - real* energy); - -#if !GMX_QMMM_GAMESS -// Stub definitions to make compilation succeed when not configured -// for GAMESS support. In that case, the module gives a fatal error -// when the initialization function is called, so there is no need to -// issue fatal errors here, because that introduces problems with -// tools suggesting and prohibiting noreturn attributes. - -void F77_FUNC(inigms, IMIGMS)(){}; -// NOLINTNEXTLINE(readability-named-parameter) -void F77_FUNC(grads, GRADS)(const int*, real*, const int*, const real*, real*, real*, real*, real*){}; -#endif - -void init_gamess(const t_commrec* cr, t_QMrec* qm, t_MMrec* mm) -{ - /* it works hopelessly complicated :-) - * first a file is written. Then the standard gamess input/output - * routine is called (no system()!) to set up all fortran arrays. - * this routine writes a punch file, like in a normal gamess run. - * via this punch file the other games routines, needed for gradient - * and energy evaluations are called. This setup works fine for - * dynamics simulations. 7-6-2002 (London) - */ - int i, j; - FILE* out; - char periodic_system[37][3] = { "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ", "O ", "F ", - "Ne", "Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar", "K ", - "Ca", "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", - "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr" }; - if (!GMX_QMMM_GAMESS) - { - gmx_fatal(FARGS, - "Cannot call GAMESS unless linked against it. Use cmake " - "-DGMX_QMMM_PROGRAM=GAMESS, and ensure that linking will work correctly."); - } - - if (PAR(cr)) - { - - if (MASTER(cr)) - { - out = fopen("FOR009", "w"); - /* of these options I am not completely sure.... the overall - * preformance on more than 4 cpu's is rather poor at the moment. - */ - fprintf(out, "memory 48000000\nPARALLEL IOMODE SCREENED\n"); - fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n", qm->nelectrons, - qm->multiplicity); - for (i = 0; i < qm->nrQMatoms; i++) - { -#ifdef DOUBLE - fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n", i / 2., i / 3., i / 4., - qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]); -#else - fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n", i / 2., i / 3., i / 4., - qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]); -#endif - } - if (mm->nrMMatoms) - { - for (j = i; j < i + 2; j++) - { -#ifdef DOUBLE - fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n", j / 5., j / 6., j / 7., 1.0); -#else - fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n", j / 5., j / 6., j / 7., 2.0); -#endif - } - } - fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n", - eQMbasis_names[qm->QMbasis], eQMmethod_names[qm->QMmethod]); /* see enum.h */ - fclose(out); - } - gmx_barrier(cr); - F77_FUNC(inigms, IMIGMS)(); - } - else /* normal serial run */ - - { - out = fopen("FOR009", "w"); - /* of these options I am not completely sure.... the overall - * preformance on more than 4 cpu's is rather poor at the moment. - */ - fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n", qm->nelectrons, - qm->multiplicity); - for (i = 0; i < qm->nrQMatoms; i++) - { -#ifdef DOUBLE - fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n", i / 2., i / 3., i / 4., - qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]); -#else - fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n", i / 2., i / 3., i / 4., - qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]); -#endif - } - if (mm->nrMMatoms) - { - for (j = i; j < i + 2; j++) - { -#ifdef DOUBLE - fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n", j / 5., j / 6., j / 7., 1.0); -#else - fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n", j / 5., j / 6., j / 7., 2.0); -#endif - } - } - fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n", eQMbasis_names[qm->QMbasis], - eQMmethod_names[qm->QMmethod]); /* see enum.h */ - F77_FUNC(inigms, IMIGMS)(); - } -} - -real call_gamess(const t_QMrec* qm, const t_MMrec* mm, rvec f[], rvec fshift[]) -{ - /* do the actual QMMM calculation using GAMESS-UK. In this - * implementation (3-2001) a system call is made to the GAMESS-UK - * binary. Now we are working to get the electron integral, SCF, and - * gradient routines linked directly - */ - int i, j; - real QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy = 0; - - snew(qmcrd, 3 * (qm->nrQMatoms)); - snew(mmcrd, 3 * (mm->nrMMatoms)); - snew(qmgrad, 3 * (qm->nrQMatoms)); - snew(mmgrad, 3 * (mm->nrMMatoms)); - - /* copy the data from qr into the arrays that are going to be used - * in the fortran routines of gamess - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - qmcrd[DIM * i + j] = 1 / BOHR2NM * qm->xQM[i][j]; - } - } - for (i = 0; i < mm->nrMMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - mmcrd[DIM * i + j] = 1 / BOHR2NM * mm->xMM[i][j]; - } - } - for (i = 0; i < 3 * qm->nrQMatoms; i += 3) - { - fprintf(stderr, "%8.5f, %8.5f, %8.5f\n", qmcrd[i], qmcrd[i + 1], qmcrd[i + 2]); - } - - F77_FUNC(grads, GRADS) - (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mm->MMcharges, mmcrd, qmgrad, mmgrad, &energy); - - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i][j] = HARTREE_BOHR2MD * qmgrad[3 * i + j]; - fshift[i][j] = HARTREE_BOHR2MD * qmgrad[3 * i + j]; - } - } - for (i = 0; i < mm->nrMMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i][j] = HARTREE_BOHR2MD * mmgrad[3 * i + j]; - fshift[i][j] = HARTREE_BOHR2MD * mmgrad[3 * i + j]; - } - } - /* convert a.u to kJ/mol */ - QMener = energy * HARTREE2KJ * AVOGADRO; - return (QMener); -} - -#pragma GCC diagnostic pop diff --git a/src/gromacs/mdlib/qm_gamess.h b/src/gromacs/mdlib/qm_gamess.h deleted file mode 100644 index 81c2b4c40e..0000000000 --- a/src/gromacs/mdlib/qm_gamess.h +++ /dev/null @@ -1,60 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2018,2019, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#ifndef GMX_MDLIB_QMGAMESS_H -#define GMX_MDLIB_QMGAMESS_H - -#include "gromacs/mdlib/qmmm.h" - -/*! \brief - * Initialize gamess datastructures. - * - * \param[in] cr Commrec datastructure. - * \param[in] qm QM forcerec. - * \param[in] mm MM part of forcerec. - */ -void init_gamess(const t_commrec* cr, t_QMrec* qm, t_MMrec* mm); - -/*! \brief - * Run calculation with Gamess. - * - * \param[in] qm QM part of forcerec. - * \param[in] mm MM part of forcerec. - * \param[in] f Force vector. - * \param[in] fshift Force shift vector. - */ -real call_gamess(const t_QMrec* qm, const t_MMrec* mm, rvec f[], rvec fshift[]); - - -#endif diff --git a/src/gromacs/mdlib/qm_gaussian.cpp b/src/gromacs/mdlib/qm_gaussian.cpp deleted file mode 100644 index 5e648e59d8..0000000000 --- a/src/gromacs/mdlib/qm_gaussian.cpp +++ /dev/null @@ -1,894 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gmxpre.h" - -#include "qm_gaussian.h" - -#include "config.h" - -#include -#include -#include - -#include - -#include "gromacs/fileio/confio.h" -#include "gromacs/gmxlib/network.h" -#include "gromacs/gmxlib/nrnb.h" -#include "gromacs/math/units.h" -#include "gromacs/math/vec.h" -#include "gromacs/mdlib/force.h" -#include "gromacs/mdlib/qmmm.h" -#include "gromacs/mdtypes/md_enums.h" -#include "gromacs/utility/cstringutil.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/smalloc.h" - -// When not built in a configuration with QMMM support, much of this -// code is unreachable by design. Tell clang not to warn about it. -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wmissing-noreturn" - -/* TODO: this should be made thread-safe */ - -/* Gaussian interface routines */ - -void init_gaussian(t_QMrec* qm) -{ - ivec basissets[eQMbasisNR] = { { 0, 3, 0 }, { 0, 3, 0 }, /*added for double sto-3g entry in names.c*/ - { 5, 0, 0 }, { 5, 0, 1 }, { 5, 0, 11 }, { 5, 6, 0 }, - { 1, 6, 0 }, { 1, 6, 1 }, { 1, 6, 11 }, { 4, 6, 0 } }; - char* buf = nullptr; - int i; - - if (!GMX_QMMM_GAUSSIAN) - { - gmx_fatal(FARGS, - "Cannot call GAUSSIAN unless linked against it. Use cmake " - "-DGMX_QMMM_PROGRAM=GAUSSIAN, and ensure that linking will work correctly."); - } - - /* using the ivec above to convert the basis read form the mdp file - * in a human readable format into some numbers for the gaussian - * route. This is necessary as we are using non standard routes to - * do SH. - */ - - /* per layer we make a new subdir for integral file, checkpoint - * files and such. These dirs are stored in the QMrec for - * convenience - */ - - - if (!qm->nQMcpus) /* this we do only once per layer - * as we call g01 externally - */ - - { - for (i = 0; i < DIM; i++) - { - qm->SHbasis[i] = basissets[qm->QMbasis][i]; - } - - /* init gradually switching on of the SA */ - qm->SAstep = 0; - /* we read the number of cpus and environment from the environment - * if set. - */ - buf = getenv("GMX_QM_GAUSSIAN_NCPUS"); - if (buf) - { - sscanf(buf, "%d", &qm->nQMcpus); - } - else - { - qm->nQMcpus = 1; - } - fprintf(stderr, "number of CPUs for gaussian = %d\n", qm->nQMcpus); - buf = getenv("GMX_QM_GAUSSIAN_MEMORY"); - if (buf) - { - sscanf(buf, "%d", &qm->QMmem); - } - else - { - qm->QMmem = 50000000; - } - fprintf(stderr, "memory for gaussian = %d\n", qm->QMmem); - buf = getenv("GMX_QM_ACCURACY"); - if (buf) - { - sscanf(buf, "%d", &qm->accuracy); - } - else - { - qm->accuracy = 8; - } - fprintf(stderr, "accuracy in l510 = %d\n", qm->accuracy); - - buf = getenv("GMX_QM_CPMCSCF"); - if (buf) - { - sscanf(buf, "%d", &i); - qm->cpmcscf = (i != 0); - } - else - { - qm->cpmcscf = FALSE; - } - if (qm->cpmcscf) - { - fprintf(stderr, "using cp-mcscf in l1003\n"); - } - else - { - fprintf(stderr, "NOT using cp-mcscf in l1003\n"); - } - buf = getenv("GMX_QM_SA_STEP"); - if (buf) - { - sscanf(buf, "%d", &qm->SAstep); - } - else - { - /* init gradually switching on of the SA */ - qm->SAstep = 0; - } - /* we read the number of cpus and environment from the environment - * if set. - */ - fprintf(stderr, "Level of SA at start = %d\n", qm->SAstep); - /* gaussian settings on the system */ - buf = getenv("GMX_QM_GAUSS_DIR"); - - if (buf) - { - qm->gauss_dir = gmx_strdup(buf); - } - else - { - gmx_fatal(FARGS, "no $GMX_QM_GAUSS_DIR, check gaussian manual\n"); - } - - buf = getenv("GMX_QM_GAUSS_EXE"); - if (buf) - { - qm->gauss_exe = gmx_strdup(buf); - } - else - { - gmx_fatal(FARGS, "no $GMX_QM_GAUSS_EXE set, check gaussian manual\n"); - } - buf = getenv("GMX_QM_MODIFIED_LINKS_DIR"); - if (buf) - { - qm->devel_dir = gmx_strdup(buf); - } - else - { - gmx_fatal(FARGS, - "no $GMX_QM_MODIFIED_LINKS_DIR, this is were the modified links reside.\n"); - } - - /* if(fr->bRF){*/ - /* reactionfield, file is needed using gaussian */ - /* rffile=fopen("rf.dat","w");*/ - /* fprintf(rffile,"%f %f\n",fr->epsilon_r,fr->rcoulomb/BOHR2NM);*/ - /* fclose(rffile);*/ - /* }*/ - } - fprintf(stderr, "gaussian initialised...\n"); -} - - -static void write_gaussian_SH_input(int step, gmx_bool swap, const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm) -{ - int i; - bool bSA = (qm->SAstep > 0); - FILE* out = fopen("input.com", "w"); - /* write the route */ - fprintf(out, "%s", "%scr=input\n"); - fprintf(out, "%s", "%rwf=input\n"); - fprintf(out, "%s", "%int=input\n"); - fprintf(out, "%s", "%d2e=input\n"); - /* if(step) - * fprintf(out,"%s","%nosave\n"); - */ - fprintf(out, "%s", "%chk=input\n"); - fprintf(out, "%s%d\n", "%mem=", qm->QMmem); - fprintf(out, "%s%3d\n", "%nprocshare=", qm->nQMcpus); - - /* use the versions of - * l301 that computes the interaction between MM and QM atoms. - * l510 that can punch the CI coefficients - * l701 that can do gradients on MM atoms - */ - - /* local version */ - fprintf(out, "%s%s%s", "%subst l510 ", qm->devel_dir, "/l510\n"); - fprintf(out, "%s%s%s", "%subst l301 ", qm->devel_dir, "/l301\n"); - fprintf(out, "%s%s%s", "%subst l701 ", qm->devel_dir, "/l701\n"); - - fprintf(out, "%s%s%s", "%subst l1003 ", qm->devel_dir, "/l1003\n"); - fprintf(out, "%s%s%s", "%subst l9999 ", qm->devel_dir, "/l9999\n"); - /* print the nonstandard route - */ - fprintf(out, "%s", "#P nonstd\n 1/18=10,20=1,38=1/1;\n"); - fprintf(out, "%s", " 2/9=110,15=1,17=6,18=5,40=1/2;\n"); - if (mm->nrMMatoms) - { - fprintf(out, " 3/5=%d,6=%d,7=%d,25=1,32=1,43=1,94=-2/1,2,3;\n", qm->SHbasis[0], - qm->SHbasis[1], qm->SHbasis[2]); /*basisset stuff */ - } - else - { - fprintf(out, " 3/5=%d,6=%d,7=%d,25=1,32=1,43=0,94=-2/1,2,3;\n", qm->SHbasis[0], - qm->SHbasis[1], qm->SHbasis[2]); /*basisset stuff */ - } - /* development */ - if (step + 1) /* fetch initial guess from check point file */ - { /* hack, to alyays read from chk file!!!!! */ - fprintf(out, "%s%d,%s%d%s", " 4/5=1,7=6,17=", qm->CASelectrons, "18=", qm->CASorbitals, - "/1,5;\n"); - } - else /* generate the first checkpoint file */ - { - fprintf(out, "%s%d,%s%d%s", " 4/5=0,7=6,17=", qm->CASelectrons, "18=", qm->CASorbitals, - "/1,5;\n"); - } - /* the rest of the input depends on where the system is on the PES - */ - if (swap && bSA) /* make a slide to the other surface */ - { - if (qm->CASorbitals > 6) /* use direct and no full diag */ - { - fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6,97=100/10;\n"); - } - else - { - if (qm->cpmcscf) - { - fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6,97=100/10;\n", qm->accuracy); - if (mm->nrMMatoms > 0) - { - fprintf(out, " 7/7=1,16=-2,30=1/1;\n"); - } - fprintf(out, " 11/31=1,42=1,45=1/1;\n"); - fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1,97=100/3;\n"); - fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n"); - } - else - { - fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6,97=100/10;\n", qm->accuracy); - fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n"); - } - } - } - else if (bSA) /* do a "state-averaged" CAS calculation */ - { - if (qm->CASorbitals > 6) /* no full diag */ - { - fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6/10;\n"); - } - else - { - if (qm->cpmcscf) - { - fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6/10;\n", qm->accuracy); - if (mm->nrMMatoms > 0) - { - fprintf(out, " 7/7=1,16=-2,30=1/1;\n"); - } - fprintf(out, " 11/31=1,42=1,45=1/1;\n"); - fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1/3;\n"); - fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n"); - } - else - { - fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6/10;\n", qm->accuracy); - fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n"); - } - } - } - else if (swap) /* do a "swapped" CAS calculation */ - { - if (qm->CASorbitals > 6) - { - fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6,97=100/10;\n"); - } - else - { - fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6,97=100/10;\n", qm->accuracy); - } - fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n"); - } - else /* do a "normal" CAS calculation */ - { - if (qm->CASorbitals > 6) - { - fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6/10;\n"); - } - else - { - fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6/10;\n", qm->accuracy); - } - fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n"); - } - fprintf(out, "\ninput-file generated by gromacs\n\n"); - fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity); - for (i = 0; i < qm->nrQMatoms; i++) - { - fprintf(out, "%3d %10.7f %10.7f %10.7f\n", qm->atomicnumberQM[i], - qm->xQM[i][XX] / BOHR2NM, qm->xQM[i][YY] / BOHR2NM, qm->xQM[i][ZZ] / BOHR2NM); - } - /* MM point charge data */ - if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) - { - fprintf(out, "\n"); - for (i = 0; i < mm->nrMMatoms; i++) - { - fprintf(out, "%10.7f %10.7f %10.7f %8.4f\n", mm->xMM[i][XX] / BOHR2NM, - mm->xMM[i][YY] / BOHR2NM, mm->xMM[i][ZZ] / BOHR2NM, mm->MMcharges[i]); - } - } - if (bSA) /* put the SA coefficients at the end of the file */ - { - fprintf(out, "\n%10.8f %10.8f\n", qm->SAstep * 0.5 / qm->SAsteps, 1 - qm->SAstep * 0.5 / qm->SAsteps); - fprintf(stderr, "State Averaging level = %d/%d\n", qm->SAstep, qm->SAsteps); - } - fprintf(out, "\n"); - fclose(out); -} /* write_gaussian_SH_input */ - -static void write_gaussian_input(int step, const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm) -{ - int i; - - FILE* out = fopen("input.com", "w"); - /* write the route */ - - if (qm->QMmethod >= eQMmethodRHF) - { - fprintf(out, "%s", "%chk=input\n"); - } - else - { - fprintf(out, "%s", "%chk=se\n"); - } - if (qm->nQMcpus > 1) - { - fprintf(out, "%s%3d\n", "%nprocshare=", qm->nQMcpus); - } - fprintf(out, "%s%d\n", "%mem=", qm->QMmem); - fprintf(out, "%s%s%s", "%subst l701 ", qm->devel_dir, "/l701\n"); - fprintf(out, "%s%s%s", "%subst l301 ", qm->devel_dir, "/l301\n"); - fprintf(out, "%s%s%s", "%subst l9999 ", qm->devel_dir, "/l9999\n"); - if (step) - { - fprintf(out, "%s", "#T "); - } - else - { - fprintf(out, "%s", "#P "); - } - if (qm->QMmethod == eQMmethodB3LYPLAN) - { - fprintf(out, " %s", "B3LYP/GEN Pseudo=Read"); - } - else - { - fprintf(out, " %s", eQMmethod_names[qm->QMmethod]); - - if (qm->QMmethod >= eQMmethodRHF) - { - if (qm->QMmethod == eQMmethodCASSCF) - { - /* in case of cas, how many electrons and orbitals do we need? - */ - fprintf(out, "(%d,%d)", qm->CASelectrons, qm->CASorbitals); - } - fprintf(out, "/%s", eQMbasis_names[qm->QMbasis]); - } - } - if (QMMMrec->QMMMscheme == eQMMMschemenormal && mm->nrMMatoms) - { - fprintf(out, " %s", "Charge "); - } - if (step || qm->QMmethod == eQMmethodCASSCF) - { - /* fetch guess from checkpoint file, always for CASSCF */ - fprintf(out, "%s", " guess=read"); - } - fprintf(out, "\nNosymm units=bohr\n"); - - fprintf(out, "FORCE Punch=(Derivatives) "); - fprintf(out, "iop(3/33=1)\n\n"); - fprintf(out, "input-file generated by gromacs\n\n"); - fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity); - for (i = 0; i < qm->nrQMatoms; i++) - { - fprintf(out, "%3d %10.7f %10.7f %10.7f\n", qm->atomicnumberQM[i], - qm->xQM[i][XX] / BOHR2NM, qm->xQM[i][YY] / BOHR2NM, qm->xQM[i][ZZ] / BOHR2NM); - } - - /* Pseudo Potential and ECP are included here if selected (MEthod suffix LAN) */ - if (qm->QMmethod == eQMmethodB3LYPLAN) - { - fprintf(out, "\n"); - for (i = 0; i < qm->nrQMatoms; i++) - { - if (qm->atomicnumberQM[i] < 21) - { - fprintf(out, "%d ", i + 1); - } - } - fprintf(out, "\n%s\n****\n", eQMbasis_names[qm->QMbasis]); - - for (i = 0; i < qm->nrQMatoms; i++) - { - if (qm->atomicnumberQM[i] > 21) - { - fprintf(out, "%d ", i + 1); - } - } - fprintf(out, "\n%s\n****\n\n", "lanl2dz"); - - for (i = 0; i < qm->nrQMatoms; i++) - { - if (qm->atomicnumberQM[i] > 21) - { - fprintf(out, "%d ", i + 1); - } - } - fprintf(out, "\n%s\n", "lanl2dz"); - } - - - /* MM point charge data */ - if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) - { - fprintf(stderr, "nr mm atoms in gaussian.c = %d\n", mm->nrMMatoms); - fprintf(out, "\n"); - for (i = 0; i < mm->nrMMatoms; i++) - { - fprintf(out, "%10.7f %10.7f %10.7f %8.4f\n", mm->xMM[i][XX] / BOHR2NM, - mm->xMM[i][YY] / BOHR2NM, mm->xMM[i][ZZ] / BOHR2NM, mm->MMcharges[i]); - } - } - fprintf(out, "\n"); - - - fclose(out); - -} /* write_gaussian_input */ - -static real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], t_QMrec* qm, t_MMrec* mm) -{ - int i; - char buf[300]; - real QMener; - FILE* in; - - in = fopen("fort.7", "r"); - - /* (There was additional content in the file in case - * of QM optimizations / transition state search, - * which was removed. - */ - /* the next line is the energy and in the case of CAS, the energy - * difference between the two states. - */ - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } - -#if GMX_DOUBLE - sscanf(buf, "%lf\n", &QMener); -#else - sscanf(buf, "%f\n", &QMener); -#endif - /* next lines contain the gradients of the QM atoms */ - for (i = 0; i < qm->nrQMatoms; i++) - { - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } -#if GMX_DOUBLE - sscanf(buf, "%lf %lf %lf\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]); -#else - sscanf(buf, "%f %f %f\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]); -#endif - } - /* the next lines are the gradients of the MM atoms */ - if (qm->QMmethod >= eQMmethodRHF) - { - for (i = 0; i < mm->nrMMatoms; i++) - { - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } -#if GMX_DOUBLE - sscanf(buf, "%lf %lf %lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]); -#else - sscanf(buf, "%f %f %f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]); -#endif - } - } - fclose(in); - return (QMener); -} - -static real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step, t_QMrec* qm, t_MMrec* mm) -{ - int i; - char buf[300]; - real QMener, DeltaE; - FILE* in; - - in = fopen("fort.7", "r"); - /* first line is the energy and in the case of CAS, the energy - * difference between the two states. - */ - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } - -#if GMX_DOUBLE - sscanf(buf, "%lf %lf\n", &QMener, &DeltaE); -#else - sscanf(buf, "%f %f\n", &QMener, &DeltaE); -#endif - - /* switch on/off the State Averaging */ - - if (DeltaE > qm->SAoff) - { - if (qm->SAstep > 0) - { - qm->SAstep--; - } - } - else if (DeltaE < qm->SAon || (qm->SAstep > 0)) - { - if (qm->SAstep < qm->SAsteps) - { - qm->SAstep++; - } - } - - /* for debugging: */ - fprintf(stderr, "Gap = %5f,SA = %3d\n", DeltaE, static_cast(qm->SAstep > 0)); - /* next lines contain the gradients of the QM atoms */ - for (i = 0; i < qm->nrQMatoms; i++) - { - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } - -#if GMX_DOUBLE - sscanf(buf, "%lf %lf %lf\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]); -#else - sscanf(buf, "%f %f %f\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]); -#endif - } - /* the next lines, are the gradients of the MM atoms */ - - for (i = 0; i < mm->nrMMatoms; i++) - { - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } -#if GMX_DOUBLE - sscanf(buf, "%lf %lf %lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]); -#else - sscanf(buf, "%f %f %f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]); -#endif - } - - /* the next line contains the two CI eigenvector elements */ - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } - if (!step) - { - sscanf(buf, "%d", &qm->CIdim); - snew(qm->CIvec1, qm->CIdim); - snew(qm->CIvec1old, qm->CIdim); - snew(qm->CIvec2, qm->CIdim); - snew(qm->CIvec2old, qm->CIdim); - } - else - { - /* before reading in the new current CI vectors, copy the current - * CI vector into the old one. - */ - for (i = 0; i < qm->CIdim; i++) - { - qm->CIvec1old[i] = qm->CIvec1[i]; - qm->CIvec2old[i] = qm->CIvec2[i]; - } - } - /* first vector */ - for (i = 0; i < qm->CIdim; i++) - { - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } -#if GMX_DOUBLE - sscanf(buf, "%lf\n", &qm->CIvec1[i]); -#else - sscanf(buf, "%f\n", &qm->CIvec1[i]); -#endif - } - /* second vector */ - for (i = 0; i < qm->CIdim; i++) - { - if (nullptr == fgets(buf, 300, in)) - { - gmx_fatal(FARGS, "Error reading Gaussian output"); - } -#if GMX_DOUBLE - sscanf(buf, "%lf\n", &qm->CIvec2[i]); -#else - sscanf(buf, "%f\n", &qm->CIvec2[i]); -#endif - } - fclose(in); - return (QMener); -} - -static real inproduct(const real* a, const real* b, int n) -{ - int i; - real dot = 0.0; - - /* computes the inner product between two vectors (a.b), both of - * which have length n. - */ - for (i = 0; i < n; i++) - { - dot += a[i] * b[i]; - } - return (dot); -} - -static int hop(int step, t_QMrec* qm) -{ - int swap = 0; - real d11 = 0.0, d12 = 0.0, d21 = 0.0, d22 = 0.0; - - /* calculates the inproduct between the current Ci vector and the - * previous CI vector. A diabatic hop will be made if d12 and d21 - * are much bigger than d11 and d22. In that case hop returns true, - * otherwise it returns false. - */ - if (step) /* only go on if more than one step has been done */ - { - d11 = inproduct(qm->CIvec1, qm->CIvec1old, qm->CIdim); - d12 = inproduct(qm->CIvec1, qm->CIvec2old, qm->CIdim); - d21 = inproduct(qm->CIvec2, qm->CIvec1old, qm->CIdim); - d22 = inproduct(qm->CIvec2, qm->CIvec2old, qm->CIdim); - } - fprintf(stderr, "-------------------\n"); - fprintf(stderr, "d11 = %13.8f\n", d11); - fprintf(stderr, "d12 = %13.8f\n", d12); - fprintf(stderr, "d21 = %13.8f\n", d21); - fprintf(stderr, "d22 = %13.8f\n", d22); - fprintf(stderr, "-------------------\n"); - - if ((fabs(d12) > 0.5) && (fabs(d21) > 0.5)) - { - swap = 1; - } - - return (swap); -} - -static void do_gaussian(int step, char* exe) -{ - char buf[STRLEN]; - - /* make the call to the gaussian binary through system() - * The location of the binary will be picked up from the - * environment using getenv(). - */ - if (step) /* hack to prevent long inputfiles */ - { - sprintf(buf, "%s < %s > %s", exe, "input.com", "input.log"); - } - else - { - sprintf(buf, "%s < %s > %s", exe, "input.com", "input.log"); - } - fprintf(stderr, "Calling '%s'\n", buf); - if (system(buf) != 0) - { - gmx_fatal(FARGS, "Call to '%s' failed\n", buf); - } -} - -real call_gaussian(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]) -{ - /* normal gaussian jobs */ - static int step = 0; - int i, j; - real QMener = 0.0; - rvec * QMgrad, *MMgrad; - char* exe; - - snew(exe, 30); - sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe); - snew(QMgrad, qm->nrQMatoms); - snew(MMgrad, mm->nrMMatoms); - - write_gaussian_input(step, qmmm, qm, mm); - do_gaussian(step, exe); - QMener = read_gaussian_output(QMgrad, MMgrad, qm, mm); - /* put the QMMM forces in the force array and to the fshift - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i][j] = HARTREE_BOHR2MD * QMgrad[i][j]; - fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j]; - } - } - for (i = 0; i < mm->nrMMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j]; - fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j]; - } - } - QMener = QMener * HARTREE2KJ * AVOGADRO; - step++; - free(exe); - return (QMener); - -} /* call_gaussian */ - -real call_gaussian_SH(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]) -{ - /* a gaussian call routine intended for doing diabatic surface - * "sliding". See the manual for the theoretical background of this - * TSH method. - */ - static int step = 0; - int state, i, j; - real QMener = 0.0; - static gmx_bool swapped = FALSE; /* handle for identifying the current PES */ - gmx_bool swap = FALSE; /* the actual swap */ - rvec * QMgrad, *MMgrad; - char* buf; - char* exe; - - snew(exe, 30); - sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe); - /* hack to do ground state simulations */ - if (!step) - { - snew(buf, 20); - buf = getenv("GMX_QM_GROUND_STATE"); - if (buf) - { - sscanf(buf, "%d", &state); - } - else - { - state = 2; - } - if (state == 1) - { - swapped = TRUE; - } - } - /* end of hack */ - - - /* copy the QMMMrec pointer */ - snew(QMgrad, qm->nrQMatoms); - snew(MMgrad, mm->nrMMatoms); - /* at step 0 there should be no SA */ - /* if(!step) - * qr->bSA=FALSE;*/ - /* temporray set to step + 1, since there is a chk start */ - write_gaussian_SH_input(step, swapped, qmmm, qm, mm); - - do_gaussian(step, exe); - QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm); - - /* check for a surface hop. Only possible if we were already state - * averaging. - */ - if (qm->SAstep > 0) - { - if (!swapped) - { - swap = ((step != 0) && (hop(step, qm) != 0)); - swapped = swap; - } - else /* already on the other surface, so check if we go back */ - { - swap = ((step != 0) && (hop(step, qm) != 0)); - swapped = !swap; /* so swapped shoud be false again */ - } - if (swap) /* change surface, so do another call */ - { - write_gaussian_SH_input(step, swapped, qmmm, qm, mm); - do_gaussian(step, exe); - QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm); - } - } - /* add the QMMM forces to the gmx force array and fshift - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i][j] = HARTREE_BOHR2MD * QMgrad[i][j]; - fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j]; - } - } - for (i = 0; i < mm->nrMMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j]; - fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j]; - } - } - QMener = QMener * HARTREE2KJ * AVOGADRO; - fprintf(stderr, "step %5d, SA = %5d, swap = %5d\n", step, static_cast(qm->SAstep > 0), - static_cast(swapped)); - step++; - free(exe); - return (QMener); - -} /* call_gaussian_SH */ - -#pragma GCC diagnostic pop diff --git a/src/gromacs/mdlib/qm_gaussian.h b/src/gromacs/mdlib/qm_gaussian.h deleted file mode 100644 index 9c30884132..0000000000 --- a/src/gromacs/mdlib/qm_gaussian.h +++ /dev/null @@ -1,72 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2018,2019, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#ifndef GMX_MDLIB_QMGAUSSIAN_H -#define GMX_MDLIB_QMGAUSSIAN_H - -#include "gromacs/math/vectypes.h" -#include "gromacs/mdlib/qmmm.h" - -struct t_forcerec; - -/*! \brief - * Initialize gaussian datastructures. - * - * \param[in] qm QM forcerec. - */ -void init_gaussian(t_QMrec* qm); - -/*! \brief - * Call gaussian to do qm calculation. - * - * \param[in] qmmm QMMM part forcerec. - * \param[in] qm QM part of forcerec. - * \param[in] mm mm part of forcerec. - * \param[in] f force vector. - * \param[in] fshift shift of force vector. - */ -real call_gaussian(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]); - -/*! \brief - * Call gaussian SH(?) to do qm calculation. - * - * \param[in] qmmm QMMM part forcerec. - * \param[in] qm QM part of forcerec. - * \param[in] mm mm part of forcerec. - * \param[in] f force vector. - * \param[in] fshift shift of force vector. - */ -real call_gaussian_SH(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]); - -#endif diff --git a/src/gromacs/mdlib/qm_mopac.cpp b/src/gromacs/mdlib/qm_mopac.cpp deleted file mode 100644 index 5fd21f49db..0000000000 --- a/src/gromacs/mdlib/qm_mopac.cpp +++ /dev/null @@ -1,260 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gmxpre.h" - -#include "qm_mopac.h" - -#include "config.h" - -#include -#include -#include - -#include - -#include "gromacs/fileio/confio.h" -#include "gromacs/gmxlib/network.h" -#include "gromacs/gmxlib/nrnb.h" -#include "gromacs/math/units.h" -#include "gromacs/math/vec.h" -#include "gromacs/mdlib/qmmm.h" -#include "gromacs/mdtypes/md_enums.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/smalloc.h" - - -// When not built in a configuration with QMMM support, much of this -// code is unreachable by design. Tell clang not to warn about it. -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wmissing-noreturn" - -#if GMX_QMMM_MOPAC -/* mopac interface routines */ -void F77_FUNC(domldt, DOMLDT)(int* nrqmat, int labels[], char keywords[]); - -void F77_FUNC(domop, DOMOP)(int* nrqmat, - double qmcrd[], - int* nrmmat, - double mmchrg[], - double mmcrd[], - double qmgrad[], - double mmgrad[], - double* energy, - double qmcharges[]); - -#else /* GMX_QMMM_MOPAC */ -// Stub definitions to make compilation succeed when not configured -// for MOPAC support. In that case, the module gives a fatal error -// when the initialization function is called, so there is no need to -// issue fatal errors here, because that introduces problems with -// tools suggesting and prohibiting noreturn attributes. - -static void F77_FUNC(domldt, DOMLDT)(int* /*unused*/, int /*unused*/[], char /*unused*/[]) {} - -static void F77_FUNC(domop, DOMOP)(int* /*unused*/, - double /*unused*/[], - int* /*unused*/, - double /*unused*/[], - double /*unused*/[], - double /*unused*/[], - double /*unused*/[], - double* /*unused*/, - double /*unused*/[]) -{ -} - -#endif - - -void init_mopac(t_QMrec* qm) -{ - /* initializes the mopac routines ans sets up the semiempirical - * computation by calling moldat(). The inline mopac routines can - * only perform gradient operations. If one would like to optimize a - * structure or find a transition state at PM3 level, gaussian is - * used instead. - */ - char* keywords; - - if (!GMX_QMMM_MOPAC) - { - gmx_fatal(FARGS, - "Cannot call MOPAC unless linked against it. Use cmake -DGMX_QMMM_PROGRAM=MOPAC, " - "and ensure that linking will work correctly."); - } - - snew(keywords, 240); - - if (!qm->bSH) /* if rerun then grad should not be done! */ - { - sprintf(keywords, "PRECISE GEO-OK CHARGE=%d GRAD MMOK ANALYT %s\n", qm->QMcharge, - eQMmethod_names[qm->QMmethod]); - } - else - { - sprintf(keywords, "PRECISE GEO-OK CHARGE=%d SINGLET GRAD %s C.I.=(%d,%d) root=2 MECI \n", - qm->QMcharge, eQMmethod_names[qm->QMmethod], qm->CASorbitals, qm->CASelectrons / 2); - } - F77_FUNC(domldt, DOMLDT)(&qm->nrQMatoms, qm->atomicnumberQM, keywords); - fprintf(stderr, "keywords are: %s\n", keywords); - free(keywords); - -} /* init_mopac */ - -real call_mopac(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]) -{ - /* do the actual QMMM calculation using directly linked mopac subroutines - */ - double /* always double as the MOPAC routines are always compiled in - double precission! */ - *qmcrd = nullptr, - *qmchrg = nullptr, *mmcrd = nullptr, *mmchrg = nullptr, *qmgrad, *mmgrad = nullptr, - energy = 0; - int i, j; - real QMener = 0.0; - snew(qmcrd, 3 * (qm->nrQMatoms)); - snew(qmgrad, 3 * (qm->nrQMatoms)); - /* copy the data from qr into the arrays that are going to be used - * in the fortran routines of MOPAC - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - qmcrd[3 * i + j] = static_cast(qm->xQM[i][j]) * 10; - } - } - if (mm->nrMMatoms) - { - /* later we will add the point charges here. There are some - * conceptual problems with semi-empirical QM in combination with - * point charges that we need to solve first.... - */ - gmx_fatal(FARGS, - "At present only ONIOM is allowed in combination" - " with MOPAC QM subroutines\n"); - } - else - { - /* now compute the energy and the gradients. - */ - - snew(qmchrg, qm->nrQMatoms); - F77_FUNC(domop, DOMOP) - (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mmchrg, mmcrd, qmgrad, mmgrad, &energy, qmchrg); - /* add the gradients to the f[] array, and also to the fshift[]. - * the mopac gradients are in kCal/angstrom. - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i][j] = static_cast(10) * CAL2JOULE * qmgrad[3 * i + j]; - fshift[i][j] = static_cast(10) * CAL2JOULE * qmgrad[3 * i + j]; - } - } - QMener = static_cast CAL2JOULE * energy; - /* do we do something with the mulliken charges?? */ - - free(qmchrg); - } - free(qmgrad); - free(qmcrd); - return (QMener); -} - -real call_mopac_SH(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]) -{ - /* do the actual SH QMMM calculation using directly linked mopac - subroutines */ - - double /* always double as the MOPAC routines are always compiled in - double precission! */ - *qmcrd = nullptr, - *qmchrg = nullptr, *mmcrd = nullptr, *mmchrg = nullptr, *qmgrad, *mmgrad = nullptr, - energy = 0; - int i, j; - real QMener = 0.0; - - snew(qmcrd, 3 * (qm->nrQMatoms)); - snew(qmgrad, 3 * (qm->nrQMatoms)); - /* copy the data from qr into the arrays that are going to be used - * in the fortran routines of MOPAC - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - qmcrd[3 * i + j] = static_cast(qm->xQM[i][j]) * 10; - } - } - if (mm->nrMMatoms) - { - /* later we will add the point charges here. There are some - * conceptual problems with semi-empirical QM in combination with - * point charges that we need to solve first.... - */ - gmx_fatal(FARGS, "At present only ONIOM is allowed in combination with MOPAC\n"); - } - else - { - /* now compute the energy and the gradients. - */ - snew(qmchrg, qm->nrQMatoms); - - F77_FUNC(domop, DOMOP) - (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mmchrg, mmcrd, qmgrad, mmgrad, &energy, qmchrg); - /* add the gradients to the f[] array, and also to the fshift[]. - * the mopac gradients are in kCal/angstrom. - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i][j] = static_cast(10) * CAL2JOULE * qmgrad[3 * i + j]; - fshift[i][j] = static_cast(10) * CAL2JOULE * qmgrad[3 * i + j]; - } - } - QMener = static_cast CAL2JOULE * energy; - } - free(qmgrad); - free(qmcrd); - return (QMener); -} /* call_mopac_SH */ - -#pragma GCC diagnostic pop diff --git a/src/gromacs/mdlib/qm_mopac.h b/src/gromacs/mdlib/qm_mopac.h deleted file mode 100644 index 0f6eca8776..0000000000 --- a/src/gromacs/mdlib/qm_mopac.h +++ /dev/null @@ -1,67 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2018,2019, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#ifndef GMX_MDLIB_QMMOPAC_H -#define GMX_MDLIB_QMMOPAC_H - -#include "gromacs/mdlib/qmmm.h" - -/*! \brief - * Initialize mopac datastructures. - * - * \param[in] qm QM forcerec. - */ -void init_mopac(t_QMrec* qm); - -/*! \brief - * Run calculation with MOPAC. - * - * \param[in] qm QM part of forcerec. - * \param[in] mm MM part of forcerec. - * \param[in] f Force vector. - * \param[in] fshift Force shift vector. - */ -real call_mopac(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]); - -/*! \brief - * Run surface-hopping calculation with MOPAC. - * - * \param[in] qm QM part of forcerec. - * \param[in] mm MM part of forcerec. - * \param[in] f Force vector. - * \param[in] fshift Force shift vector. - */ -real call_mopac_SH(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]); - -#endif diff --git a/src/gromacs/mdlib/qm_orca.cpp b/src/gromacs/mdlib/qm_orca.cpp deleted file mode 100644 index 470b9a5b1b..0000000000 --- a/src/gromacs/mdlib/qm_orca.cpp +++ /dev/null @@ -1,369 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2008, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gmxpre.h" - -#include "qm_orca.h" - -#include "config.h" - -#include -#include -#include -#include - -#include "gromacs/fileio/confio.h" -#include "gromacs/gmxlib/network.h" -#include "gromacs/gmxlib/nrnb.h" -#include "gromacs/math/units.h" -#include "gromacs/math/vec.h" -#include "gromacs/mdlib/qmmm.h" -#include "gromacs/mdtypes/md_enums.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/smalloc.h" - -// When not built in a configuration with QMMM support, much of this -// code is unreachable by design. Tell clang not to warn about it. -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wmissing-noreturn" - -/* ORCA interface routines */ - -void init_orca(t_QMrec* qm) -{ - char* buf; - snew(buf, 200); - - if (!GMX_QMMM_ORCA) - { - gmx_fatal(FARGS, - "Cannot call ORCA unless linked against it. Use cmake -DGMX_QMMM_PROGRAM=ORCA, " - "and ensure that linking will work correctly."); - } - - /* ORCA settings on the system */ - buf = getenv("GMX_QM_ORCA_BASENAME"); - if (buf) - { - snew(qm->orca_basename, 200); - sscanf(buf, "%s", qm->orca_basename); - } - else - { - gmx_fatal(FARGS, "$GMX_QM_ORCA_BASENAME is not set\n"); - } - - /* ORCA directory on the system */ - snew(buf, 200); - buf = getenv("GMX_ORCA_PATH"); - - if (buf) - { - snew(qm->orca_dir, 200); - sscanf(buf, "%s", qm->orca_dir); - } - else - { - gmx_fatal(FARGS, "$GMX_ORCA_PATH not set, check manual\n"); - } - - fprintf(stderr, "Setting ORCA path to: %s...\n", qm->orca_dir); - fprintf(stderr, "ORCA initialised...\n\n"); - /* since we append the output to the BASENAME.out file, - we should delete an existent old out-file here. */ - sprintf(buf, "%s.out", qm->orca_basename); - remove(buf); -} - - -static void write_orca_input(const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm) -{ - int i; - FILE *pcFile, *addInputFile; - char *buf, *orcaInput, *addInputFilename, *pcFilename; - - /* write the first part of the input-file */ - snew(orcaInput, 200); - sprintf(orcaInput, "%s.inp", qm->orca_basename); - FILE* out = fopen(orcaInput, "w"); - - snew(addInputFilename, 200); - sprintf(addInputFilename, "%s.ORCAINFO", qm->orca_basename); - addInputFile = fopen(addInputFilename, "r"); - - fprintf(out, "#input-file generated by GROMACS\n"); - - fprintf(out, "!EnGrad TightSCF\n"); - - /* here we include the insertion of the additional orca-input */ - snew(buf, 200); - if (addInputFile != nullptr) - { - while (!feof(addInputFile)) - { - if (fgets(buf, 200, addInputFile) != nullptr) - { - fputs(buf, out); - } - } - } - else - { - gmx_fatal(FARGS, "No information on the calculation given in %s\n", addInputFilename); - } - - fclose(addInputFile); - - /* write charge and multiplicity */ - fprintf(out, "*xyz %2d%2d\n", qm->QMcharge, qm->multiplicity); - - /* write the QM coordinates */ - for (i = 0; i < qm->nrQMatoms; i++) - { - int atomNr; - if (qm->atomicnumberQM[i] == 0) - { - atomNr = 1; - } - else - { - atomNr = qm->atomicnumberQM[i]; - } - fprintf(out, "%3d %10.7f %10.7f %10.7f\n", atomNr, qm->xQM[i][XX] / 0.1, - qm->xQM[i][YY] / 0.1, qm->xQM[i][ZZ] / 0.1); - } - fprintf(out, "*\n"); - - /* write the MM point charge data */ - if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) - { - /* name of the point charge file */ - snew(pcFilename, 200); - sprintf(pcFilename, "%s.pc", qm->orca_basename); - fprintf(out, "%s%s%s\n", "%pointcharges \"", pcFilename, "\""); - pcFile = fopen(pcFilename, "w"); - fprintf(pcFile, "%d\n", mm->nrMMatoms); - for (i = 0; i < mm->nrMMatoms; i++) - { - fprintf(pcFile, "%8.4f %10.7f %10.7f %10.7f\n", mm->MMcharges[i], - mm->xMM[i][XX] / 0.1, mm->xMM[i][YY] / 0.1, mm->xMM[i][ZZ] / 0.1); - } - fprintf(pcFile, "\n"); - fclose(pcFile); - } - fprintf(out, "\n"); - - fclose(out); -} /* write_orca_input */ - -static real read_orca_output(rvec QMgrad[], rvec MMgrad[], const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm) -{ - int i, j; - char buf[300], orca_pcgradFilename[300], orca_engradFilename[300]; - real QMener; - FILE *pcgrad, *engrad; - int k; - - /* the energy and gradients for the QM part are stored in the engrad file - * and the gradients for the point charges are stored in the pc file. - */ - sprintf(orca_engradFilename, "%s.engrad", qm->orca_basename); - engrad = fopen(orca_engradFilename, "r"); - /* we read the energy and the gradient for the qm-atoms from the engrad file - */ - /* we can skip the first seven lines - */ - for (j = 0; j < 7; j++) - { - if (fgets(buf, 300, engrad) == nullptr) - { - gmx_fatal(FARGS, "Unexpected end of ORCA output"); - } - } - /* now comes the energy - */ - if (fgets(buf, 300, engrad) == nullptr) - { - gmx_fatal(FARGS, "Unexpected end of ORCA output"); - } -#if GMX_DOUBLE - sscanf(buf, "%lf\n", &QMener); -#else - sscanf(buf, "%f\n", &QMener); -#endif - /* we can skip the next three lines - */ - for (j = 0; j < 3; j++) - { - if (fgets(buf, 300, engrad) == nullptr) - { - gmx_fatal(FARGS, "Unexpected end of ORCA output"); - } - } - /* next lines contain the gradients of the QM atoms - * now comes the gradient, one value per line: - * (atom1 x \n atom1 y \n atom1 z \n atom2 x ... - */ - - for (i = 0; i < 3 * qm->nrQMatoms; i++) - { - k = i / 3; - if (fgets(buf, 300, engrad) == nullptr) - { - gmx_fatal(FARGS, "Unexpected end of ORCA output"); - } -#if GMX_DOUBLE - if (i % 3 == 0) - { - sscanf(buf, "%lf\n", &QMgrad[k][XX]); - } - else if (i % 3 == 1) - { - sscanf(buf, "%lf\n", &QMgrad[k][YY]); - } - else if (i % 3 == 2) - { - sscanf(buf, "%lf\n", &QMgrad[k][ZZ]); - } -#else - if (i % 3 == 0) - { - sscanf(buf, "%f\n", &QMgrad[k][XX]); - } - else if (i % 3 == 1) - { - sscanf(buf, "%f\n", &QMgrad[k][YY]); - } - else if (i % 3 == 2) - { - sscanf(buf, "%f\n", &QMgrad[k][ZZ]); - } -#endif - } - fclose(engrad); - /* write the MM point charge data - */ - if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) - { - sprintf(orca_pcgradFilename, "%s.pcgrad", qm->orca_basename); - pcgrad = fopen(orca_pcgradFilename, "r"); - - /* we read the gradient for the mm-atoms from the pcgrad file - */ - /* we can skip the first line - */ - if (fgets(buf, 300, pcgrad) == nullptr) - { - gmx_fatal(FARGS, "Unexpected end of ORCA output"); - } - for (i = 0; i < mm->nrMMatoms; i++) - { - if (fgets(buf, 300, pcgrad) == nullptr) - { - gmx_fatal(FARGS, "Unexpected end of ORCA output"); - } -#if GMX_DOUBLE - sscanf(buf, "%lf%lf%lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]); -#else - sscanf(buf, "%f%f%f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]); -#endif - } - fclose(pcgrad); - } - return (QMener); -} - -static void do_orca(char* orca_dir, char* basename) -{ - - /* make the call to the orca binary through system() - * The location of the binary is set through the - * environment. - */ - char buf[100]; - sprintf(buf, "%s/%s %s.inp >> %s.out", orca_dir, "orca", basename, basename); - fprintf(stderr, "Calling '%s'\n", buf); - if (system(buf) != 0) - { - gmx_fatal(FARGS, "Call to '%s' failed\n", buf); - } -} - -real call_orca(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]) -{ - /* normal orca jobs */ - static int step = 0; - int i, j; - real QMener; - rvec * QMgrad, *MMgrad; - char* exe; - - snew(exe, 30); - sprintf(exe, "%s", "orca"); - snew(QMgrad, qm->nrQMatoms); - snew(MMgrad, mm->nrMMatoms); - - write_orca_input(qmmm, qm, mm); - do_orca(qm->orca_dir, qm->orca_basename); - QMener = read_orca_output(QMgrad, MMgrad, qmmm, qm, mm); - /* put the QMMM forces in the force array and to the fshift - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i][j] = HARTREE_BOHR2MD * QMgrad[i][j]; - fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j]; - } - } - for (i = 0; i < mm->nrMMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - f[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j]; - fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j]; - } - } - QMener = QMener * HARTREE2KJ * AVOGADRO; - step++; - free(exe); - return (QMener); -} /* call_orca */ - -/* end of orca sub routines */ - -#pragma GCC diagnostic pop diff --git a/src/gromacs/mdlib/qm_orca.h b/src/gromacs/mdlib/qm_orca.h deleted file mode 100644 index c6a3b55917..0000000000 --- a/src/gromacs/mdlib/qm_orca.h +++ /dev/null @@ -1,44 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2017,2018,2019, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gromacs/mdlib/qmmm.h" - -#ifndef GMX_MDLIB_QM_ORCA_H -# define GMX_MDLIB_QM_ORCA_H - -void init_orca(t_QMrec* qm); - -real call_orca(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]); - -#endif diff --git a/src/gromacs/mdlib/qmmm.cpp b/src/gromacs/mdlib/qmmm.cpp deleted file mode 100644 index 508f87cef3..0000000000 --- a/src/gromacs/mdlib/qmmm.cpp +++ /dev/null @@ -1,889 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gmxpre.h" - -#include "qmmm.h" - -#include "config.h" - -#include -#include -#include -#include - -#include - -#include "gromacs/domdec/domdec_struct.h" -#include "gromacs/fileio/confio.h" -#include "gromacs/gmxlib/network.h" -#include "gromacs/gmxlib/nrnb.h" -#include "gromacs/math/functions.h" -#include "gromacs/math/units.h" -#include "gromacs/math/vec.h" -#include "gromacs/mdlib/qm_gamess.h" -#include "gromacs/mdlib/qm_gaussian.h" -#include "gromacs/mdlib/qm_mopac.h" -#include "gromacs/mdlib/qm_orca.h" -#include "gromacs/mdtypes/commrec.h" -#include "gromacs/mdtypes/forceoutput.h" -#include "gromacs/mdtypes/forcerec.h" -#include "gromacs/mdtypes/inputrec.h" -#include "gromacs/mdtypes/md_enums.h" -#include "gromacs/mdtypes/mdatom.h" -#include "gromacs/mdtypes/nblist.h" -#include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/pbc.h" -#include "gromacs/topology/mtop_lookup.h" -#include "gromacs/topology/mtop_util.h" -#include "gromacs/topology/topology.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/smalloc.h" - -// When not built in a configuration with QMMM support, much of this -// code is unreachable by design. Tell clang not to warn about it. -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wunreachable-code" -#pragma GCC diagnostic ignored "-Wmissing-noreturn" - -/* this struct and these comparison functions are needed for creating - * a QMMM input for the QM routines from the QMMM neighbor list. - */ - -typedef struct -{ - int j; - int shift; -} t_j_particle; - -static bool struct_comp(const t_j_particle& a, const t_j_particle& b) -{ - return a.j < b.j; -} - -static real call_QMroutine(const t_commrec gmx_unused* cr, - const t_QMMMrec gmx_unused* qmmm, - t_QMrec gmx_unused* qm, - t_MMrec gmx_unused* mm, - rvec gmx_unused f[], - rvec gmx_unused fshift[]) -{ - /* makes a call to the requested QM routine (qm->QMmethod) - * Note that f is actually the gradient, i.e. -f - */ - /* do a semi-empiprical calculation */ - - if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms)) - { - if (GMX_QMMM_MOPAC) - { - if (qm->bSH) - { - return call_mopac_SH(qm, mm, f, fshift); - } - else - { - return call_mopac(qm, mm, f, fshift); - } - } - else - { - gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac."); - } - } - else - { - /* do an ab-initio calculation */ - if (qm->bSH && qm->QMmethod == eQMmethodCASSCF) - { - if (GMX_QMMM_GAUSSIAN) - { - return call_gaussian_SH(qmmm, qm, mm, f, fshift); - } - else - { - gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian."); - } - } - else - { - if (GMX_QMMM_GAMESS) - { - return call_gamess(qm, mm, f, fshift); - } - else if (GMX_QMMM_GAUSSIAN) - { - return call_gaussian(qmmm, qm, mm, f, fshift); - } - else if (GMX_QMMM_ORCA) - { - return call_orca(qmmm, qm, mm, f, fshift); - } - else - { - gmx_fatal(FARGS, - "Ab-initio calculation only supported with Gamess, Gaussian or ORCA."); - } - } - } -} - -static void init_QMroutine(const t_commrec gmx_unused* cr, t_QMrec gmx_unused* qm, t_MMrec gmx_unused* mm) -{ - /* makes a call to the requested QM routine (qm->QMmethod) - */ - if (qm->QMmethod < eQMmethodRHF) - { - if (GMX_QMMM_MOPAC) - { - /* do a semi-empiprical calculation */ - init_mopac(qm); - } - else - { - gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac."); - } - } - else - { - /* do an ab-initio calculation */ - if (GMX_QMMM_GAMESS) - { - init_gamess(cr, qm, mm); - } - else if (GMX_QMMM_GAUSSIAN) - { - init_gaussian(qm); - } - else if (GMX_QMMM_ORCA) - { - init_orca(qm); - } - else - { - gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA."); - } - } -} /* init_QMroutine */ - -static void update_QMMM_coord(const rvec* x, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm) -{ - /* shifts the QM and MM particles into the central box and stores - * these shifted coordinates in the coordinate arrays of the - * QMMMrec. These coordinates are passed on the QM subroutines. - */ - int i; - - /* shift the QM atoms into the central box - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]); - } - /* also shift the MM atoms into the central box, if any - */ - for (i = 0; i < mm->nrMMatoms; i++) - { - rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]); - } -} /* update_QMMM_coord */ - -/* end of QMMM subroutines */ - -/* QMMM core routines */ - -static t_QMrec* mk_QMrec() -{ - t_QMrec* qm; - snew(qm, 1); - return qm; -} /* mk_QMrec */ - -static t_MMrec* mk_MMrec() -{ - t_MMrec* mm; - snew(mm, 1); - return mm; -} /* mk_MMrec */ - -static void init_QMrec(int grpnr, t_QMrec* qm, int nr, const int* atomarray, const gmx_mtop_t* mtop, const t_inputrec* ir) -{ - /* fills the t_QMrec struct of QM group grpnr - */ - - qm->nrQMatoms = nr; - snew(qm->xQM, nr); - snew(qm->indexQM, nr); - snew(qm->shiftQM, nr); /* the shifts */ - for (int i = 0; i < nr; i++) - { - qm->indexQM[i] = atomarray[i]; - } - - snew(qm->atomicnumberQM, nr); - int molb = 0; - for (int i = 0; i < qm->nrQMatoms; i++) - { - const t_atom& atom = mtopGetAtomParameters(mtop, qm->indexQM[i], &molb); - qm->nelectrons += mtop->atomtypes.atomnumber[atom.type]; - qm->atomicnumberQM[i] = mtop->atomtypes.atomnumber[atom.type]; - } - - qm->QMcharge = ir->opts.QMcharge[grpnr]; - qm->multiplicity = ir->opts.QMmult[grpnr]; - qm->nelectrons -= ir->opts.QMcharge[grpnr]; - - qm->QMmethod = ir->opts.QMmethod[grpnr]; - qm->QMbasis = ir->opts.QMbasis[grpnr]; - /* trajectory surface hopping setup (Gaussian only) */ - qm->bSH = ir->opts.bSH[grpnr]; - qm->CASorbitals = ir->opts.CASorbitals[grpnr]; - qm->CASelectrons = ir->opts.CASelectrons[grpnr]; - qm->SAsteps = ir->opts.SAsteps[grpnr]; - qm->SAon = ir->opts.SAon[grpnr]; - qm->SAoff = ir->opts.SAoff[grpnr]; - /* hack to prevent gaussian from reinitializing all the time */ - qm->nQMcpus = 0; /* number of CPU's to be used by g01, is set - * upon initializing gaussian - * (init_gaussian() - */ - /* print the current layer to allow users to check their input */ - fprintf(stderr, "Layer %d\nnr of QM atoms %d\n", grpnr, nr); - fprintf(stderr, "QMlevel: %s/%s\n\n", eQMmethod_names[qm->QMmethod], eQMbasis_names[qm->QMbasis]); -} /* init_QMrec */ - -static t_QMrec* copy_QMrec(t_QMrec* qm) -{ - /* copies the contents of qm into a new t_QMrec struct */ - t_QMrec* qmcopy; - int i; - - qmcopy = mk_QMrec(); - qmcopy->nrQMatoms = qm->nrQMatoms; - snew(qmcopy->xQM, qmcopy->nrQMatoms); - snew(qmcopy->indexQM, qmcopy->nrQMatoms); - snew(qmcopy->atomicnumberQM, qm->nrQMatoms); - snew(qmcopy->shiftQM, qmcopy->nrQMatoms); /* the shifts */ - for (i = 0; i < qmcopy->nrQMatoms; i++) - { - qmcopy->shiftQM[i] = qm->shiftQM[i]; - qmcopy->indexQM[i] = qm->indexQM[i]; - qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i]; - } - qmcopy->nelectrons = qm->nelectrons; - qmcopy->multiplicity = qm->multiplicity; - qmcopy->QMcharge = qm->QMcharge; - qmcopy->nelectrons = qm->nelectrons; - qmcopy->QMmethod = qm->QMmethod; - qmcopy->QMbasis = qm->QMbasis; - /* trajectory surface hopping setup (Gaussian only) */ - qmcopy->bSH = qm->bSH; - qmcopy->CASorbitals = qm->CASorbitals; - qmcopy->CASelectrons = qm->CASelectrons; - qmcopy->SAsteps = qm->SAsteps; - qmcopy->SAon = qm->SAon; - qmcopy->SAoff = qm->SAoff; - - /* Gaussian init. variables */ - qmcopy->nQMcpus = qm->nQMcpus; - for (i = 0; i < DIM; i++) - { - qmcopy->SHbasis[i] = qm->SHbasis[i]; - } - qmcopy->QMmem = qm->QMmem; - qmcopy->accuracy = qm->accuracy; - qmcopy->cpmcscf = qm->cpmcscf; - qmcopy->SAstep = qm->SAstep; - - return (qmcopy); - -} /*copy_QMrec */ - -#if GMX_QMMM - -t_QMMMrec* mk_QMMMrec() -{ - t_QMMMrec* qr; - - snew(qr, 1); - - return qr; - -} /* mk_QMMMrec */ - -#else /* GMX_QMMM */ - -t_QMMMrec* mk_QMMMrec() -{ - gmx_incons("Compiled without QMMM"); -} /* mk_QMMMrec */ -#endif - -void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr) -{ - /* we put the atomsnumbers of atoms that belong to the QMMM group in - * an array that will be copied later to QMMMrec->indexQM[..]. Also - * it will be used to create an QMMMrec->bQMMM index array that - * simply contains true/false for QM and MM (the other) atoms. - */ - - t_QMMMrec* qr; - t_MMrec* mm; - - if (!GMX_QMMM) - { - gmx_incons("Compiled without QMMM"); - } - - if (ir->cutoff_scheme != ecutsGROUP) - { - gmx_fatal(FARGS, "QMMM is currently only supported with cutoff-scheme=group"); - } - if (!EI_DYNAMICS(ir->eI)) - { - gmx_fatal(FARGS, "QMMM is only supported with dynamics"); - } - - /* issue a fatal if the user wants to run with more than one node */ - if (PAR(cr)) - { - gmx_fatal(FARGS, "QM/MM does not work in parallel, use a single rank instead\n"); - } - - /* Make a local copy of the QMMMrec */ - qr = fr->qr; - - /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are - * QM/not QM. We first set all elemenst at false. Afterwards we use - * the qm_arr (=MMrec->indexQM) to changes the elements - * corresponding to the QM atoms at TRUE. */ - - qr->QMMMscheme = ir->QMMMscheme; - - /* we take the possibility into account that a user has - * defined more than one QM group: - */ - /* an ugly work-around in case there is only one group In this case - * the whole system is treated as QM. Otherwise the second group is - * always the rest of the total system and is treated as MM. - */ - - /* small problem if there is only QM.... so no MM */ - - int numQmmmGroups = ir->opts.ngQM; - - if (qr->QMMMscheme == eQMMMschemeoniom) - { - qr->nrQMlayers = numQmmmGroups; - } - else - { - qr->nrQMlayers = 1; - } - - /* there are numQmmmGroups groups of QM atoms. In case of multiple QM groups - * I assume that the users wants to do ONIOM. However, maybe it - * should also be possible to define more than one QM subsystem with - * independent neighbourlists. I have to think about - * that.. 11-11-2003 - */ - std::vector qmmmAtoms; - snew(qr->qm, numQmmmGroups); - for (int i = 0; i < numQmmmGroups; i++) - { - /* new layer */ - if (qr->QMMMscheme == eQMMMschemeoniom) - { - /* add the atoms to the bQMMM array - */ - - /* I assume that users specify the QM groups from small to - * big(ger) in the mdp file - */ - qr->qm[i] = mk_QMrec(); - /* store QM atoms in this layer in the QMrec and initialise layer - */ - init_QMrec(i, qr->qm[i], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir); - } - } - if (qr->QMMMscheme != eQMMMschemeoniom) - { - - /* standard QMMM, all layers are merged together so there is one QM - * subsystem and one MM subsystem. - * Also we set the charges to zero in mtop to prevent the innerloops - * from doubly counting the electostatic QM MM interaction - * TODO: Consider doing this in grompp instead. - */ - - qr->qm[0] = mk_QMrec(); - /* store QM atoms in the QMrec and initialise - */ - init_QMrec(0, qr->qm[0], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir); - - /* MM rec creation */ - mm = mk_MMrec(); - mm->scalefactor = ir->scalefactor; - mm->nrMMatoms = (mtop->natoms) - (qr->qm[0]->nrQMatoms); /* rest of the atoms */ - qr->mm = mm; - } - else /* ONIOM */ - { /* MM rec creation */ - mm = mk_MMrec(); - mm->scalefactor = ir->scalefactor; - mm->nrMMatoms = 0; - qr->mm = mm; - } - - /* these variables get updated in the update QMMMrec */ - - if (qr->nrQMlayers == 1) - { - /* with only one layer there is only one initialisation - * needed. Multilayer is a bit more complicated as it requires - * re-initialisation at every step of the simulation. This is due - * to the use of COMMON blocks in the fortran QM subroutines. - */ - if (qr->qm[0]->QMmethod < eQMmethodRHF) - { - if (GMX_QMMM_MOPAC) - { - /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */ - init_mopac(qr->qm[0]); - } - else - { - gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac."); - } - } - else - { - /* ab initio calculation requested (gamess/gaussian/ORCA) */ - if (GMX_QMMM_GAMESS) - { - init_gamess(cr, qr->qm[0], qr->mm); - } - else if (GMX_QMMM_GAUSSIAN) - { - init_gaussian(qr->qm[0]); - } - else if (GMX_QMMM_ORCA) - { - init_orca(qr->qm[0]); - } - else - { - gmx_fatal(FARGS, - "Ab-initio calculation only supported with Gamess, Gaussian or ORCA."); - } - } - } -} /* init_QMMMrec */ - -void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, const t_mdatoms* md, const matrix box) -{ - /* updates the coordinates of both QM atoms and MM atoms and stores - * them in the QMMMrec. - * - * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple - * ns needs to be fixed! - */ - int mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift; - t_j_particle *mm_j_particles = nullptr, *qm_i_particles = nullptr; - t_QMMMrec* qr; - t_nblist* QMMMlist; - rvec dx; - ivec crd; - t_QMrec* qm; - t_MMrec* mm; - t_pbc pbc; - int* parallelMMarray = nullptr; - - if (!GMX_QMMM) - { - gmx_incons("Compiled without QMMM"); - } - - /* every cpu has this array. On every processor we fill this array - * with 1's and 0's. 1's indicate the atoms is a QM atom on the - * current cpu in a later stage these arrays are all summed. indexes - * > 0 indicate the atom is a QM atom. Every node therefore knows - * whcih atoms are part of the QM subsystem. - */ - /* copy some pointers */ - qr = fr->qr; - mm = qr->mm; - QMMMlist = fr->QMMMlist; - - /* init_pbc(box); needs to be called first, see pbc.h */ - ivec null_ivec; - clear_ivec(null_ivec); - set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : null_ivec, FALSE, box); - /* only in standard (normal) QMMM we need the neighbouring MM - * particles to provide a electric field of point charges for the QM - * atoms. - */ - if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */ - { - /* we NOW create/update a number of QMMMrec entries: - * - * 1) the shiftQM, containing the shifts of the QM atoms - * - * 2) the indexMM array, containing the index of the MM atoms - * - * 3) the shiftMM, containing the shifts of the MM atoms - * - * 4) the shifted coordinates of the MM atoms - * - * the shifts are used for computing virial of the QM/MM particles. - */ - qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */ - snew(qm_i_particles, QMMMlist->nri); - if (QMMMlist->nri) - { - qm_i_particles[0].shift = XYZ2IS(0, 0, 0); - for (i = 0; i < QMMMlist->nri; i++) - { - qm_i_particles[i].j = QMMMlist->iinr[i]; - - if (i) - { - qm_i_particles[i].shift = - pbc_dx_aiuc(&pbc, x[QMMMlist->iinr[0]], x[QMMMlist->iinr[i]], dx); - } - /* However, since nri >= nrQMatoms, we do a quicksort, and throw - * out double, triple, etc. entries later, as we do for the MM - * list too. - */ - - /* compute the shift for the MM j-particles with respect to - * the QM i-particle and store them. - */ - - crd[0] = IS2X(QMMMlist->shift[i]) + IS2X(qm_i_particles[i].shift); - crd[1] = IS2Y(QMMMlist->shift[i]) + IS2Y(qm_i_particles[i].shift); - crd[2] = IS2Z(QMMMlist->shift[i]) + IS2Z(qm_i_particles[i].shift); - is = XYZ2IS(crd[0], crd[1], crd[2]); - for (j = QMMMlist->jindex[i]; j < QMMMlist->jindex[i + 1]; j++) - { - if (mm_nr >= mm_max) - { - mm_max += 1000; - srenew(mm_j_particles, mm_max); - } - - mm_j_particles[mm_nr].j = QMMMlist->jjnr[j]; - mm_j_particles[mm_nr].shift = is; - mm_nr++; - } - } - - /* quicksort QM and MM shift arrays and throw away multiple entries */ - - - std::sort(qm_i_particles, qm_i_particles + QMMMlist->nri, struct_comp); - /* The mm_j_particles argument to qsort is not allowed to be nullptr */ - if (mm_nr > 0) - { - std::sort(mm_j_particles, mm_j_particles + mm_nr, struct_comp); - } - /* remove multiples in the QM shift array, since in init_QMMM() we - * went through the atom numbers from 0 to md.nr, the order sorted - * here matches the one of QMindex already. - */ - j = 0; - for (i = 0; i < QMMMlist->nri; i++) - { - if (i == 0 || qm_i_particles[i].j != qm_i_particles[i - 1].j) - { - qm_i_particles[j++] = qm_i_particles[i]; - } - } - mm_nr_new = 0; - /* Remove double entries for the MM array. - * Also remove mm atoms that have no charges! - * actually this is already done in the ns.c - */ - for (i = 0; i < mm_nr; i++) - { - if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i - 1].j) - && !md->bQM[mm_j_particles[i].j] - && ((md->chargeA[mm_j_particles[i].j] != 0.0_real) - || (md->chargeB && (md->chargeB[mm_j_particles[i].j] != 0.0_real)))) - { - mm_j_particles[mm_nr_new++] = mm_j_particles[i]; - } - } - mm_nr = mm_nr_new; - /* store the data retrieved above into the QMMMrec - */ - k = 0; - /* Keep the compiler happy, - * shift will always be set in the loop for i=0 - */ - shift = 0; - for (i = 0; i < qm->nrQMatoms; i++) - { - /* not all qm particles might have appeared as i - * particles. They might have been part of the same charge - * group for instance. - */ - if (qm->indexQM[i] == qm_i_particles[k].j) - { - shift = qm_i_particles[k++].shift; - } - /* use previous shift, assuming they belong the same charge - * group anyway, - */ - - qm->shiftQM[i] = shift; - } - } - /* parallel excecution */ - if (PAR(cr)) - { - snew(parallelMMarray, 2 * (md->nr)); - /* only MM particles have a 1 at their atomnumber. The second part - * of the array contains the shifts. Thus: - * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM - * step or not. p[i+md->nr] is the shift of atomnumber i. - */ - for (i = 0; i < 2 * (md->nr); i++) - { - parallelMMarray[i] = 0; - } - - for (i = 0; i < mm_nr; i++) - { - parallelMMarray[mm_j_particles[i].j] = 1; - parallelMMarray[mm_j_particles[i].j + (md->nr)] = mm_j_particles[i].shift; - } - gmx_sumi(md->nr, parallelMMarray, cr); - mm_nr = 0; - - mm_max = 0; - for (i = 0; i < md->nr; i++) - { - if (parallelMMarray[i]) - { - if (mm_nr >= mm_max) - { - mm_max += 1000; - srenew(mm->indexMM, mm_max); - srenew(mm->shiftMM, mm_max); - } - mm->indexMM[mm_nr] = i; - mm->shiftMM[mm_nr++] = parallelMMarray[i + md->nr] / parallelMMarray[i]; - } - } - mm->nrMMatoms = mm_nr; - free(parallelMMarray); - } - /* serial execution */ - else - { - mm->nrMMatoms = mm_nr; - srenew(mm->shiftMM, mm_nr); - srenew(mm->indexMM, mm_nr); - for (i = 0; i < mm_nr; i++) - { - mm->indexMM[i] = mm_j_particles[i].j; - mm->shiftMM[i] = mm_j_particles[i].shift; - } - } - /* (re) allocate memory for the MM coordiate array. The QM - * coordinate array was already allocated in init_QMMM, and is - * only (re)filled in the update_QMMM_coordinates routine - */ - srenew(mm->xMM, mm->nrMMatoms); - /* now we (re) fill the array that contains the MM charges with - * the forcefield charges. If requested, these charges will be - * scaled by a factor - */ - srenew(mm->MMcharges, mm->nrMMatoms); - for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */ - { - mm->MMcharges[i] = md->chargeA[mm->indexMM[i]] * mm->scalefactor; - } - /* the next routine fills the coordinate fields in the QMMM rec of - * both the qunatum atoms and the MM atoms, using the shifts - * calculated above. - */ - - update_QMMM_coord(x, fr, qr->qm[0], qr->mm); - free(qm_i_particles); - free(mm_j_particles); - } - else /* ONIOM */ /* ????? */ - { - mm->nrMMatoms = 0; - /* do for each layer */ - for (j = 0; j < qr->nrQMlayers; j++) - { - qm = qr->qm[j]; - qm->shiftQM[0] = XYZ2IS(0, 0, 0); - for (i = 1; i < qm->nrQMatoms; i++) - { - qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]], dx); - } - update_QMMM_coord(x, fr, qm, mm); - } - } -} /* update_QMMM_rec */ - -real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qr) -{ - real QMener = 0.0; - /* a selection for the QM package depending on which is requested - * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now - * it works through defines.... Not so nice yet - */ - t_QMrec *qm, *qm2; - t_MMrec* mm = nullptr; - rvec * forces = nullptr, *fshift = nullptr, *forces2 = nullptr, - *fshift2 = nullptr; /* needed for multilayer ONIOM */ - int i, j, k; - - if (!GMX_QMMM) - { - gmx_incons("Compiled without QMMM"); - } - - /* make a local copy the QMMMrec pointer - */ - mm = qr->mm; - - /* now different procedures are carried out for one layer ONION and - * normal QMMM on one hand and multilayer oniom on the other - */ - gmx::ArrayRef fMM = forceWithShiftForces->force(); - gmx::ArrayRef fshiftMM = forceWithShiftForces->shiftForces(); - if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1) - { - qm = qr->qm[0]; - snew(forces, (qm->nrQMatoms + mm->nrMMatoms)); - snew(fshift, (qm->nrQMatoms + mm->nrMMatoms)); - QMener = call_QMroutine(cr, qr, qm, mm, forces, fshift); - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - fMM[qm->indexQM[i]][j] -= forces[i][j]; - fshiftMM[qm->shiftQM[i]][j] += fshift[i][j]; - } - } - for (i = 0; i < mm->nrMMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - fMM[mm->indexMM[i]][j] -= forces[qm->nrQMatoms + i][j]; - fshiftMM[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms + i][j]; - } - } - free(forces); - free(fshift); - } - else /* Multi-layer ONIOM */ - { - for (i = 0; i < qr->nrQMlayers - 1; i++) /* last layer is special */ - { - qm = qr->qm[i]; - qm2 = copy_QMrec(qr->qm[i + 1]); - - qm2->nrQMatoms = qm->nrQMatoms; - - for (j = 0; j < qm2->nrQMatoms; j++) - { - for (k = 0; k < DIM; k++) - { - qm2->xQM[j][k] = qm->xQM[j][k]; - } - qm2->indexQM[j] = qm->indexQM[j]; - qm2->atomicnumberQM[j] = qm->atomicnumberQM[j]; - qm2->shiftQM[j] = qm->shiftQM[j]; - } - - qm2->QMcharge = qm->QMcharge; - /* this layer at the higher level of theory */ - srenew(forces, qm->nrQMatoms); - srenew(fshift, qm->nrQMatoms); - /* we need to re-initialize the QMroutine every step... */ - init_QMroutine(cr, qm, mm); - QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift); - - /* this layer at the lower level of theory */ - srenew(forces2, qm->nrQMatoms); - srenew(fshift2, qm->nrQMatoms); - init_QMroutine(cr, qm2, mm); - QMener -= call_QMroutine(cr, qr, qm2, mm, forces2, fshift2); - /* E = E1high-E1low The next layer includes the current layer at - * the lower level of theory, which provides + E2low - * this is similar for gradients - */ - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - fMM[qm->indexQM[i]][j] -= (forces[i][j] - forces2[i][j]); - fshiftMM[qm->shiftQM[i]][j] += (fshift[i][j] - fshift2[i][j]); - } - } - free(qm2); - } - /* now the last layer still needs to be done: */ - qm = qr->qm[qr->nrQMlayers - 1]; /* C counts from 0 */ - init_QMroutine(cr, qm, mm); - srenew(forces, qm->nrQMatoms); - srenew(fshift, qm->nrQMatoms); - QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift); - for (i = 0; i < qm->nrQMatoms; i++) - { - for (j = 0; j < DIM; j++) - { - fMM[qm->indexQM[i]][j] -= forces[i][j]; - fshiftMM[qm->shiftQM[i]][j] += fshift[i][j]; - } - } - free(forces); - free(fshift); - free(forces2); - free(fshift2); - } - return (QMener); -} /* calculate_QMMM */ - -#pragma GCC diagnostic pop diff --git a/src/gromacs/mdlib/qmmm.h b/src/gromacs/mdlib/qmmm.h deleted file mode 100644 index 0eda21d216..0000000000 --- a/src/gromacs/mdlib/qmmm.h +++ /dev/null @@ -1,152 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2008, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS 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 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS 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 GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#ifndef GMX_MDLIB_QMMM_H -#define GMX_MDLIB_QMMM_H - -#include "config.h" - -#include - -#include "gromacs/math/vectypes.h" -#include "gromacs/mdlib/tgroup.h" -#include "gromacs/utility/arrayref.h" - -#define GMX_QMMM (GMX_QMMM_MOPAC || GMX_QMMM_GAMESS || GMX_QMMM_GAUSSIAN || GMX_QMMM_ORCA) - -struct gmx_localtop_t; -struct gmx_mtop_t; -struct t_commrec; -struct t_forcerec; -struct t_inputrec; -struct t_mdatoms; -struct t_QMMMrec; - -namespace gmx -{ -class ForceWithShiftForces; -} - -typedef struct -{ - int nrQMatoms; /* total nr of QM atoms */ - rvec* xQM; /* shifted to center of box */ - int* indexQM; /* atom i = atom indexQM[i] in mdrun */ - int* atomicnumberQM; /* atomic numbers of QM atoms */ - real* QMcharges; /* atomic charges of QM atoms(ONIOM) */ - int* shiftQM; - int QMcharge; /* charge of the QM system */ - int multiplicity; /* multipicity (no of unpaired eln) */ - int QMmethod; /* see enums.h for all methods */ - int QMbasis; /* see enums.h for all bases */ - int nelectrons; /* total number of elecs in QM region*/ - /* Gaussian specific stuff */ - int nQMcpus; /* no. of CPUs used for the QM calc. */ - int QMmem; /* memory for the gaussian calc. */ - int accuracy; /* convergence criterium (E(-x)) */ - gmx_bool cpmcscf; /* using cpmcscf(l1003)*/ - char* gauss_dir; - char* gauss_exe; - char* devel_dir; - char* orca_basename; /* basename for I/O with orca */ - char* orca_dir; /* directory for ORCA */ - /* Surface hopping stuff */ - gmx_bool bSH; /* surface hopping (diabatic only) */ - real SAon; /* at which energy gap the SA starts */ - real SAoff; /* at which energy gap the SA stops */ - int SAsteps; /* stepwise switchinng on the SA */ - int SAstep; /* current state of SA */ - int CIdim; - real* CIvec1; - real* CIvec2; - real* CIvec1old; - real* CIvec2old; - ivec SHbasis; - int CASelectrons; - int CASorbitals; -} t_QMrec; - -typedef struct -{ - int nrMMatoms; /* nr of MM atoms, updated every step*/ - rvec* xMM; /* shifted to center of box */ - int* indexMM; /* atom i = atom indexMM[I] in mdrun */ - real* MMcharges; /* MM point charges in std QMMM calc.*/ - int* shiftMM; - int* MMatomtype; /* only important for semi-emp. */ - real scalefactor; -} t_MMrec; - - -typedef struct t_QMMMrec -{ - int QMMMscheme; /* ONIOM (multi-layer) or normal */ - int nrQMlayers; /* number of QM layers (total layers +1 (MM)) */ - t_QMrec** qm; /* atoms and run params for each QM group */ - t_MMrec* mm; /* there can only be one MM subsystem ! */ -} t_QMMMrec; - -void atomic_number(int nr, char*** atomtype, int* nucnum); - -t_QMMMrec* mk_QMMMrec(); -/* allocates memory for QMMMrec */ - -void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr); - -/* init_QMMMrec initializes the QMMM record. From - * topology->atoms.atomname and topology->atoms.atomtype the atom - * names and types are read; from inputrec->QMcharge - * resp. inputrec->QMmult the nelecs and multiplicity are determined - * and md->cQMMM gives numbers of the MM and QM atoms - */ -void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, const t_mdatoms* md, const matrix box); - -/* update_QMMMrec fills the MM stuff in QMMMrec. The MM atoms are - * taken froom the neighbourlists of the QM atoms. In a QMMM run this - * routine should be called at every step, since it updates the MM - * elements of the t_QMMMrec struct. - */ -real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qmmm); - -/* QMMM computes the QM forces. This routine makes either function - * calls to gmx QM routines (derived from MOPAC7 (semi-emp.) and MPQC - * (ab initio)) or generates input files for an external QM package - * (listed in QMMMrec.QMpackage). The binary of the QM package is - * called by system(). - */ - -#endif diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 72a788c362..a700d2b673 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -79,7 +79,6 @@ #include "gromacs/mdlib/force_flags.h" #include "gromacs/mdlib/forcerec.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" -#include "gromacs/mdlib/qmmm.h" #include "gromacs/mdlib/update.h" #include "gromacs/mdlib/vsite.h" #include "gromacs/mdtypes/commrec.h" @@ -1520,12 +1519,6 @@ void do_force(FILE* fplog, } } - /* update QMMMrec, if necessary */ - if (fr->bQMMM) - { - update_QMMMrec(cr, fr, as_rvec_array(x.unpaddedArrayRef().data()), mdatoms, box); - } - // TODO Force flags should include haveFreeEnergyWork for this domain if (ddUsesGpuDirectCommunication && (domainWork.haveCpuBondedWork || domainWork.haveFreeEnergyWork)) { diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 65312657d8..4348a151de 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -97,9 +97,9 @@ #include "gromacs/mdlib/md_support.h" #include "gromacs/mdlib/mdatoms.h" #include "gromacs/mdlib/membed.h" -#include "gromacs/mdlib/qmmm.h" #include "gromacs/mdlib/sighandler.h" #include "gromacs/mdlib/stophandler.h" +#include "gromacs/mdlib/tgroup.h" #include "gromacs/mdlib/updategroups.h" #include "gromacs/mdlib/vsite.h" #include "gromacs/mdrun/mdmodules.h" diff --git a/src/gromacs/mdtypes/forcerec.h b/src/gromacs/mdtypes/forcerec.h index 874307c7ea..f093cc2fe1 100644 --- a/src/gromacs/mdtypes/forcerec.h +++ b/src/gromacs/mdtypes/forcerec.h @@ -252,13 +252,6 @@ struct t_forcerec */ int n_tpi = 0; - /* QMMM stuff */ - gmx_bool bQMMM = FALSE; - struct t_QMMMrec* qr = nullptr; - - /* QM-MM neighborlists */ - struct t_nblist* QMMMlist = nullptr; - /* Limit for printing large forces, negative is don't print */ real print_force = 0; -- 2.22.0