From 483a7a112f7a2b5e1b87534280d5a77e90fdecd9 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Wed, 9 Nov 2016 22:25:53 +0100 Subject: [PATCH] Add mdp option check with decoupled modes When atoms involved in an angle with constrained bonds have very different masses, there can be very weakly coupled dynamics modes. Default mdp settings are often not sufficiently accurate to obtain equipartitioning. This change adds a grompp check for this issue. Part of #2071. Change-Id: I64323154c38abe23b8d38d9d539d2a713a5c35e0 --- src/gromacs/gmxpreprocess/grompp.cpp | 173 ++++++++++++++++++++++++++++++++++- 1 file changed, 172 insertions(+), 1 deletion(-) diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index cd4ee3af4b..be33e421c5 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017, 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. @@ -73,6 +73,7 @@ #include "gromacs/math/vec.h" #include "gromacs/mdlib/calc_verletbuf.h" #include "gromacs/mdlib/compute_io.h" +#include "gromacs/mdlib/constr.h" #include "gromacs/mdlib/genborn.h" #include "gromacs/mdlib/perf_est.h" #include "gromacs/mdtypes/inputrec.h" @@ -1421,6 +1422,174 @@ static void checkForUnboundAtoms(const gmx_mtop_t *mtop, } } +/*! \brief Checks if there are decoupled modes in moleculetype \p molt. + * + * The specific decoupled modes this routine check for are angle modes + * where the two bonds are constrained and the atoms a both ends are only + * involved in a single constraint; the mass of the two atoms needs to + * differ by more than \p massFactorThreshold. + */ +static bool haveDecoupledModeInMol(const gmx_moltype_t *molt, + const t_iparams *iparams, + real massFactorThreshold) +{ + if (molt->ilist[F_CONSTR].nr == 0 && molt->ilist[F_CONSTRNC].nr == 0) + { + return false; + } + + const t_atom * atom = molt->atoms.atom; + + int numFlexibleConstraints; + t_blocka atomToConstraints = make_at2con(0, molt->atoms.nr, + molt->ilist, iparams, + FALSE, + &numFlexibleConstraints); + + bool haveDecoupledMode = false; + for (int ftype = 0; ftype < F_NRE; ftype++) + { + if (interaction_function[ftype].flags & IF_ATYPE) + { + const int nral = NRAL(ftype); + const t_ilist *il = &molt->ilist[ftype]; + for (int i = 0; i < il->nr; i += 1 + nral) + { + /* Here we check for the mass difference between the atoms + * at both ends of the angle, that the atoms at the ends have + * 1 contraint and the atom in the middle at least 3; we check + * that the 3 atoms are linked by constraints below. + * We check for at least three constraints for the middle atom, + * since with only the two bonds in the angle, we have 3-atom + * molecule, which has much more thermal exhange in this single + * angle mode than molecules with more atoms. + * Note that this check also catches molecules where more atoms + * are connected to one or more atoms in the angle, but by + * bond potentials instead of angles. But such cases will not + * occur in "normal" molecules and it doesn't hurt running + * those with higher accuracy settings as well. + */ + int a0 = il->iatoms[1 + i]; + int a1 = il->iatoms[1 + i + 1]; + int a2 = il->iatoms[1 + i + 2]; + if ((atom[a0].m > atom[a2].m*massFactorThreshold || + atom[a2].m > atom[a0].m*massFactorThreshold) && + atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1 && + atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1 && + atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3) + { + int constraint0 = atomToConstraints.a[atomToConstraints.index[a0]]; + int constraint2 = atomToConstraints.a[atomToConstraints.index[a2]]; + + bool foundAtom0 = false; + bool foundAtom2 = false; + for (int conIndex = atomToConstraints.index[a1]; conIndex < atomToConstraints.index[a1 + 1]; conIndex++) + { + if (atomToConstraints.a[conIndex] == constraint0) + { + foundAtom0 = true; + } + if (atomToConstraints.a[conIndex] == constraint2) + { + foundAtom2 = true; + } + } + if (foundAtom0 && foundAtom2) + { + haveDecoupledMode = true; + } + } + } + } + } + + done_blocka(&atomToConstraints); + + return haveDecoupledMode; +} + +/*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes. + * + * When decoupled modes are present and the accuray in insufficient, + * this routine issues a warning if the accuracy is insufficient. + */ +static void checkDecoupledModeAccuracy(const gmx_mtop_t *mtop, + const t_inputrec *ir, + warninp_t wi) +{ + /* We only have issues with decoupled modes with normal MD. + * With stochastic dynamics equipartitioning is enforced strongly. + */ + if (!EI_MD(ir->eI)) + { + return; + } + + /* When atoms of very different mass are involved in an angle potential + * and both bonds in the angle are constrained, the dynamic modes in such + * angles have very different periods and significant energy exchange + * takes several nanoseconds. Thus even a small amount of error in + * different algorithms can lead to violation of equipartitioning. + * The parameters below are mainly based on an all-atom chloroform model + * with all bonds constrained. Equipartitioning is off by more than 1% + * (need to run 5-10 ns) when the difference in mass between H and Cl + * is more than a factor 13 and the accuracy is less than the thresholds + * given below. This has been verified on some other molecules. + * + * Note that the buffer and shake parameters have unit length and + * energy/time, respectively, so they will "only" work correctly + * for atomistic force fields using MD units. + */ + const real massFactorThreshold = 13.0; + const real bufferToleranceThreshold = 1e-4; + const int lincsIterationThreshold = 2; + const int lincsOrderThreshold = 4; + const real shakeToleranceThreshold = 0.005*ir->delta_t; + + bool lincsWithSufficientTolerance = (ir->eConstrAlg == econtLINCS && ir->nLincsIter >= lincsIterationThreshold && ir->nProjOrder >= lincsOrderThreshold); + bool shakeWithSufficientTolerance = (ir->eConstrAlg == econtSHAKE && ir->shake_tol <= 1.1*shakeToleranceThreshold); + if (ir->cutoff_scheme == ecutsVERLET && + ir->verletbuf_tol <= 1.1*bufferToleranceThreshold && + (lincsWithSufficientTolerance || shakeWithSufficientTolerance)) + { + return; + } + + bool haveDecoupledMode = false; + for (int mt = 0; mt < mtop->nmoltype; mt++) + { + if (haveDecoupledModeInMol(&mtop->moltype[mt], mtop->ffparams.iparams, + massFactorThreshold)) + { + haveDecoupledMode = true; + } + } + + if (haveDecoupledMode) + { + char modeMessage[STRLEN]; + sprintf(modeMessage, "There are atoms at both ends of an angle, connected by constraints and with masses that differ by more than a factor of %g. This means that there are likely dynamic modes that are only very weakly coupled.", + massFactorThreshold); + char buf[STRLEN]; + if (ir->cutoff_scheme == ecutsVERLET) + { + sprintf(buf, "%s To ensure good equipartitioning, you need to either not use constraints on all bonds (but, if possible, only on bonds involving hydrogens) or use integrator = %s or decrease one or more tolerances: verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order >= %d or SHAKE tolerance <= %g", + modeMessage, + ei_names[eiSD1], + bufferToleranceThreshold, + lincsIterationThreshold, lincsOrderThreshold, + shakeToleranceThreshold); + } + else + { + sprintf(buf, "%s To ensure good equipartitioning, we suggest to switch to the %s cutoff-scheme, since that allows for better control over the Verlet buffer size and thus over the energy drift.", + modeMessage, + ecutscheme_names[ecutsVERLET]); + } + warning(wi, buf); + } +} + static void set_verlet_buffer(const gmx_mtop_t *mtop, t_inputrec *ir, real buffer_temp, @@ -1889,6 +2058,8 @@ int gmx_grompp(int argc, char *argv[]) check_bonds_timestep(sys, ir->delta_t, wi); } + checkDecoupledModeAccuracy(sys, ir, wi); + if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps) { warning_note(wi, "Zero-step energy minimization will alter the coordinates before calculating the energy. If you just want the energy of a single point, try zero-step MD (with unconstrained_start = yes). To do multiple single-point energy evaluations of different configurations of the same topology, use mdrun -rerun."); -- 2.11.4.GIT