src/gromacs/mdlib/updategroups.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /* \internal \file
  36  *
  37  * \brief Implements functions for generating update groups
  38  *
  39  * \author Berk Hess <hess@kth.se>
  40  * \ingroup module_mdlib
  41  */
  42
  43 #include "gmxpre.h"
  44
  45 #include "updategroups.h"
  46
  47 #include <cmath>
  48
  49 #include "gromacs/mdlib/constr.h"
  50 #include "gromacs/mdtypes/inputrec.h"
  51 #include "gromacs/topology/block.h"
  52 #include "gromacs/topology/idef.h"
  53 #include "gromacs/topology/ifunc.h"
  54 #include "gromacs/topology/topology.h"
  55
  56 namespace gmx
  57 {
  58
  59 /*! \brief Returns whether \p moltype contains flexible constraints */
  60 static bool hasFlexibleConstraints(const gmx_moltype_t            &moltype,
  61                                    gmx::ArrayRef<const t_iparams>  iparams)
  62 {
  63     for (auto &ilist : extractILists(moltype.ilist, IF_CONSTRAINT))
  64     {
  65         if (ilist.functionType != F_SETTLE)
  66         {
  67             for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
  68             {
  69                 if (isConstraintFlexible(iparams.data(), ilist.iatoms[i]))
  70                 {
  71                     return true;
  72                 }
  73             }
  74         }
  75     }
  76
  77     return false;
  78 }
  79
  80 /*! \brief Returns whether moltype has incompatible vsites.
  81  *
  82  * For simplicity the only compatible vsites are linear 2 or 3 atom sites
  83  * that are constructed in between the 2 or 3 contructing atoms,
  84  */
  85 static bool hasIncompatibleVsites(const gmx_moltype_t            &moltype,
  86                                   gmx::ArrayRef<const t_iparams>  iparams)
  87 {
  88     bool hasIncompatibleVsites = false;
  89
  90     for (auto &ilist : extractILists(moltype.ilist, IF_VSITE))
  91     {
  92         if (ilist.functionType == F_VSITE2 || ilist.functionType == F_VSITE3)
  93         {
  94             for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
  95             {
  96                 const t_iparams &iparam = iparams[ilist.iatoms[i]];
  97                 real             coeffMin;
  98                 real             coeffSum;
  99                 if (ilist.functionType == F_VSITE2)
 100                 {
 101                     coeffMin = iparam.vsite.a;
 102                     coeffSum = iparam.vsite.a;
 103                 }
 104                 else
 105                 {
 106                     coeffMin = std::min(iparam.vsite.a, iparam.vsite.b);
 107                     coeffSum = iparam.vsite.a + iparam.vsite.b;
 108                 }
 109                 if (coeffMin < 0 || coeffSum > 1)
 110                 {
 111                     hasIncompatibleVsites = true;
 112                     break;
 113                 }
 114             }
 115         }
 116         else
 117         {
 118             hasIncompatibleVsites = true;
 119             break;
 120         }
 121     }
 122
 123     return hasIncompatibleVsites;
 124 }
 125
 126 /*! \brief Returns a merged list with constraints of all types */
 127 static InteractionList jointConstraintList(const gmx_moltype_t &moltype)
 128 {
 129     InteractionList   ilistCombined;
 130     std::vector<int> &iatoms = ilistCombined.iatoms;
 131
 132     for (auto &ilist : extractILists(moltype.ilist, IF_CONSTRAINT))
 133     {
 134         if (ilist.functionType == F_SETTLE)
 135         {
 136             for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
 137             {
 138                 iatoms.push_back(-1);
 139                 iatoms.push_back(ilist.iatoms[i + 1]);
 140                 iatoms.push_back(ilist.iatoms[i + 2]);
 141                 iatoms.push_back(-1);
 142                 iatoms.push_back(ilist.iatoms[i + 1]);
 143                 iatoms.push_back(ilist.iatoms[i + 3]);
 144                 iatoms.push_back(-1);
 145                 iatoms.push_back(ilist.iatoms[i + 2]);
 146                 iatoms.push_back(ilist.iatoms[i + 3]);
 147             }
 148         }
 149         else
 150         {
 151             GMX_RELEASE_ASSERT(NRAL(ilist.functionType) == 2, "Can only handle two-atom non-SETTLE constraints");
 152
 153             iatoms.insert(iatoms.end(),
 154                           ilist.iatoms.begin(), ilist.iatoms.end());
 155         }
 156     }
 157
 158     return ilistCombined;
 159 }
 160
 161 /*! \brief Struct for returning an atom range */
 162 struct AtomIndexExtremes
 163 {
 164     int minAtom; //!< The minimum atom index
 165     int maxAtom; //!< The maximum atom index
 166 };
 167
 168 /*! \brief Returns the range of constructing atom for vsite with atom index \p a */
 169 static AtomIndexExtremes
 170 vsiteConstructRange(int                  a,
 171                     const gmx_moltype_t &moltype)
 172 {
 173     AtomIndexExtremes extremes = { -1, -1 };
 174
 175     for (auto &ilist : extractILists(moltype.ilist, IF_VSITE))
 176     {
 177         for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
 178         {
 179             if (ilist.iatoms[i + 1] == a)
 180             {
 181                 extremes.minAtom = ilist.iatoms[i + 2];
 182                 extremes.maxAtom = ilist.iatoms[i + 2];
 183                 for (size_t j = i + 3; j < i + ilistStride(ilist); j++)
 184                 {
 185                     extremes.minAtom = std::min(extremes.minAtom, ilist.iatoms[j]);
 186                     extremes.maxAtom = std::max(extremes.maxAtom, ilist.iatoms[j]);
 187                 }
 188                 return extremes;
 189             }
 190         }
 191     }
 192
 193     GMX_RELEASE_ASSERT(false, "If a is a vsite, we should have found constructing atoms");
 194
 195     return extremes;
 196 }
 197
 198 /*! \brief Returns the range of atoms constrained to atom \p a (including \p a itself) */
 199 static AtomIndexExtremes
 200 constraintAtomRange(int                    a,
 201                     const t_blocka        &at2con,
 202                     const InteractionList &ilistConstraints)
 203 {
 204     AtomIndexExtremes extremes = { a, a };
 205
 206     for (int i = at2con.index[a]; i < at2con.index[a + 1]; i++)
 207     {
 208         for (int j = 0; j < 2; j++)
 209         {
 210             int atomJ        = ilistConstraints.iatoms[at2con.a[i]*3 + 1 + j];
 211             extremes.minAtom = std::min(extremes.minAtom, atomJ);
 212             extremes.maxAtom = std::max(extremes.maxAtom, atomJ);
 213         }
 214     }
 215
 216     return extremes;
 217 }
 218
 219 /*! \brief Returns a list that tells whether atoms in \p moltype are vsites */
 220 static std::vector<bool> buildIsParticleVsite(const gmx_moltype_t &moltype)
 221 {
 222     std::vector<bool> isVsite(moltype.atoms.nr);
 223
 224     for (auto &ilist : extractILists(moltype.ilist, IF_VSITE))
 225     {
 226         for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
 227         {
 228             int vsiteAtom      = ilist.iatoms[i + 1];
 229             isVsite[vsiteAtom] = true;
 230         }
 231     }
 232
 233     return isVsite;
 234 }
 235
 236 /*! \brief Returns the size of the update group starting at \p firstAtom or 0 when criteria (see updategroups.h) are not met */
 237 static int detectGroup(int                    firstAtom,
 238                        const gmx_moltype_t   &moltype,
 239                        const t_blocka        &at2con,
 240                        const InteractionList &ilistConstraints)
 241 {
 242     /* We should be using moltype.atoms.atom[].ptype for checking whether
 243      * a particle is a vsite. But the test code can't fill t_atoms,
 244      * because it uses C pointers which get double freed.
 245      */
 246     std::vector<bool> isParticleVsite = buildIsParticleVsite(moltype);
 247
 248     /* A non-vsite atom without constraints is an update group by itself */
 249     if (!isParticleVsite[firstAtom] &&
 250         at2con.index[firstAtom + 1] - at2con.index[firstAtom] == 0)
 251     {
 252         return 1;
 253     }
 254
 255     /* A (potential) update group starts at firstAtom and should consist
 256      * of two or more atoms and possibly vsites. At least one atom should
 257      * have constraints with all other atoms and vsites should have all
 258      * constructing atoms inside the group. Here we increase lastAtom until
 259      * the criteria are fulfilled or exit when criteria cannot be met.
 260      */
 261     int numAtomsWithConstraints = 0;
 262     int maxConstraintsPerAtom   = 0;
 263     int lastAtom                = firstAtom;
 264     int a                       = firstAtom;
 265     while (a <= lastAtom)
 266     {
 267         if (isParticleVsite[a])
 268         {
 269             AtomIndexExtremes extremes = vsiteConstructRange(a, moltype);
 270             if (extremes.minAtom < firstAtom)
 271             {
 272                 /* A constructing atom is outside the group,
 273                  * we can not use update groups.
 274                  */
 275                 return 0;
 276             }
 277             lastAtom = std::max(lastAtom, extremes.maxAtom);
 278         }
 279         else
 280         {
 281             int numConstraints = at2con.index[a + 1] - at2con.index[a];
 282             if (numConstraints == 0)
 283             {
 284                 /* We can not have unconstrained atoms in an update group */
 285                 return 0;
 286             }
 287             /* This atom has at least one constraint.
 288              * Check whether all constraints are within the group
 289              * and whether we need to extend the group.
 290              */
 291             numAtomsWithConstraints += 1;
 292             maxConstraintsPerAtom    = std::max(maxConstraintsPerAtom, numConstraints);
 293
 294             AtomIndexExtremes extremes =
 295                 constraintAtomRange(a, at2con, ilistConstraints);
 296             if (extremes.minAtom < firstAtom)
 297             {
 298                 /* Constraint to atom outside the "group" */
 299                 return 0;
 300             }
 301             lastAtom = std::max(lastAtom, extremes.maxAtom);
 302         }
 303
 304         a++;
 305     }
 306
 307     /* lastAtom might be followed by a vsite that is constructed from atoms
 308      * with index <= lastAtom. Handle this case.
 309      */
 310     if (lastAtom + 1 < moltype.atoms.nr &&
 311         isParticleVsite[lastAtom + 1])
 312     {
 313         AtomIndexExtremes extremes = vsiteConstructRange(lastAtom + 1, moltype);
 314         if (extremes.minAtom < firstAtom)
 315         {
 316             /* Constructing atom precedes the group */
 317             return 0;
 318         }
 319         else if (extremes.maxAtom <= lastAtom)
 320         {
 321             /* All constructing atoms are in the group, add the vsite to the group */
 322             lastAtom++;
 323         }
 324         else if (extremes.minAtom <= lastAtom)
 325         {
 326             /* Some, but not all constructing atoms are in the group */
 327             return 0;
 328         }
 329     }
 330
 331     GMX_RELEASE_ASSERT(maxConstraintsPerAtom < numAtomsWithConstraints,
 332                        "We have checked that atoms are only constrained to atoms within the group,"
 333                        "so each atom should have fewer constraints than the number of atoms");
 334     /* Check that at least one atom is constrained to all others */
 335     if (maxConstraintsPerAtom != numAtomsWithConstraints - 1)
 336     {
 337         return 0;
 338     }
 339
 340     return lastAtom - firstAtom + 1;
 341 }
 342
 343 /*! \brief Returns a list of update groups for \p moltype */
 344 static RangePartitioning
 345 makeUpdateGroups(const gmx_moltype_t            &moltype,
 346                  gmx::ArrayRef<const t_iparams>  iparams)
 347 {
 348     RangePartitioning groups;
 349
 350     /* Update groups are not compatible with flexible constraints.
 351      * Without dynamics the flexible constraints are ignored,
 352      * but since performance for EM/NM is less critical, we do not
 353      * use update groups to keep the code here simpler.
 354      */
 355     if (hasFlexibleConstraints(moltype, iparams) ||
 356         hasIncompatibleVsites(moltype, iparams))
 357     {
 358         return groups;
 359     }
 360
 361     /* Combine all constraint ilists into a single one */
 362     InteractionList constraintsCombined = jointConstraintList(moltype);
 363     t_ilist         ilistsCombined[F_NRE];
 364     ilistsCombined[F_CONSTR].nr         = constraintsCombined.iatoms.size();
 365     ilistsCombined[F_CONSTR].iatoms     = constraintsCombined.iatoms.data();
 366     ilistsCombined[F_CONSTRNC].nr       = 0;
 367     /* We "include" flexible constraints, but none are present (checked above) */
 368     t_blocka             at2con         = make_at2con(moltype.atoms.nr,
 369                                                       ilistsCombined, iparams.data(),
 370                                                       FlexibleConstraintTreatment::Include);
 371
 372     bool satisfiesCriteria = true;
 373
 374     int  firstAtom         = 0;
 375     while (satisfiesCriteria && firstAtom < moltype.atoms.nr)
 376     {
 377         int numAtomsInGroup =
 378             detectGroup(firstAtom, moltype, at2con, constraintsCombined);
 379
 380         if (numAtomsInGroup == 0)
 381         {
 382             satisfiesCriteria = false;
 383         }
 384         else
 385         {
 386             groups.appendBlock(numAtomsInGroup);
 387         }
 388         firstAtom += numAtomsInGroup;
 389     }
 390
 391     if (!satisfiesCriteria)
 392     {
 393         /* Make groups empty, to signal not satisfying the criteria */
 394         groups.clear();
 395     }
 396
 397     done_blocka(&at2con);
 398
 399     return groups;
 400 }
 401
 402 std::vector<RangePartitioning> makeUpdateGroups(const gmx_mtop_t &mtop)
 403 {
 404     std::vector<RangePartitioning> updateGroups;
 405
 406     bool                           systemSatisfiesCriteria = true;
 407     for (const gmx_moltype_t &moltype : mtop.moltype)
 408     {
 409         updateGroups.push_back(makeUpdateGroups(moltype, mtop.ffparams.iparams));
 410
 411         if (updateGroups.back().numBlocks() == 0)
 412         {
 413             systemSatisfiesCriteria = false;
 414         }
 415     }
 416
 417     if (!systemSatisfiesCriteria)
 418     {
 419         updateGroups.clear();
 420     }
 421
 422     return updateGroups;
 423 }
 424
 425 /*! \brief Returns the maximum update group radius for \p moltype */
 426 static real
 427 computeMaxUpdateGroupRadius(const gmx_moltype_t            &moltype,
 428                             gmx::ArrayRef<const t_iparams>  iparams,
 429                             const RangePartitioning        &updateGroups)
 430 {
 431     GMX_RELEASE_ASSERT(!hasFlexibleConstraints(moltype, iparams),
 432                        "Flexible constraints are not supported here");
 433
 434     const InteractionList &settles = moltype.ilist[F_SETTLE];
 435
 436     t_blocka               at2con =
 437         make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include);
 438
 439     const int  stride = 1 + NRAL(F_CONSTR);
 440     const real half   = 0.5;
 441
 442     real       maxRadius = 0;
 443     for (int group = 0; group < updateGroups.numBlocks(); group++)
 444     {
 445         if (updateGroups.block(group).size() == 1)
 446         {
 447             /* Single atom group, radius is zero */
 448             continue;
 449         }
 450
 451         /* Find the atom maxAtom with the maximum number of constraints */
 452         int maxNumConstraints = 0;
 453         int maxAtom           = -1;
 454         for (int a : updateGroups.block(group))
 455         {
 456             int numConstraints = at2con.index[a + 1] - at2con.index[a];
 457             if (numConstraints > maxNumConstraints)
 458             {
 459                 maxNumConstraints = numConstraints;
 460                 maxAtom           = a;
 461             }
 462         }
 463         GMX_ASSERT(maxAtom >= 0 || settles.size() > 0,
 464                    "We should have at least two atoms in the group with constraints");
 465         if (maxAtom < 0)
 466         {
 467             continue;
 468         }
 469
 470         for (int i = at2con.index[maxAtom]; i < at2con.index[maxAtom + 1]; i++)
 471         {
 472             int conIndex = at2con.a[i]*stride;
 473             int iparamsIndex;
 474             if (conIndex < moltype.ilist[F_CONSTR].size())
 475             {
 476                 iparamsIndex = moltype.ilist[F_CONSTR].iatoms[conIndex];
 477             }
 478             else
 479             {
 480                 iparamsIndex = moltype.ilist[F_CONSTRNC].iatoms[conIndex - moltype.ilist[F_CONSTR].size()];
 481             }
 482             /* Here we take the maximum constraint length of the A and B
 483              * topology, which assumes lambda is between 0 and 1 for
 484              * free-energy runs.
 485              */
 486             real constraintLength = std::max(iparams[iparamsIndex].constr.dA,
 487                                              iparams[iparamsIndex].constr.dB);
 488             if (at2con.index[maxAtom + 1] == at2con.index[maxAtom] + 1)
 489             {
 490                 /* Single constraint: use the distance from the midpoint */
 491                 maxRadius = std::max(maxRadius, half*constraintLength);
 492             }
 493             else
 494             {
 495                 /* Multiple constraints: use the distance from the central
 496                  * atom. We can do better than this if we would assume
 497                  * that the angles between bonds do not stretch a lot.
 498                  */
 499                 maxRadius = std::max(maxRadius, constraintLength);
 500             }
 501         }
 502     }
 503
 504     for (int i = 0; i < settles.size(); i += 1 + NRAL(F_SETTLE))
 505     {
 506         const real dOH   = iparams[settles.iatoms[i]].settle.doh;
 507         const real dHH   = iparams[settles.iatoms[i]].settle.dhh;
 508         /* Compute distance^2 from center of geometry to O and H */
 509         const real dCO2  = (4*dOH*dOH - dHH*dHH)/9;
 510         const real dCH2  = (dOH*dOH + 2*dHH*dHH)/9;
 511         const real dCAny = std::sqrt(std::max(dCO2, dCH2));
 512         maxRadius        = std::max(maxRadius, dCAny);
 513     }
 514
 515     done_blocka(&at2con);
 516
 517     return maxRadius;
 518 }
 519
 520 real computeMaxUpdateGroupRadius(const gmx_mtop_t                       &mtop,
 521                                  gmx::ArrayRef<const RangePartitioning>  updateGroups)
 522 {
 523     if (updateGroups.empty())
 524     {
 525         return 0;
 526     }
 527
 528     GMX_RELEASE_ASSERT(static_cast<size_t>(updateGroups.size()) == mtop.moltype.size(),
 529                        "We need one update group entry per moleculetype");
 530
 531     real maxRadius = 0;
 532
 533     for (size_t moltype = 0; moltype < mtop.moltype.size(); moltype++)
 534     {
 535         maxRadius
 536             = std::max(maxRadius,
 537                        computeMaxUpdateGroupRadius(mtop.moltype[moltype],
 538                                                    mtop.ffparams.iparams,
 539                                                    updateGroups[moltype]));
 540     }
 541
 542     return maxRadius;
 543 }
 544
 545 } // namespace gmx