From: Viacheslav Bolnykh Date: Fri, 5 Jan 2018 17:16:44 +0000 (+0100) Subject: Implemented changes to preprocessor to work with MiMiC QM/MM X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=61e7a6a126acd5ed6397bd90a10fa6897b094e5d;p=alexxy%2Fgromacs.git Implemented changes to preprocessor to work with MiMiC QM/MM Added new integrator "mimic" and updated usage of QMMM-grps to be used with MiMiC QM/MM. Updated code zeros charges within QM region and changes bonds to CONNBOND and updates exclusion lists to exclude non-bonded interactions between quantum atoms Refs #2308 Change-Id: I34e798ae0c54457dad8f9a8f220482468db34708 --- diff --git a/CMakeLists.txt b/CMakeLists.txt index 5782ea0a2e..d06b9416ac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -200,6 +200,9 @@ gmx_dependent_option( ON GMX_MPI) mark_as_advanced(GMX_MPI_IN_PLACE) + +option(GMX_MIMIC "Enable MiMiC QM/MM interface (CPMD is required)" OFF) + option(GMX_FAHCORE "Build a library with mdrun functionality" OFF) mark_as_advanced(GMX_FAHCORE) diff --git a/cmake/gmxManageMimic.cmake b/cmake/gmxManageMimic.cmake new file mode 100644 index 0000000000..c0c98980bf --- /dev/null +++ b/cmake/gmxManageMimic.cmake @@ -0,0 +1,39 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright (c) 2017,2018, 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. + +if(GMX_MIMIC) + find_package(MimicCommLib REQUIRED) + include_directories(SYSTEM ${MimicCommLib_INCLUDE_DIRS}) + list(APPEND GMX_COMMON_LIBRARIES Upstream::mimiccomm) +endif() diff --git a/src/config.h.cmakein b/src/config.h.cmakein index e03dd3c4d4..0beccbcedc 100644 --- a/src/config.h.cmakein +++ b/src/config.h.cmakein @@ -393,6 +393,9 @@ /* Build using clang analyzer */ #cmakedefine01 GMX_CLANG_ANALYZER +/* Use MiMiC QM/MM interface */ +#cmakedefine01 GMX_MIMIC + /*! \endcond */ #endif diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 29b0208b03..5e40fe03a9 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -121,6 +121,7 @@ enum tpxv { tpxv_AcceleratedWeightHistogram, /**< sampling with accelerated weight histogram method (AWH) */ tpxv_RemoveImplicitSolvation, /**< removed support for implicit solvation */ tpxv_PullPrevStepCOMAsReference, /**< Enabled using the COM of the pull group of the last frame as reference for PBC */ + tpxv_MimicQMMM, /**< Inroduced support for MiMiC QM/MM interface */ tpxv_Count /**< the total number of tpxv versions */ }; @@ -1617,7 +1618,7 @@ static void do_inputrec(t_fileio *fio, t_inputrec *ir, gmx_bool bRead, snew(ir->opts.SAoff, ir->opts.ngQM); snew(ir->opts.SAsteps, ir->opts.ngQM); } - if (ir->opts.ngQM > 0) + if (ir->opts.ngQM > 0 && ir->bQMMM) { gmx_fio_ndo_int(fio, ir->opts.QMmethod, ir->opts.ngQM); gmx_fio_ndo_int(fio, ir->opts.QMbasis, ir->opts.ngQM); diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 5357d64b1f..2485d5a2e2 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -2146,10 +2146,15 @@ int gmx_grompp(int argc, char *argv[]) } else { - generate_qmexcl(&sys, ir, wi); + generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL); } } + if (ir->eI == eiMimic) + { + generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC); + } + if (ftp2bSet(efTRN, NFILE, fnm)) { if (bVerbose) diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index aff817aa61..d7720fda07 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -3576,57 +3576,75 @@ void do_index(const char* mdparin, const char *ndx, auto qmGroupNames = gmx::splitString(is->QMMM); auto qmMethods = gmx::splitString(is->QMmethod); auto qmBasisSets = gmx::splitString(is->QMbasis); - if (qmMethods.size() != qmGroupNames.size() || - qmBasisSets.size() != qmGroupNames.size()) + if (ir->eI != eiMimic) + { + if (qmMethods.size() != qmGroupNames.size() || + qmBasisSets.size() != qmGroupNames.size()) + { + gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets" + " and %zu methods\n", qmGroupNames.size(), + qmBasisSets.size(), qmMethods.size()); + } + /* group rest, if any, is always MM! */ + do_numbering(natoms, groups, qmGroupNames, grps, gnames, egcQMMM, + restnm, egrptpALL_GENREST, bVerbose, wi); + nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/ + ir->opts.ngQM = qmGroupNames.size(); + snew(ir->opts.QMmethod, nr); + snew(ir->opts.QMbasis, nr); + for (i = 0; i < nr; i++) + { + /* input consists of strings: RHF CASSCF PM3 .. These need to be + * converted to the corresponding enum in names.c + */ + ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(), + eQMmethodNR, + eQMmethod_names); + ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(), + eQMbasisNR, + eQMbasis_names); + + } + auto qmMultiplicities = gmx::splitString(is->QMmult); + auto qmCharges = gmx::splitString(is->QMcharge); + auto qmbSH = gmx::splitString(is->bSH); + snew(ir->opts.QMmult, nr); + snew(ir->opts.QMcharge, nr); + snew(ir->opts.bSH, nr); + convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult); + convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge); + convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH); + + auto CASelectrons = gmx::splitString(is->CASelectrons); + auto CASorbitals = gmx::splitString(is->CASorbitals); + snew(ir->opts.CASelectrons, nr); + snew(ir->opts.CASorbitals, nr); + convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons); + convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals); + + auto SAon = gmx::splitString(is->SAon); + auto SAoff = gmx::splitString(is->SAoff); + auto SAsteps = gmx::splitString(is->SAsteps); + snew(ir->opts.SAon, nr); + snew(ir->opts.SAoff, nr); + snew(ir->opts.SAsteps, nr); + convertInts(wi, SAon, "SAon", ir->opts.SAon); + convertInts(wi, SAoff, "SAoff", ir->opts.SAoff); + convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps); + } + else { - gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets" - " and %zu methods\n", qmGroupNames.size(), - qmBasisSets.size(), qmMethods.size()); + /* MiMiC */ + if (qmGroupNames.size() > 1) + { + gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported"); + } + /* group rest, if any, is always MM! */ + do_numbering(natoms, groups, qmGroupNames, grps, gnames, egcQMMM, + restnm, egrptpALL_GENREST, bVerbose, wi); + + ir->opts.ngQM = qmGroupNames.size(); } - /* group rest, if any, is always MM! */ - do_numbering(natoms, groups, qmGroupNames, grps, gnames, egcQMMM, - restnm, egrptpALL_GENREST, bVerbose, wi); - nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/ - ir->opts.ngQM = qmGroupNames.size(); - snew(ir->opts.QMmethod, nr); - snew(ir->opts.QMbasis, nr); - for (i = 0; i < nr; i++) - { - /* input consists of strings: RHF CASSCF PM3 .. These need to be - * converted to the corresponding enum in names.c - */ - ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(), eQMmethodNR, - eQMmethod_names); - ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(), eQMbasisNR, - eQMbasis_names); - - } - auto qmMultiplicities = gmx::splitString(is->QMmult); - auto qmCharges = gmx::splitString(is->QMcharge); - auto qmbSH = gmx::splitString(is->bSH); - snew(ir->opts.QMmult, nr); - snew(ir->opts.QMcharge, nr); - snew(ir->opts.bSH, nr); - convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult); - convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge); - convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH); - - auto CASelectrons = gmx::splitString(is->CASelectrons); - auto CASorbitals = gmx::splitString(is->CASorbitals); - snew(ir->opts.CASelectrons, nr); - snew(ir->opts.CASorbitals, nr); - convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons); - convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals); - - auto SAon = gmx::splitString(is->SAon); - auto SAoff = gmx::splitString(is->SAoff); - auto SAsteps = gmx::splitString(is->SAsteps); - snew(ir->opts.SAon, nr); - snew(ir->opts.SAoff, nr); - snew(ir->opts.SAsteps, nr); - convertInts(wi, SAon, "SAon", ir->opts.SAon); - convertInts(wi, SAoff, "SAoff", ir->opts.SAoff); - convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps); /* end of QMMM input */ diff --git a/src/gromacs/gmxpreprocess/tests/readir.cpp b/src/gromacs/gmxpreprocess/tests/readir.cpp index e4645b9c1f..5bff568412 100644 --- a/src/gromacs/gmxpreprocess/tests/readir.cpp +++ b/src/gromacs/gmxpreprocess/tests/readir.cpp @@ -209,5 +209,11 @@ TEST_F(GetIrTest, ImplicitSolventYesWorks) EXPECT_DEATH_IF_SUPPORTED(runTest(inputMdpFile), "Invalid enum"); } +TEST_F(GetIrTest, HandlesMimic) +{ + const char *inputMdpFile[] = {"integrator = mimic", "QMMM-grps = QMatoms"}; + runTest(joinStrings(inputMdpFile, "\n")); +} + } // namespace test } // namespace gmx diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesMimic.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesMimic.xml new file mode 100644 index 0000000000..55c0ed7e3b --- /dev/null +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesMimic.xml @@ -0,0 +1,321 @@ + + + + false + +; VARIOUS PREPROCESSING OPTIONS +; Preprocessor information: use cpp syntax. +; e.g.: -I/home/joe/doe -I/home/mary/roe +include = +; e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive) +define = + +; RUN CONTROL PARAMETERS +integrator = mimic +; Start time and timestep in ps +tinit = 0 +dt = 0.001 +nsteps = 0 +; For exact run continuation or redoing part of a run +init-step = 0 +; Part index is updated automatically on checkpointing (keeps files separate) +simulation-part = 1 +; mode for center of mass motion removal +comm-mode = Linear +; number of steps for center of mass motion removal +nstcomm = 100 +; group(s) for center of mass motion removal +comm-grps = + +; LANGEVIN DYNAMICS OPTIONS +; Friction coefficient (amu/ps) and random seed +bd-fric = 0 +ld-seed = -1 + +; ENERGY MINIMIZATION OPTIONS +; Force tolerance and initial step-size +emtol = 10 +emstep = 0.01 +; Max number of iterations in relax-shells +niter = 20 +; Step size (ps^2) for minimization of flexible constraints +fcstep = 0 +; Frequency of steepest descents steps when doing CG +nstcgsteep = 1000 +nbfgscorr = 10 + +; TEST PARTICLE INSERTION OPTIONS +rtpi = 0.05 + +; OUTPUT CONTROL OPTIONS +; Output frequency for coords (x), velocities (v) and forces (f) +nstxout = 0 +nstvout = 0 +nstfout = 0 +; Output frequency for energies to log file and energy file +nstlog = 1000 +nstcalcenergy = 100 +nstenergy = 1000 +; Output frequency and precision for .xtc file +nstxout-compressed = 0 +compressed-x-precision = 1000 +; This selects the subset of atoms for the compressed +; trajectory file. You can select multiple groups. By +; default, all atoms will be written. +compressed-x-grps = +; Selection of energy groups +energygrps = + +; NEIGHBORSEARCHING PARAMETERS +; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups) +cutoff-scheme = Verlet +; nblist update frequency +nstlist = 10 +; ns algorithm (simple or grid) +ns-type = Grid +; Periodic boundary conditions: xyz, no, xy +pbc = xyz +periodic-molecules = no +; Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom, +; a value of -1 means: use rlist +verlet-buffer-tolerance = 0.005 +; nblist cut-off +rlist = 1 +; long-range cut-off for switched potentials + +; OPTIONS FOR ELECTROSTATICS AND VDW +; Method for doing electrostatics +coulombtype = Cut-off +coulomb-modifier = Potential-shift-Verlet +rcoulomb-switch = 0 +rcoulomb = 1 +; Relative dielectric constant for the medium and the reaction field +epsilon-r = 1 +epsilon-rf = 0 +; Method for doing Van der Waals +vdw-type = Cut-off +vdw-modifier = Potential-shift-Verlet +; cut-off lengths +rvdw-switch = 0 +rvdw = 1 +; Apply long range dispersion corrections for Energy and Pressure +DispCorr = No +; Extension of the potential lookup tables beyond the cut-off +table-extension = 1 +; Separate tables between energy group pairs +energygrp-table = +; Spacing for the PME/PPPM FFT grid +fourierspacing = 0.12 +; FFT grid size, when a value is 0 fourierspacing will be used +fourier-nx = 0 +fourier-ny = 0 +fourier-nz = 0 +; EWALD/PME/PPPM parameters +pme-order = 4 +ewald-rtol = 1e-05 +ewald-rtol-lj = 0.001 +lj-pme-comb-rule = Geometric +ewald-geometry = 3d +epsilon-surface = 0 +implicit-solvent = no + +; OPTIONS FOR WEAK COUPLING ALGORITHMS +; Temperature coupling +tcoupl = No +nsttcouple = -1 +nh-chain-length = 10 +print-nose-hoover-chain-variables = no +; Groups to couple separately +tc-grps = +; Time constant (ps) and reference temperature (K) +tau-t = +ref-t = +; pressure coupling +pcoupl = No +pcoupltype = Isotropic +nstpcouple = -1 +; Time constant (ps), compressibility (1/bar) and reference P (bar) +tau-p = 1 +compressibility = +ref-p = +; Scaling of reference coordinates, No, All or COM +refcoord-scaling = No + +; OPTIONS FOR QMMM calculations +QMMM = no +; Groups treated Quantum Mechanically +QMMM-grps = QMatoms +; QM method +QMmethod = +; QMMM scheme +QMMMscheme = normal +; QM basisset +QMbasis = +; QM charge +QMcharge = +; QM multiplicity +QMmult = +; Surface Hopping +SH = +; CAS space options +CASorbitals = +CASelectrons = +SAon = +SAoff = +SAsteps = +; Scale factor for MM charges +MMChargeScaleFactor = 1 + +; SIMULATED ANNEALING +; Type of annealing for each temperature group (no/single/periodic) +annealing = +; Number of time points to use for specifying annealing in each group +annealing-npoints = +; List of times at the annealing points for each group +annealing-time = +; Temp. at each annealing point, for each group. +annealing-temp = + +; GENERATE VELOCITIES FOR STARTUP RUN +gen-vel = no +gen-temp = 300 +gen-seed = -1 + +; OPTIONS FOR BONDS +constraints = none +; Type of constraint algorithm +constraint-algorithm = Lincs +; Do not constrain the start configuration +continuation = no +; Use successive overrelaxation to reduce the number of shake iterations +Shake-SOR = no +; Relative tolerance of shake +shake-tol = 0.0001 +; Highest order in the expansion of the constraint coupling matrix +lincs-order = 4 +; Number of iterations in the final step of LINCS. 1 is fine for +; normal simulations, but use 2 to conserve energy in NVE runs. +; For energy minimization with constraints it should be 4 to 8. +lincs-iter = 1 +; Lincs will write a warning to the stderr if in one step a bond +; rotates over more degrees than +lincs-warnangle = 30 +; Convert harmonic bonds to morse potentials +morse = no + +; ENERGY GROUP EXCLUSIONS +; Pairs of energy groups for which all non-bonded interactions are excluded +energygrp-excl = + +; WALLS +; Number of walls, type, atom types, densities and box-z scale factor for Ewald +nwall = 0 +wall-type = 9-3 +wall-r-linpot = -1 +wall-atomtype = +wall-density = +wall-ewald-zfac = 3 + +; COM PULLING +pull = no + +; AWH biasing +awh = no + +; ENFORCED ROTATION +; Enforced rotation: No or Yes +rotation = no + +; Group to display and/or manipulate in interactive MD session +IMD-group = + +; NMR refinement stuff +; Distance restraints type: No, Simple or Ensemble +disre = No +; Force weighting of pairs in one distance restraint: Conservative or Equal +disre-weighting = Conservative +; Use sqrt of the time averaged times the instantaneous violation +disre-mixed = no +disre-fc = 1000 +disre-tau = 0 +; Output frequency for pair distances to energy file +nstdisreout = 100 +; Orientation restraints: No or Yes +orire = no +; Orientation restraints force constant and tau for time averaging +orire-fc = 0 +orire-tau = 0 +orire-fitgrp = +; Output frequency for trace(SD) and S to energy file +nstorireout = 100 + +; Free energy variables +free-energy = no +couple-moltype = +couple-lambda0 = vdw-q +couple-lambda1 = vdw-q +couple-intramol = no +init-lambda = -1 +init-lambda-state = -1 +delta-lambda = 0 +nstdhdl = 50 +fep-lambdas = +mass-lambdas = +coul-lambdas = +vdw-lambdas = +bonded-lambdas = +restraint-lambdas = +temperature-lambdas = +calc-lambda-neighbors = 1 +init-lambda-weights = +dhdl-print-energy = no +sc-alpha = 0 +sc-power = 1 +sc-r-power = 6 +sc-sigma = 0.3 +sc-coul = no +separate-dhdl-file = yes +dhdl-derivatives = yes +dh_hist_size = 0 +dh_hist_spacing = 0.1 + +; Non-equilibrium MD stuff +acc-grps = +accelerate = +freezegrps = +freezedim = +cos-acceleration = 0 +deform = + +; simulated tempering variables +simulated-tempering = no +simulated-tempering-scaling = geometric +sim-temp-low = 300 +sim-temp-high = 300 + +; Ion/water position swapping for computational electrophysiology setups +; Swap positions along direction: no, X, Y, Z +swapcoords = no +adress = no + +; User defined thingies +user1-grps = +user2-grps = +userint1 = 0 +userint2 = 0 +userint3 = 0 +userint4 = 0 +userreal1 = 0 +userreal2 = 0 +userreal3 = 0 +userreal4 = 0 +; Electric fields +; Format for electric-field-x, etc. is: four real variables: +; amplitude (V/nm), frequency omega (1/ps), time for the pulse peak (ps), +; and sigma (ps) width of the pulse. Omega = 0 means static field, +; sigma = 0 means no pulse, leaving the field to be a cosine function. +electric-field-x = 0 0 0 0 +electric-field-y = 0 0 0 0 +electric-field-z = 0 0 0 0 + + diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index e753c04d7d..e07b0f7d20 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -1040,9 +1040,22 @@ char **do_top(bool bVerbose, return title; } - +/*! \brief + * Generate exclusion lists for QM/MM. + * + * This routine updates the exclusion lists for QM atoms in order to include all other QM + * atoms of this molecule. Moreover, this routine replaces bonds between QM atoms with + * CONNBOND and, when MiMiC is not used, removes bonded interactions between QM and link atoms. + * Finally, in case if MiMiC QM/MM is used - charges of QM atoms are set to 0 + * + * @param molt molecule type with QM atoms + * @param grpnr group informatio + * @param ir input record + * @param wi warning handler + * @param qmmmMode QM/MM mode switch: original/MiMiC + */ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *grpnr, - t_inputrec *ir, warninp_t wi) + t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode) { /* This routine expects molt->ilist to be of size F_NRE and ordered. */ @@ -1077,7 +1090,9 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr } if ((grpnr ? grpnr[i] : 0) == j) { - qm_arr[qm_nr++] = i; + qm_arr[qm_nr++] = i; + molt->atoms.atom[i].q = 0.0; + molt->atoms.atom[i].qB = 0.0; } } } @@ -1173,7 +1188,17 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr numQmAtoms++; } } - bexcl = (numQmAtoms >= nratoms - 1); + + /* MiMiC treats link atoms as quantum atoms - therefore + * we do not need do additional exclusions here */ + if (qmmmMode == GmxQmmmMode::GMX_QMMM_MIMIC) + { + bexcl = numQmAtoms == nratoms; + } + else + { + bexcl = (numQmAtoms >= nratoms - 1); + } if (bexcl && ftype == F_SETTLE) { @@ -1204,32 +1229,35 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr * linkatoms interaction with the QMatoms and would be counted * twice. */ - for (int i = 0; i < F_NRE; i++) + if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC) { - if (IS_CHEMBOND(i)) + for (int i = 0; i < F_NRE; i++) { - int j = 0; - while (j < molt->ilist[i].size()) + if (IS_CHEMBOND(i)) { - int a1 = molt->ilist[i].iatoms[j+1]; - int a2 = molt->ilist[i].iatoms[j+2]; - if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2])) + int j = 0; + while (j < molt->ilist[i].size()) { - if (link_nr >= link_max) - { - link_max += 10; - srenew(link_arr, link_max); - } - if (bQMMM[a1]) + int a1 = molt->ilist[i].iatoms[j + 1]; + int a2 = molt->ilist[i].iatoms[j + 2]; + if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2])) { - link_arr[link_nr++] = a2; - } - else - { - link_arr[link_nr++] = a1; + if (link_nr >= link_max) + { + link_max += 10; + srenew(link_arr, link_max); + } + if (bQMMM[a1]) + { + link_arr[link_nr++] = a2; + } + else + { + link_arr[link_nr++] = a1; + } } + j += 3; } - j += 3; } } } @@ -1238,9 +1266,13 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr { blink[i] = FALSE; } - for (int i = 0; i < link_nr; i++) + + if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC) { - blink[link_arr[i]] = TRUE; + for (int i = 0; i < link_nr; i++) + { + blink[link_arr[i]] = TRUE; + } } /* creating the exclusion block for the QM atoms. Each QM atom has * as excluded elements all the other QMatoms (and itself). @@ -1327,15 +1359,17 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr free(blink); } /* generate_qmexcl */ -void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi) +void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode) { /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered. */ - unsigned char *grpnr; - int mol, nat_mol, nr_mol_with_qm_atoms = 0; - gmx_molblock_t *molb; - bool bQMMM; + unsigned char *grpnr; + int mol, nat_mol, nr_mol_with_qm_atoms = 0; + gmx_molblock_t *molb; + bool bQMMM; + int index_offset = 0; + int qm_nr = 0; grpnr = sys->groups.grpnr[egcQMMM]; @@ -1348,11 +1382,13 @@ void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi) bQMMM = FALSE; for (int i = 0; i < nat_mol; i++) { - if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM) + if ((grpnr ? grpnr[i] : 0) < (ir->opts.ngQM)) { - bQMMM = TRUE; + bQMMM = TRUE; + qm_nr++; } } + if (bQMMM) { nr_mol_with_qm_atoms++; @@ -1379,25 +1415,39 @@ void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi) sys->molblock[mb+1].nmol -= 1; } - /* Add a moltype for the QMMM molecule */ - sys->moltype.push_back(sys->moltype[molb->type]); + /* Create a copy of a moltype for a molecule + * containing QM atoms and append it in the end of the list + */ + std::vector temp(sys->moltype.size()); + for (size_t i = 0; i < sys->moltype.size(); ++i) + { + copy_moltype(&sys->moltype[i], &temp[i]); + } + sys->moltype.resize(sys->moltype.size() + 1); + for (size_t i = 0; i < temp.size(); ++i) + { + copy_moltype(&temp[i], &sys->moltype[i]); + } + copy_moltype(&sys->moltype[molb->type], &sys->moltype.back()); /* Copy the exclusions to a new array, since this is the only * thing that needs to be modified for QMMM. */ - copy_blocka(&sys->moltype[molb->type ].excls, + copy_blocka(&sys->moltype[molb->type].excls, &sys->moltype.back().excls); /* Set the molecule type for the QMMM molblock */ molb->type = sys->moltype.size() - 1; } - generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, wi); + generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, wi, qmmmMode); } if (grpnr) { grpnr += nat_mol; } + index_offset += nat_mol; } } - if (nr_mol_with_qm_atoms > 1) + if (qmmmMode == GmxQmmmMode::GMX_QMMM_ORIGINAL && + nr_mol_with_qm_atoms > 1) { /* generate a warning is there are QM atoms in different * topologies. In this case it is not possible at this stage to diff --git a/src/gromacs/gmxpreprocess/topio.h b/src/gromacs/gmxpreprocess/topio.h index 8d8ce58575..fd9c620cd5 100644 --- a/src/gromacs/gmxpreprocess/topio.h +++ b/src/gromacs/gmxpreprocess/topio.h @@ -48,6 +48,7 @@ struct gmx_mtop_t; struct t_gromppopts; struct t_inputrec; struct warninp; +enum struct GmxQmmmMode; typedef warninp *warninp_t; double check_mol(const gmx_mtop_t *mtop, warninp_t wi); @@ -73,6 +74,6 @@ char **do_top(bool bVerbose, warninp_t wi); /* This routine expects sys->molt[m].ilist to be of size F_NRE and ordered. */ -void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi); +void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode); #endif diff --git a/src/gromacs/mdtypes/md_enums.cpp b/src/gromacs/mdtypes/md_enums.cpp index de4e0c362d..3b03645e3f 100644 --- a/src/gromacs/mdtypes/md_enums.cpp +++ b/src/gromacs/mdtypes/md_enums.cpp @@ -56,8 +56,6 @@ const char *yesno_names[BOOL_NR+1] = "no", "yes", nullptr }; - - const char *ens_names[ensNR+1] = { "Grid", "Simple", nullptr @@ -65,7 +63,7 @@ const char *ens_names[ensNR+1] = const char *ei_names[eiNR+1] = { - "md", "steep", "cg", "bd", "sd2 - removed", "nm", "l-bfgs", "tpi", "tpic", "sd", "md-vv", "md-vv-avek", nullptr + "md", "steep", "cg", "bd", "sd2 - removed", "nm", "l-bfgs", "tpi", "tpic", "sd", "md-vv", "md-vv-avek", "mimic", nullptr }; const char *ecutscheme_names[ecutsNR+1] = { diff --git a/src/gromacs/mdtypes/md_enums.h b/src/gromacs/mdtypes/md_enums.h index af55870361..f5cf68e68e 100644 --- a/src/gromacs/mdtypes/md_enums.h +++ b/src/gromacs/mdtypes/md_enums.h @@ -229,16 +229,18 @@ extern const char *ens_names[ensNR+1]; * and the half step kinetic energy for temperature control */ enum { - eiMD, eiSteep, eiCG, eiBD, eiSD2_REMOVED, eiNM, eiLBFGS, eiTPI, eiTPIC, eiSD1, eiVV, eiVVAK, eiNR + eiMD, eiSteep, eiCG, eiBD, eiSD2_REMOVED, eiNM, eiLBFGS, eiTPI, eiTPIC, eiSD1, eiVV, eiVVAK, eiMimic, eiNR }; //! Name of the integrator algorithm extern const char *ei_names[eiNR+1]; //! Macro returning integrator string #define EI(e) enum_name(e, eiNR, ei_names) +//! Do we use MiMiC QM/MM? +#define EI_MIMIC(e) ((e) == eiMimic) //! Do we use velocity Verlet #define EI_VV(e) ((e) == eiVV || (e) == eiVVAK) //! Do we use molecular dynamics -#define EI_MD(e) ((e) == eiMD || EI_VV(e)) +#define EI_MD(e) ((e) == eiMD || EI_VV(e) || EI_MIMIC(e)) //! Do we use stochastic dynamics #define EI_SD(e) ((e) == eiSD1) //! Do we use any stochastic integrator @@ -638,4 +640,10 @@ enum gmx_nblist_interaction_type //! String corresponding to interactions in neighborlist code extern const char *gmx_nblist_interaction_names[GMX_NBLIST_INTERACTION_NR+1]; + +//! \brief QM/MM mode +enum struct GmxQmmmMode { + GMX_QMMM_ORIGINAL, + GMX_QMMM_MIMIC +}; #endif /* GMX_MDTYPES_MD_ENUMS_H */ diff --git a/src/gromacs/topology/atoms.cpp b/src/gromacs/topology/atoms.cpp index 7e92cf676b..fe3c390cd1 100644 --- a/src/gromacs/topology/atoms.cpp +++ b/src/gromacs/topology/atoms.cpp @@ -226,7 +226,12 @@ t_atoms *copy_t_atoms(const t_atoms *src) dst->atomtypeB[i] = src->atomtypeB[i]; } } - dst->nres = src->nres; + dst->haveBState = src->haveBState; + dst->haveCharge = src->haveCharge; + dst->haveMass = src->haveMass; + dst->havePdbInfo = src->havePdbInfo; + dst->haveType = src->haveType; + dst->nres = src->nres; for (i = 0; (i < src->nres); i++) { dst->resinfo[i] = src->resinfo[i]; diff --git a/src/gromacs/topology/block.cpp b/src/gromacs/topology/block.cpp index b2351bd31c..b371453154 100644 --- a/src/gromacs/topology/block.cpp +++ b/src/gromacs/topology/block.cpp @@ -311,3 +311,18 @@ void pr_blocka(FILE *fp, int indent, const char *title, const t_blocka *block, g } } } + +void copy_block(const t_block *src, t_block *dst) +{ + dst->nr = src->nr; + /* Workaround for inconsistent handling of nalloc_index in + * other parts of the code. Often nalloc_index and nalloc_a + * are not set. + */ + dst->nalloc_index = std::max(src->nalloc_index, dst->nr + 1); + snew(dst->index, dst->nalloc_index); + for (int i = 0; i < dst->nr+1; ++i) + { + dst->index[i] = src->index[i]; + } +} diff --git a/src/gromacs/topology/block.h b/src/gromacs/topology/block.h index 3eb9d520f6..b8a47faf00 100644 --- a/src/gromacs/topology/block.h +++ b/src/gromacs/topology/block.h @@ -225,6 +225,8 @@ void done_blocka(t_blocka *block); void copy_blocka(const t_blocka *src, t_blocka *dest); +void copy_block(const t_block *src, t_block *dst); + void stupid_fill_block(t_block *grp, int natom, gmx_bool bOneIndexGroup); /* Fill a block structure with numbers identical to the index * (0, 1, 2, .. natom-1) diff --git a/src/gromacs/topology/idef.cpp b/src/gromacs/topology/idef.cpp index ac8fade143..1fde04c3ff 100644 --- a/src/gromacs/topology/idef.cpp +++ b/src/gromacs/topology/idef.cpp @@ -488,3 +488,16 @@ void done_idef(t_idef *idef) delete idef->cmap_grid; init_idef(idef); } + +void copy_ilist(const t_ilist *src, t_ilist *dst) +{ + dst->nr = src->nr; + dst->nr_nonperturbed = src->nr_nonperturbed; + dst->nalloc = src->nalloc; + + snew(dst->iatoms, dst->nr); + for (int i = 0; i < dst->nr; ++i) + { + dst->iatoms[i] = src->iatoms[i]; + } +} diff --git a/src/gromacs/topology/idef.h b/src/gromacs/topology/idef.h index 62a9fa6f0d..df9b6ffd97 100644 --- a/src/gromacs/topology/idef.h +++ b/src/gromacs/topology/idef.h @@ -392,4 +392,6 @@ void init_idef(t_idef *idef); */ void done_idef(t_idef *idef); +void copy_ilist(const t_ilist *src, t_ilist *dst); + #endif diff --git a/src/gromacs/topology/topology.cpp b/src/gromacs/topology/topology.cpp index 58f8490bbb..2dc6509752 100644 --- a/src/gromacs/topology/topology.cpp +++ b/src/gromacs/topology/topology.cpp @@ -639,3 +639,18 @@ int getGroupType(const gmx_groups_t *group, int type, int atom) { return (group->grpnr[type] ? group->grpnr[type][atom] : 0); } + +void copy_moltype(const gmx_moltype_t *src, gmx_moltype_t *dst) +{ + dst->name = src->name; + copy_blocka(&src->excls, &dst->excls); + copy_block(&src->cgs, &dst->cgs); + t_atoms *atomsCopy = copy_t_atoms(&src->atoms); + dst->atoms = *atomsCopy; + sfree(atomsCopy); + + for (int i = 0; i < F_NRE; ++i) + { + dst->ilist[i] = src->ilist[i]; + } +} diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index e4bfcf1696..3b2f1a8c52 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -227,4 +227,6 @@ void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1, //! Deleter for gmx_localtop_t, needed until it has a proper destructor. using ExpandedTopologyPtr = gmx::unique_cptr; +void copy_moltype(const gmx_moltype_t *src, gmx_moltype_t *dst); + #endif