src/gromacs/gmxpreprocess/vsite_parm.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "vsite_parm.h"
  41
  42 #include <cmath>
  43 #include <cstdio>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/gmxpreprocess/add_par.h"
  49 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  50 #include "gromacs/gmxpreprocess/grompp_impl.h"
  51 #include "gromacs/gmxpreprocess/notset.h"
  52 #include "gromacs/gmxpreprocess/toputil.h"
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/mdtypes/md_enums.h"
  57 #include "gromacs/topology/ifunc.h"
  58 #include "gromacs/topology/topology.h"
  59 #include "gromacs/utility/arrayref.h"
  60 #include "gromacs/utility/cstringutil.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/gmxassert.h"
  63 #include "gromacs/utility/logger.h"
  64 #include "gromacs/utility/smalloc.h"
  65 #include "gromacs/utility/strconvert.h"
  66
  67 #include "hackblock.h"
  68 #include "resall.h"
  69
  70 /*! \internal \brief
  71  * Data type to store information about bonded interactions for virtual sites.
  72  */
  73 class VsiteBondedInteraction
  74 {
  75 public:
  76     //! Constructor initializes datastructure.
  77     VsiteBondedInteraction(gmx::ArrayRef<const int> atomIndex, real parameterValue) :
  78         parameterValue_(parameterValue)
  79     {
  80         GMX_RELEASE_ASSERT(atomIndex.size() <= atomIndex_.size(),
  81                            "Cannot add more atom indices than maximum number");
  82         auto atomIndexIt = atomIndex_.begin();
  83         for (const auto index : atomIndex)
  84         {
  85             *atomIndexIt++ = index;
  86         }
  87     }
  88     /*!@{*/
  89     //! Access the individual elements set for the vsite bonded interaction.
  90     const int& ai() const { return atomIndex_[0]; }
  91     const int& aj() const { return atomIndex_[1]; }
  92     const int& ak() const { return atomIndex_[2]; }
  93     const int& al() const { return atomIndex_[3]; }
  94
  95     const real& parameterValue() const { return parameterValue_; }
  96     /*!@}*/
  97
  98 private:
  99     //! The distance value for this bonded interaction.
 100     real parameterValue_;
 101     //! Array of atom indices
 102     std::array<int, 4> atomIndex_;
 103 };
 104
 105 /*! \internal \brief
 106  * Stores information about single virtual site bonded parameter.
 107  */
 108 struct VsiteBondParameter
 109 {
 110     //! Constructor initializes datastructure.
 111     VsiteBondParameter(int ftype, const InteractionOfType& vsiteInteraction) :
 112         ftype_(ftype),
 113         vsiteInteraction_(vsiteInteraction)
 114     {
 115     }
 116     //! Function type for virtual site.
 117     int ftype_;
 118     /*! \brief
 119      * Interaction type data.
 120      *
 121      * The datastructure should never be used in a case where the InteractionType
 122      * used to construct it might go out of scope before this object, as it would cause
 123      * the reference to type_ to dangle.
 124      */
 125     const InteractionOfType& vsiteInteraction_;
 126 };
 127
 128 /*! \internal \brief
 129  * Helper type for conversion of bonded parameters to virtual sites.
 130  */
 131 struct Atom2VsiteBond
 132 {
 133     //! Function type for conversion.
 134     int ftype;
 135     //! The vsite parameters in a list.
 136     std::vector<VsiteBondParameter> vSiteBondedParameters;
 137 };
 138
 139 //! Convenience type def for virtual site connections.
 140 using Atom2VsiteConnection = std::vector<int>;
 141
 142 static int vsite_bond_nrcheck(int ftype)
 143 {
 144     int nrcheck;
 145
 146     if ((interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT | IF_ATYPE)) || (ftype == F_IDIHS))
 147     {
 148         nrcheck = NRAL(ftype);
 149     }
 150     else
 151     {
 152         nrcheck = 0;
 153     }
 154
 155     return nrcheck;
 156 }
 157
 158 static void enter_bonded(int nratoms, std::vector<VsiteBondedInteraction>* bondeds, const InteractionOfType& type)
 159 {
 160     GMX_RELEASE_ASSERT(nratoms == type.atoms().ssize(), "Size of atom array must match");
 161     bondeds->emplace_back(VsiteBondedInteraction(type.atoms(), type.c0()));
 162 }
 163
 164 /*! \internal \brief
 165  * Wraps the datastructures for the different vsite bondeds.
 166  */
 167 struct AllVsiteBondedInteractions
 168 {
 169     //! Bond vsites.
 170     std::vector<VsiteBondedInteraction> bonds;
 171     //! Angle vsites.
 172     std::vector<VsiteBondedInteraction> angles;
 173     //! Dihedral vsites.
 174     std::vector<VsiteBondedInteraction> dihedrals;
 175 };
 176
 177 static AllVsiteBondedInteractions createVsiteBondedInformation(int                      nrat,
 178                                                                gmx::ArrayRef<const int> atoms,
 179                                                                gmx::ArrayRef<const Atom2VsiteBond> at2vb)
 180 {
 181     AllVsiteBondedInteractions allVsiteBondeds;
 182     for (int k = 0; k < nrat; k++)
 183     {
 184         for (auto& vsite : at2vb[atoms[k]].vSiteBondedParameters)
 185         {
 186             int                      ftype   = vsite.ftype_;
 187             const InteractionOfType& type    = vsite.vsiteInteraction_;
 188             int                      nrcheck = vsite_bond_nrcheck(ftype);
 189             /* abuse nrcheck to see if we're adding bond, angle or idih */
 190             switch (nrcheck)
 191             {
 192                 case 2: enter_bonded(nrcheck, &allVsiteBondeds.bonds, type); break;
 193                 case 3: enter_bonded(nrcheck, &allVsiteBondeds.angles, type); break;
 194                 case 4: enter_bonded(nrcheck, &allVsiteBondeds.dihedrals, type); break;
 195             }
 196         }
 197     }
 198     return allVsiteBondeds;
 199 }
 200
 201 static std::vector<Atom2VsiteBond> make_at2vsitebond(int natoms, gmx::ArrayRef<InteractionsOfType> plist)
 202 {
 203     bool* bVSI;
 204
 205     std::vector<Atom2VsiteBond> at2vb(natoms);
 206
 207     snew(bVSI, natoms);
 208     for (int ftype = 0; (ftype < F_NRE); ftype++)
 209     {
 210         if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
 211         {
 212             for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
 213             {
 214                 gmx::ArrayRef<const int> atoms = plist[ftype].interactionTypes[i].atoms();
 215                 for (int j = 0; j < NRAL(ftype); j++)
 216                 {
 217                     bVSI[atoms[j]] = TRUE;
 218                 }
 219             }
 220         }
 221     }
 222
 223     for (int ftype = 0; (ftype < F_NRE); ftype++)
 224     {
 225         int nrcheck = vsite_bond_nrcheck(ftype);
 226         if (nrcheck > 0)
 227         {
 228             for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
 229             {
 230                 gmx::ArrayRef<const int> aa = plist[ftype].interactionTypes[i].atoms();
 231                 for (int j = 0; j < nrcheck; j++)
 232                 {
 233                     if (bVSI[aa[j]])
 234                     {
 235                         at2vb[aa[j]].vSiteBondedParameters.emplace_back(
 236                                 ftype, plist[ftype].interactionTypes[i]);
 237                     }
 238                 }
 239             }
 240         }
 241     }
 242     sfree(bVSI);
 243
 244     return at2vb;
 245 }
 246
 247 static std::vector<Atom2VsiteConnection> make_at2vsitecon(int natoms, gmx::ArrayRef<InteractionsOfType> plist)
 248 {
 249     std::vector<bool>                 bVSI(natoms);
 250     std::vector<Atom2VsiteConnection> at2vc(natoms);
 251
 252     for (int ftype = 0; (ftype < F_NRE); ftype++)
 253     {
 254         if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
 255         {
 256             for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
 257             {
 258                 gmx::ArrayRef<const int> atoms = plist[ftype].interactionTypes[i].atoms();
 259                 for (int j = 0; j < NRAL(ftype); j++)
 260                 {
 261                     bVSI[atoms[j]] = TRUE;
 262                 }
 263             }
 264         }
 265     }
 266
 267     for (int ftype = 0; (ftype < F_NRE); ftype++)
 268     {
 269         if (interaction_function[ftype].flags & IF_CONSTRAINT)
 270         {
 271             for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
 272             {
 273                 int ai = plist[ftype].interactionTypes[i].ai();
 274                 int aj = plist[ftype].interactionTypes[i].aj();
 275                 if (bVSI[ai] && bVSI[aj])
 276                 {
 277                     /* Store forward direction */
 278                     at2vc[ai].emplace_back(aj);
 279                     /* Store backward direction */
 280                     at2vc[aj].emplace_back(ai);
 281                 }
 282             }
 283         }
 284     }
 285     return at2vc;
 286 }
 287
 288 /* for debug */
 289 static void print_bad(FILE*                                       fp,
 290                       gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 291                       gmx::ArrayRef<const VsiteBondedInteraction> angles,
 292                       gmx::ArrayRef<const VsiteBondedInteraction> idihs)
 293 {
 294     if (!bonds.empty())
 295     {
 296         fprintf(fp, "bonds:");
 297         for (const auto& bond : bonds)
 298         {
 299             fprintf(fp, " %d-%d (%g)", bond.ai() + 1, bond.aj() + 1, bond.parameterValue());
 300         }
 301         fprintf(fp, "\n");
 302     }
 303     if (!angles.empty())
 304     {
 305         fprintf(fp, "angles:");
 306         for (const auto& angle : angles)
 307         {
 308             fprintf(fp, " %d-%d-%d (%g)", angle.ai() + 1, angle.aj() + 1, angle.ak() + 1,
 309                     angle.parameterValue());
 310         }
 311         fprintf(fp, "\n");
 312     }
 313     if (!idihs.empty())
 314     {
 315         fprintf(fp, "idihs:");
 316         for (const auto& idih : idihs)
 317         {
 318             fprintf(fp, " %d-%d-%d-%d (%g)", idih.ai() + 1, idih.aj() + 1, idih.ak() + 1,
 319                     idih.al() + 1, idih.parameterValue());
 320         }
 321         fprintf(fp, "\n");
 322     }
 323 }
 324
 325 static void printInteractionOfType(FILE* fp, int ftype, int i, const InteractionOfType& type)
 326 {
 327     static int pass       = 0;
 328     static int prev_ftype = NOTSET;
 329     static int prev_i     = NOTSET;
 330
 331     if ((ftype != prev_ftype) || (i != prev_i))
 332     {
 333         pass       = 0;
 334         prev_ftype = ftype;
 335         prev_i     = i;
 336     }
 337     fprintf(fp, "(%d) plist[%s].param[%d]", pass, interaction_function[ftype].name, i);
 338     gmx::ArrayRef<const real> forceParam = type.forceParam();
 339     for (int j = 0; j < NRFP(ftype); j++)
 340     {
 341         fprintf(fp, ".c[%d]=%g ", j, forceParam[j]);
 342     }
 343     fprintf(fp, "\n");
 344     pass++;
 345 }
 346
 347 static real get_bond_length(gmx::ArrayRef<const VsiteBondedInteraction> bonds, t_iatom ai, t_iatom aj)
 348 {
 349     real bondlen;
 350
 351     bondlen = NOTSET;
 352     for (const auto& bond : bonds)
 353     {
 354         /* check both ways */
 355         if (((ai == bond.ai()) && (aj == bond.aj())) || ((ai == bond.aj()) && (aj == bond.ai())))
 356         {
 357             bondlen = bond.parameterValue(); /* note: parameterValue might be NOTSET */
 358             break;
 359         }
 360     }
 361     return bondlen;
 362 }
 363
 364 static real get_angle(gmx::ArrayRef<const VsiteBondedInteraction> angles, t_iatom ai, t_iatom aj, t_iatom ak)
 365 {
 366     real angle;
 367
 368     angle = NOTSET;
 369     for (const auto& ang : angles)
 370     {
 371         /* check both ways */
 372         if (((ai == ang.ai()) && (aj == ang.aj()) && (ak == ang.ak()))
 373             || ((ai == ang.ak()) && (aj == ang.aj()) && (ak == ang.ai())))
 374         {
 375             angle = DEG2RAD * ang.parameterValue();
 376             break;
 377         }
 378     }
 379     return angle;
 380 }
 381
 382 static const char* get_atomtype_name_AB(t_atom* atom, PreprocessingAtomTypes* atypes)
 383 {
 384     const char* name = atypes->atomNameFromAtomType(atom->type);
 385
 386     /* When using the decoupling option, atom types are changed
 387      * to decoupled for the non-bonded interactions, but the virtual
 388      * sites constructions should be based on the "normal" interactions.
 389      * So we return the state B atom type name if the state A atom
 390      * type is the decoupled one. We should actually check for the atom
 391      * type number, but that's not passed here. So we check for
 392      * the decoupled atom type name. This should not cause trouble
 393      * as this code is only used for topologies with v-sites without
 394      * parameters generated by pdb2gmx.
 395      */
 396     if (strcmp(name, "decoupled") == 0)
 397     {
 398         name = atypes->atomNameFromAtomType(atom->typeB);
 399     }
 400
 401     return name;
 402 }
 403
 404 static bool calc_vsite3_param(PreprocessingAtomTypes*                     atypes,
 405                               InteractionOfType*                          vsite,
 406                               t_atoms*                                    at,
 407                               gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 408                               gmx::ArrayRef<const VsiteBondedInteraction> angles)
 409 {
 410     /* i = virtual site          |    ,k
 411      * j = 1st bonded heavy atom | i-j
 412      * k,l = 2nd bonded atoms    |    `l
 413      */
 414
 415     bool bXH3, bError;
 416     real bjk, bjl, a = -1, b = -1;
 417     /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
 418      * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 419     bXH3 = ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MNH", 3))
 420             && (gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->al()], atypes),
 421                                           "MNH", 3)))
 422            || ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes),
 423                                           "MCH3", 4))
 424                && (gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->al()], atypes),
 425                                              "MCH3", 4)));
 426
 427     bjk    = get_bond_length(bonds, vsite->aj(), vsite->ak());
 428     bjl    = get_bond_length(bonds, vsite->aj(), vsite->al());
 429     bError = (bjk == NOTSET) || (bjl == NOTSET);
 430     if (bXH3)
 431     {
 432         /* now we get some XH2/XH3 group specific construction */
 433         /* note: we call the heavy atom 'C' and the X atom 'N' */
 434         real bMM, bCM, bCN, bNH, aCNH, dH, rH, dM, rM;
 435         int  aN;
 436
 437         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 438         bError = bError || (bjk != bjl);
 439
 440         /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
 441         aN = std::max(vsite->ak(), vsite->al()) + 1;
 442
 443         /* get common bonds */
 444         bMM    = get_bond_length(bonds, vsite->ak(), vsite->al());
 445         bCM    = bjk;
 446         bCN    = get_bond_length(bonds, vsite->aj(), aN);
 447         bError = bError || (bMM == NOTSET) || (bCN == NOTSET);
 448
 449         /* calculate common things */
 450         rM = 0.5 * bMM;
 451         dM = std::sqrt(gmx::square(bCM) - gmx::square(rM));
 452
 453         /* are we dealing with the X atom? */
 454         if (vsite->ai() == aN)
 455         {
 456             /* this is trivial */
 457             a = b = 0.5 * bCN / dM;
 458         }
 459         else
 460         {
 461             /* get other bondlengths and angles: */
 462             bNH    = get_bond_length(bonds, aN, vsite->ai());
 463             aCNH   = get_angle(angles, vsite->aj(), aN, vsite->ai());
 464             bError = bError || (bNH == NOTSET) || (aCNH == NOTSET);
 465
 466             /* calculate */
 467             dH = bCN - bNH * std::cos(aCNH);
 468             rH = bNH * std::sin(aCNH);
 469
 470             a = 0.5 * (dH / dM + rH / rM);
 471             b = 0.5 * (dH / dM - rH / rM);
 472         }
 473     }
 474     else
 475     {
 476         gmx_fatal(FARGS,
 477                   "calc_vsite3_param not implemented for the general case "
 478                   "(atom %d)",
 479                   vsite->ai() + 1);
 480     }
 481     vsite->setForceParameter(0, a);
 482     vsite->setForceParameter(1, b);
 483
 484     return bError;
 485 }
 486
 487 static bool calc_vsite3fd_param(InteractionOfType*                          vsite,
 488                                 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 489                                 gmx::ArrayRef<const VsiteBondedInteraction> angles)
 490 {
 491     /* i = virtual site          |    ,k
 492      * j = 1st bonded heavy atom | i-j
 493      * k,l = 2nd bonded atoms    |    `l
 494      */
 495
 496     bool bError;
 497     real bij, bjk, bjl, aijk, aijl, rk, rl;
 498
 499     bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 500     bjk  = get_bond_length(bonds, vsite->aj(), vsite->ak());
 501     bjl  = get_bond_length(bonds, vsite->aj(), vsite->al());
 502     aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 503     aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 504     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET);
 505
 506     rk = bjk * std::sin(aijk);
 507     rl = bjl * std::sin(aijl);
 508     vsite->setForceParameter(0, rk / (rk + rl));
 509     vsite->setForceParameter(1, -bij);
 510
 511     return bError;
 512 }
 513
 514 static bool calc_vsite3fad_param(InteractionOfType*                          vsite,
 515                                  gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 516                                  gmx::ArrayRef<const VsiteBondedInteraction> angles)
 517 {
 518     /* i = virtual site          |
 519      * j = 1st bonded heavy atom | i-j
 520      * k = 2nd bonded heavy atom |    `k-l
 521      * l = 3d bonded heavy atom  |
 522      */
 523
 524     bool bSwapParity, bError;
 525     real bij, aijk;
 526
 527     bSwapParity = (vsite->c1() == -1);
 528
 529     bij    = get_bond_length(bonds, vsite->ai(), vsite->aj());
 530     aijk   = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 531     bError = (bij == NOTSET) || (aijk == NOTSET);
 532
 533     vsite->setForceParameter(1, bij);
 534     vsite->setForceParameter(0, RAD2DEG * aijk);
 535
 536     if (bSwapParity)
 537     {
 538         vsite->setForceParameter(0, 360 - vsite->c0());
 539     }
 540
 541     return bError;
 542 }
 543
 544 static bool calc_vsite3out_param(PreprocessingAtomTypes*                     atypes,
 545                                  InteractionOfType*                          vsite,
 546                                  t_atoms*                                    at,
 547                                  gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 548                                  gmx::ArrayRef<const VsiteBondedInteraction> angles)
 549 {
 550     /* i = virtual site          |    ,k
 551      * j = 1st bonded heavy atom | i-j
 552      * k,l = 2nd bonded atoms    |    `l
 553      * NOTE: i is out of the j-k-l plane!
 554      */
 555
 556     bool bXH3, bError, bSwapParity;
 557     real bij, bjk, bjl, aijk, aijl, akjl, pijk, pijl, a, b, c;
 558
 559     /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
 560      * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 561     bXH3 = ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MNH", 3))
 562             && (gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->al()], atypes),
 563                                           "MNH", 3)))
 564            || ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes),
 565                                           "MCH3", 4))
 566                && (gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->al()], atypes),
 567                                              "MCH3", 4)));
 568
 569     /* check if construction parity must be swapped */
 570     bSwapParity = (vsite->c1() == -1);
 571
 572     bjk    = get_bond_length(bonds, vsite->aj(), vsite->ak());
 573     bjl    = get_bond_length(bonds, vsite->aj(), vsite->al());
 574     bError = (bjk == NOTSET) || (bjl == NOTSET);
 575     if (bXH3)
 576     {
 577         /* now we get some XH3 group specific construction */
 578         /* note: we call the heavy atom 'C' and the X atom 'N' */
 579         real bMM, bCM, bCN, bNH, aCNH, dH, rH, rHx, rHy, dM, rM;
 580         int  aN;
 581
 582         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 583         bError = bError || (bjk != bjl);
 584
 585         /* the X atom (C or N) in the XH3 group is the first after the masses: */
 586         aN = std::max(vsite->ak(), vsite->al()) + 1;
 587
 588         /* get all bondlengths and angles: */
 589         bMM    = get_bond_length(bonds, vsite->ak(), vsite->al());
 590         bCM    = bjk;
 591         bCN    = get_bond_length(bonds, vsite->aj(), aN);
 592         bNH    = get_bond_length(bonds, aN, vsite->ai());
 593         aCNH   = get_angle(angles, vsite->aj(), aN, vsite->ai());
 594         bError = bError || (bMM == NOTSET) || (bCN == NOTSET) || (bNH == NOTSET) || (aCNH == NOTSET);
 595
 596         /* calculate */
 597         dH = bCN - bNH * std::cos(aCNH);
 598         rH = bNH * std::sin(aCNH);
 599         /* we assume the H's are symmetrically distributed */
 600         rHx = rH * std::cos(DEG2RAD * 30);
 601         rHy = rH * std::sin(DEG2RAD * 30);
 602         rM  = 0.5 * bMM;
 603         dM  = std::sqrt(gmx::square(bCM) - gmx::square(rM));
 604         a   = 0.5 * ((dH / dM) - (rHy / rM));
 605         b   = 0.5 * ((dH / dM) + (rHy / rM));
 606         c   = rHx / (2 * dM * rM);
 607     }
 608     else
 609     {
 610         /* this is the general construction */
 611
 612         bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 613         aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 614         aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 615         akjl = get_angle(angles, vsite->ak(), vsite->aj(), vsite->al());
 616         bError = bError || (bij == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET) || (akjl == NOTSET);
 617
 618         pijk = std::cos(aijk) * bij;
 619         pijl = std::cos(aijl) * bij;
 620         a = (pijk + (pijk * std::cos(akjl) - pijl) * std::cos(akjl) / gmx::square(std::sin(akjl))) / bjk;
 621         b = (pijl + (pijl * std::cos(akjl) - pijk) * std::cos(akjl) / gmx::square(std::sin(akjl))) / bjl;
 622         c = -std::sqrt(gmx::square(bij)
 623                        - (gmx::square(pijk) - 2 * pijk * pijl * std::cos(akjl) + gmx::square(pijl))
 624                                  / gmx::square(std::sin(akjl)))
 625             / (bjk * bjl * std::sin(akjl));
 626     }
 627
 628     vsite->setForceParameter(0, a);
 629     vsite->setForceParameter(1, b);
 630     if (bSwapParity)
 631     {
 632         vsite->setForceParameter(2, -c);
 633     }
 634     else
 635     {
 636         vsite->setForceParameter(2, c);
 637     }
 638     return bError;
 639 }
 640
 641 static bool calc_vsite4fd_param(InteractionOfType*                          vsite,
 642                                 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 643                                 gmx::ArrayRef<const VsiteBondedInteraction> angles,
 644                                 const gmx::MDLogger&                        logger)
 645 {
 646     /* i = virtual site          |    ,k
 647      * j = 1st bonded heavy atom | i-j-m
 648      * k,l,m = 2nd bonded atoms  |    `l
 649      */
 650
 651     bool bError;
 652     real bij, bjk, bjl, bjm, aijk, aijl, aijm, akjm, akjl;
 653     real pk, pl, pm, cosakl, cosakm, sinakl, sinakm, cl, cm;
 654
 655     bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 656     bjk  = get_bond_length(bonds, vsite->aj(), vsite->ak());
 657     bjl  = get_bond_length(bonds, vsite->aj(), vsite->al());
 658     bjm  = get_bond_length(bonds, vsite->aj(), vsite->am());
 659     aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 660     aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 661     aijm = get_angle(angles, vsite->ai(), vsite->aj(), vsite->am());
 662     akjm = get_angle(angles, vsite->ak(), vsite->aj(), vsite->am());
 663     akjl = get_angle(angles, vsite->ak(), vsite->aj(), vsite->al());
 664     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET) || (aijk == NOTSET)
 665              || (aijl == NOTSET) || (aijm == NOTSET) || (akjm == NOTSET) || (akjl == NOTSET);
 666
 667     if (!bError)
 668     {
 669         pk     = bjk * std::sin(aijk);
 670         pl     = bjl * std::sin(aijl);
 671         pm     = bjm * std::sin(aijm);
 672         cosakl = (std::cos(akjl) - std::cos(aijk) * std::cos(aijl)) / (std::sin(aijk) * std::sin(aijl));
 673         cosakm = (std::cos(akjm) - std::cos(aijk) * std::cos(aijm)) / (std::sin(aijk) * std::sin(aijm));
 674         if (cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1)
 675         {
 676             GMX_LOG(logger.warning)
 677                     .asParagraph()
 678                     .appendTextFormatted(
 679                             "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f",
 680                             vsite->ai() + 1, RAD2DEG * aijk, RAD2DEG * aijl, RAD2DEG * aijm);
 681             gmx_fatal(FARGS,
 682                       "invalid construction in calc_vsite4fd for atom %d: "
 683                       "cosakl=%f, cosakm=%f\n",
 684                       vsite->ai() + 1, cosakl, cosakm);
 685         }
 686         sinakl = std::sqrt(1 - gmx::square(cosakl));
 687         sinakm = std::sqrt(1 - gmx::square(cosakm));
 688
 689         /* note: there is a '+' because of the way the sines are calculated */
 690         cl = -pk / (pl * cosakl - pk + pl * sinakl * (pm * cosakm - pk) / (pm * sinakm));
 691         cm = -pk / (pm * cosakm - pk + pm * sinakm * (pl * cosakl - pk) / (pl * sinakl));
 692
 693         vsite->setForceParameter(0, cl);
 694         vsite->setForceParameter(1, cm);
 695         vsite->setForceParameter(2, -bij);
 696     }
 697
 698     return bError;
 699 }
 700
 701
 702 static bool calc_vsite4fdn_param(InteractionOfType*                          vsite,
 703                                  gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 704                                  gmx::ArrayRef<const VsiteBondedInteraction> angles,
 705                                  const gmx::MDLogger&                        logger)
 706 {
 707     /* i = virtual site          |    ,k
 708      * j = 1st bonded heavy atom | i-j-m
 709      * k,l,m = 2nd bonded atoms  |    `l
 710      */
 711
 712     bool bError;
 713     real bij, bjk, bjl, bjm, aijk, aijl, aijm;
 714     real pk, pl, pm, a, b;
 715
 716     bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 717     bjk  = get_bond_length(bonds, vsite->aj(), vsite->ak());
 718     bjl  = get_bond_length(bonds, vsite->aj(), vsite->al());
 719     bjm  = get_bond_length(bonds, vsite->aj(), vsite->am());
 720     aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 721     aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 722     aijm = get_angle(angles, vsite->ai(), vsite->aj(), vsite->am());
 723
 724     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET)
 725              || (aijk == NOTSET) || (aijl == NOTSET) || (aijm == NOTSET);
 726
 727     if (!bError)
 728     {
 729
 730         /* Calculate component of bond j-k along the direction i-j */
 731         pk = -bjk * std::cos(aijk);
 732
 733         /* Calculate component of bond j-l along the direction i-j */
 734         pl = -bjl * std::cos(aijl);
 735
 736         /* Calculate component of bond j-m along the direction i-j */
 737         pm = -bjm * std::cos(aijm);
 738
 739         if (fabs(pl) < 1000 * GMX_REAL_MIN || fabs(pm) < 1000 * GMX_REAL_MIN)
 740         {
 741             GMX_LOG(logger.warning)
 742                     .asParagraph()
 743                     .appendTextFormatted(
 744                             "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f",
 745                             vsite->ai() + 1, RAD2DEG * aijk, RAD2DEG * aijl, RAD2DEG * aijm);
 746             gmx_fatal(FARGS,
 747                       "invalid construction in calc_vsite4fdn for atom %d: "
 748                       "pl=%f, pm=%f\n",
 749                       vsite->ai() + 1, pl, pm);
 750         }
 751
 752         a = pk / pl;
 753         b = pk / pm;
 754
 755         vsite->setForceParameter(0, a);
 756         vsite->setForceParameter(1, b);
 757         vsite->setForceParameter(2, bij);
 758     }
 759
 760     return bError;
 761 }
 762
 763
 764 int set_vsites(bool                              bVerbose,
 765                t_atoms*                          atoms,
 766                PreprocessingAtomTypes*           atypes,
 767                gmx::ArrayRef<InteractionsOfType> plist,
 768                const gmx::MDLogger&              logger)
 769 {
 770     int  ftype;
 771     int  nvsite, nrset;
 772     bool bFirst, bERROR;
 773
 774     bFirst = TRUE;
 775     nvsite = 0;
 776
 777     /* Make a reverse list to avoid ninteractions^2 operations */
 778     std::vector<Atom2VsiteBond> at2vb = make_at2vsitebond(atoms->nr, plist);
 779
 780     for (ftype = 0; (ftype < F_NRE); ftype++)
 781     {
 782         if (interaction_function[ftype].flags & IF_VSITE)
 783         {
 784             nvsite += plist[ftype].size();
 785
 786             if (ftype == F_VSITEN)
 787             {
 788                 /* We don't calculate parameters for VSITEN */
 789                 continue;
 790             }
 791
 792             nrset = 0;
 793             int i = 0;
 794             for (auto& param : plist[ftype].interactionTypes)
 795             {
 796                 /* check if all parameters are set */
 797                 bool                      bSet       = true;
 798                 gmx::ArrayRef<const real> forceParam = param.forceParam();
 799                 for (int j = 0; (j < NRFP(ftype)) && bSet; j++)
 800                 {
 801                     bSet = forceParam[j] != NOTSET;
 802                 }
 803
 804                 if (debug)
 805                 {
 806                     fprintf(debug, "bSet=%s ", gmx::boolToString(bSet));
 807                     printInteractionOfType(debug, ftype, i, plist[ftype].interactionTypes[i]);
 808                 }
 809                 if (!bSet)
 810                 {
 811                     if (bVerbose && bFirst)
 812                     {
 813                         GMX_LOG(logger.info)
 814                                 .asParagraph()
 815                                 .appendTextFormatted("Calculating parameters for virtual sites");
 816                         bFirst = FALSE;
 817                     }
 818
 819                     nrset++;
 820                     /* now set the vsite parameters: */
 821                     AllVsiteBondedInteractions allVsiteBondeds =
 822                             createVsiteBondedInformation(NRAL(ftype), param.atoms(), at2vb);
 823                     if (debug)
 824                     {
 825                         fprintf(debug,
 826                                 "Found %zu bonds, %zu angles and %zu idihs "
 827                                 "for virtual site %d (%s)\n",
 828                                 allVsiteBondeds.bonds.size(), allVsiteBondeds.angles.size(),
 829                                 allVsiteBondeds.dihedrals.size(), param.ai() + 1,
 830                                 interaction_function[ftype].longname);
 831                         print_bad(debug, allVsiteBondeds.bonds, allVsiteBondeds.angles,
 832                                   allVsiteBondeds.dihedrals);
 833                     } /* debug */
 834                     switch (ftype)
 835                     {
 836                         case F_VSITE3:
 837                             bERROR = calc_vsite3_param(atypes, &param, atoms, allVsiteBondeds.bonds,
 838                                                        allVsiteBondeds.angles);
 839                             break;
 840                         case F_VSITE3FD:
 841                             bERROR = calc_vsite3fd_param(&param, allVsiteBondeds.bonds,
 842                                                          allVsiteBondeds.angles);
 843                             break;
 844                         case F_VSITE3FAD:
 845                             bERROR = calc_vsite3fad_param(&param, allVsiteBondeds.bonds,
 846                                                           allVsiteBondeds.angles);
 847                             break;
 848                         case F_VSITE3OUT:
 849                             bERROR = calc_vsite3out_param(atypes, &param, atoms, allVsiteBondeds.bonds,
 850                                                           allVsiteBondeds.angles);
 851                             break;
 852                         case F_VSITE4FD:
 853                             bERROR = calc_vsite4fd_param(&param, allVsiteBondeds.bonds,
 854                                                          allVsiteBondeds.angles, logger);
 855                             break;
 856                         case F_VSITE4FDN:
 857                             bERROR = calc_vsite4fdn_param(&param, allVsiteBondeds.bonds,
 858                                                           allVsiteBondeds.angles, logger);
 859                             break;
 860                         default:
 861                             gmx_fatal(FARGS,
 862                                       "Automatic parameter generation not supported "
 863                                       "for %s atom %d",
 864                                       interaction_function[ftype].longname, param.ai() + 1);
 865                             bERROR = TRUE;
 866                     } /* switch */
 867                     if (bERROR)
 868                     {
 869                         gmx_fatal(FARGS,
 870                                   "Automatic parameter generation not supported "
 871                                   "for %s atom %d for this bonding configuration",
 872                                   interaction_function[ftype].longname, param.ai() + 1);
 873                     }
 874                 } /* if bSet */
 875                 i++;
 876             }
 877         } /* if IF_VSITE */
 878     }
 879     return nvsite;
 880 }
 881
 882 void set_vsites_ptype(bool bVerbose, gmx_moltype_t* molt, const gmx::MDLogger& logger)
 883 {
 884     int ftype, i;
 885
 886     if (bVerbose)
 887     {
 888         GMX_LOG(logger.info)
 889                 .asParagraph()
 890                 .appendTextFormatted("Setting particle type to V for virtual sites");
 891     }
 892     for (ftype = 0; ftype < F_NRE; ftype++)
 893     {
 894         InteractionList* il = &molt->ilist[ftype];
 895         if (interaction_function[ftype].flags & IF_VSITE)
 896         {
 897             const int                nra = interaction_function[ftype].nratoms;
 898             const int                nrd = il->size();
 899             gmx::ArrayRef<const int> ia  = il->iatoms;
 900
 901             if (debug && nrd)
 902             {
 903                 GMX_LOG(logger.info)
 904                         .asParagraph()
 905                         .appendTextFormatted("doing %d %s virtual sites", (nrd / (nra + 1)),
 906                                              interaction_function[ftype].longname);
 907             }
 908
 909             for (i = 0; (i < nrd);)
 910             {
 911                 /* The virtual site */
 912                 int avsite                     = ia[i + 1];
 913                 molt->atoms.atom[avsite].ptype = eptVSite;
 914
 915                 i += nra + 1;
 916             }
 917         }
 918     }
 919 }
 920
 921 /*! \brief
 922  *  Convenience typedef for linking function type to parameter numbers.
 923  *
 924  *  The entries in this datastructure are valid if the particle participates in
 925  *  a virtual site interaction and has a valid vsite function type other than VSITEN.
 926  *  \todo Change to remove empty constructor when gmx::compat::optional is available.
 927  */
 928 class VsiteAtomMapping
 929 {
 930 public:
 931     //! Only construct with all information in place or nothing
 932     VsiteAtomMapping(int functionType, int interactionIndex) :
 933         functionType_(functionType),
 934         interactionIndex_(interactionIndex)
 935     {
 936     }
 937     VsiteAtomMapping() : functionType_(-1), interactionIndex_(-1) {}
 938     //! Get function type.
 939     const int& functionType() const { return functionType_; }
 940     //! Get parameter number.
 941     const int& interactionIndex() const { return interactionIndex_; }
 942
 943 private:
 944     //! Function type for the linked parameter.
 945     int functionType_;
 946     //! The linked parameter.
 947     int interactionIndex_;
 948 };
 949
 950 static void check_vsite_constraints(gmx::ArrayRef<InteractionsOfType> plist,
 951                                     int                               cftype,
 952                                     const int                         vsite_type[],
 953                                     const gmx::MDLogger&              logger)
 954 {
 955     int n = 0;
 956     for (const auto& param : plist[cftype].interactionTypes)
 957     {
 958         gmx::ArrayRef<const int> atoms = param.atoms();
 959         for (int k = 0; k < 2; k++)
 960         {
 961             int atom = atoms[k];
 962             if (vsite_type[atom] != NOTSET)
 963             {
 964                 GMX_LOG(logger.info)
 965                         .asParagraph()
 966                         .appendTextFormatted(
 967                                 "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)",
 968                                 param.ai() + 1, param.aj() + 1, atom + 1);
 969                 n++;
 970             }
 971         }
 972     }
 973     if (n)
 974     {
 975         gmx_fatal(FARGS, "There were %d virtual sites involved in constraints", n);
 976     }
 977 }
 978
 979 static void clean_vsite_bonds(gmx::ArrayRef<InteractionsOfType>     plist,
 980                               gmx::ArrayRef<const VsiteAtomMapping> pindex,
 981                               int                                   cftype,
 982                               const int                             vsite_type[],
 983                               const gmx::MDLogger&                  logger)
 984 {
 985     int                 ftype, nOut;
 986     int                 nconverted, nremoved;
 987     int                 oatom, at1, at2;
 988     bool                bKeep, bRemove, bAllFD;
 989     InteractionsOfType* ps;
 990
 991     if (cftype == F_CONNBONDS)
 992     {
 993         return;
 994     }
 995
 996     ps         = &(plist[cftype]);
 997     nconverted = 0;
 998     nremoved   = 0;
 999     nOut       = 0;
1000     for (auto bond = ps->interactionTypes.begin(); bond != ps->interactionTypes.end();)
1001     {
1002         int        vsnral      = 0;
1003         const int* first_atoms = nullptr;
1004
1005         bKeep   = false;
1006         bRemove = false;
1007         bAllFD  = true;
1008         /* check if all virtual sites are constructed from the same atoms */
1009         int                      nvsite = 0;
1010         gmx::ArrayRef<const int> atoms  = bond->atoms();
1011         for (int k = 0; (k < 2) && !bKeep && !bRemove; k++)
1012         {
1013             /* for all atoms in the bond */
1014             int atom = atoms[k];
1015             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1016             {
1017                 nvsite++;
1018                 bool bThisFD  = ((pindex[atom].functionType() == F_VSITE3FD)
1019                                 || (pindex[atom].functionType() == F_VSITE3FAD)
1020                                 || (pindex[atom].functionType() == F_VSITE4FD)
1021                                 || (pindex[atom].functionType() == F_VSITE4FDN));
1022                 bool bThisOUT = ((pindex[atom].functionType() == F_VSITE3OUT)
1023                                  && ((interaction_function[cftype].flags & IF_CONSTRAINT) != 0U));
1024                 bAllFD        = bAllFD && bThisFD;
1025                 if (bThisFD || bThisOUT)
1026                 {
1027                     oatom = atoms[1 - k]; /* the other atom */
1028                     if (vsite_type[oatom] == NOTSET
1029                         && oatom
1030                                    == plist[pindex[atom].functionType()]
1031                                               .interactionTypes[pindex[atom].interactionIndex()]
1032                                               .aj())
1033                     {
1034                         /* if the other atom isn't a vsite, and it is AI */
1035                         bRemove = true;
1036                         if (bThisOUT)
1037                         {
1038                             nOut++;
1039                         }
1040                     }
1041                 }
1042                 if (!bRemove)
1043                 {
1044                     /* TODO This fragment, and corresponding logic in
1045                        clean_vsite_angles and clean_vsite_dihs should
1046                        be refactored into a common function */
1047                     if (nvsite == 1)
1048                     {
1049                         /* if this is the first vsite we encounter then
1050                            store construction atoms */
1051                         /* TODO This would be nicer to implement with
1052                            a C++ "vector view" class" with an
1053                            STL-container-like interface. */
1054                         vsnral      = NRAL(pindex[atom].functionType()) - 1;
1055                         first_atoms = plist[pindex[atom].functionType()]
1056                                               .interactionTypes[pindex[atom].interactionIndex()]
1057                                               .atoms()
1058                                               .data()
1059                                       + 1;
1060                     }
1061                     else
1062                     {
1063                         GMX_ASSERT(vsnral != 0, "nvsite > 1 must have vsnral != 0");
1064                         GMX_ASSERT(first_atoms != nullptr,
1065                                    "nvsite > 1 must have first_atoms != NULL");
1066                         /* if it is not the first then
1067                            check if this vsite is constructed from the same atoms */
1068                         if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1069                         {
1070                             for (int m = 0; (m < vsnral) && !bKeep; m++)
1071                             {
1072                                 const int* atoms;
1073
1074                                 bool bPresent = false;
1075                                 atoms         = plist[pindex[atom].functionType()]
1076                                                 .interactionTypes[pindex[atom].interactionIndex()]
1077                                                 .atoms()
1078                                                 .data()
1079                                         + 1;
1080                                 for (int n = 0; (n < vsnral) && !bPresent; n++)
1081                                 {
1082                                     if (atoms[m] == first_atoms[n])
1083                                     {
1084                                         bPresent = true;
1085                                     }
1086                                 }
1087                                 if (!bPresent)
1088                                 {
1089                                     bKeep = true;
1090                                 }
1091                             }
1092                         }
1093                         else
1094                         {
1095                             bKeep = true;
1096                         }
1097                     }
1098                 }
1099             }
1100         }
1101
1102         if (bRemove)
1103         {
1104             bKeep = false;
1105         }
1106         else
1107         {
1108             /* if we have no virtual sites in this bond, keep it */
1109             if (nvsite == 0)
1110             {
1111                 bKeep = true;
1112             }
1113
1114             /* TODO This loop and the corresponding loop in
1115                check_vsite_angles should be refactored into a common
1116                function */
1117             /* check if all non-vsite atoms are used in construction: */
1118             bool bFirstTwo = true;
1119             for (int k = 0; (k < 2) && !bKeep; k++) /* for all atoms in the bond */
1120             {
1121                 int atom = atoms[k];
1122                 if (vsite_type[atom] == NOTSET)
1123                 {
1124                     bool bUsed = false;
1125                     for (int m = 0; (m < vsnral) && !bUsed; m++)
1126                     {
1127                         GMX_ASSERT(first_atoms != nullptr,
1128                                    "If we've seen a vsite before, we know what its first atom "
1129                                    "index was");
1130
1131                         if (atom == first_atoms[m])
1132                         {
1133                             bUsed     = true;
1134                             bFirstTwo = bFirstTwo && m < 2;
1135                         }
1136                     }
1137                     if (!bUsed)
1138                     {
1139                         bKeep = true;
1140                     }
1141                 }
1142             }
1143
1144             if (!(bAllFD && bFirstTwo))
1145             {
1146                 /* Two atom bonded interactions include constraints.
1147                  * We need to remove constraints between vsite pairs that have
1148                  * a fixed distance due to being constructed from the same
1149                  * atoms, since this can be numerically unstable.
1150                  */
1151                 for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1152                 {
1153                     at1           = first_atoms[m];
1154                     at2           = first_atoms[(m + 1) % vsnral];
1155                     bool bPresent = false;
1156                     for (ftype = 0; ftype < F_NRE; ftype++)
1157                     {
1158                         if (interaction_function[ftype].flags & IF_CONSTRAINT)
1159                         {
1160                             for (auto entry = plist[ftype].interactionTypes.begin();
1161                                  (entry != plist[ftype].interactionTypes.end()) && !bPresent; entry++)
1162                             {
1163                                 /* all constraints until one matches */
1164                                 bPresent = (((entry->ai() == at1) && (entry->aj() == at2))
1165                                             || ((entry->ai() == at2) && (entry->aj() == at1)));
1166                             }
1167                         }
1168                     }
1169                     if (!bPresent)
1170                     {
1171                         bKeep = true;
1172                     }
1173                 }
1174             }
1175         }
1176
1177         if (bKeep)
1178         {
1179             ++bond;
1180         }
1181         else if (IS_CHEMBOND(cftype))
1182         {
1183             plist[F_CONNBONDS].interactionTypes.emplace_back(*bond);
1184             bond = ps->interactionTypes.erase(bond);
1185             nconverted++;
1186         }
1187         else
1188         {
1189             bond = ps->interactionTypes.erase(bond);
1190             nremoved++;
1191         }
1192     }
1193
1194     if (nremoved)
1195     {
1196         GMX_LOG(logger.info)
1197                 .asParagraph()
1198                 .appendTextFormatted("Removed   %4d %15ss with virtual sites, %zu left", nremoved,
1199                                      interaction_function[cftype].longname, ps->size());
1200     }
1201     if (nconverted)
1202     {
1203         GMX_LOG(logger.info)
1204                 .asParagraph()
1205                 .appendTextFormatted(
1206                         "Converted %4d %15ss with virtual sites to connections, %zu left",
1207                         nconverted, interaction_function[cftype].longname, ps->size());
1208     }
1209     if (nOut)
1210     {
1211         GMX_LOG(logger.info)
1212                 .asParagraph()
1213                 .appendTextFormatted(
1214                         "Warning: removed %d %ss with vsite with %s construction\n"
1215                         "         This vsite construction does not guarantee constant "
1216                         "bond-length\n"
1217                         "         If the constructions were generated by pdb2gmx ignore "
1218                         "this warning",
1219                         nOut, interaction_function[cftype].longname,
1220                         interaction_function[F_VSITE3OUT].longname);
1221     }
1222 }
1223
1224 static void clean_vsite_angles(gmx::ArrayRef<InteractionsOfType>         plist,
1225                                gmx::ArrayRef<VsiteAtomMapping>           pindex,
1226                                int                                       cftype,
1227                                const int                                 vsite_type[],
1228                                gmx::ArrayRef<const Atom2VsiteConnection> at2vc,
1229                                const gmx::MDLogger&                      logger)
1230 {
1231     int                 atom, at1, at2;
1232     InteractionsOfType* ps;
1233
1234     ps          = &(plist[cftype]);
1235     int oldSize = ps->size();
1236     for (auto angle = ps->interactionTypes.begin(); angle != ps->interactionTypes.end();)
1237     {
1238         int        vsnral      = 0;
1239         const int* first_atoms = nullptr;
1240
1241         bool bKeep    = false;
1242         bool bAll3FAD = true;
1243         /* check if all virtual sites are constructed from the same atoms */
1244         int                      nvsite = 0;
1245         gmx::ArrayRef<const int> atoms  = angle->atoms();
1246         for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1247         {
1248             int atom = atoms[k];
1249             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1250             {
1251                 nvsite++;
1252                 bAll3FAD = bAll3FAD && (pindex[atom].functionType() == F_VSITE3FAD);
1253                 if (nvsite == 1)
1254                 {
1255                     /* store construction atoms of first vsite */
1256                     vsnral      = NRAL(pindex[atom].functionType()) - 1;
1257                     first_atoms = plist[pindex[atom].functionType()]
1258                                           .interactionTypes[pindex[atom].interactionIndex()]
1259                                           .atoms()
1260                                           .data()
1261                                   + 1;
1262                 }
1263                 else
1264                 {
1265                     GMX_ASSERT(vsnral != 0,
1266                                "If we've seen a vsite before, we know how many constructing atoms "
1267                                "it had");
1268                     GMX_ASSERT(
1269                             first_atoms != nullptr,
1270                             "If we've seen a vsite before, we know what its first atom index was");
1271                     /* check if this vsite is constructed from the same atoms */
1272                     if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1273                     {
1274                         for (int m = 0; (m < vsnral) && !bKeep; m++)
1275                         {
1276                             const int* subAtoms;
1277
1278                             bool bPresent = false;
1279                             subAtoms      = plist[pindex[atom].functionType()]
1280                                                .interactionTypes[pindex[atom].interactionIndex()]
1281                                                .atoms()
1282                                                .data()
1283                                        + 1;
1284                             for (int n = 0; (n < vsnral) && !bPresent; n++)
1285                             {
1286                                 if (subAtoms[m] == first_atoms[n])
1287                                 {
1288                                     bPresent = true;
1289                                 }
1290                             }
1291                             if (!bPresent)
1292                             {
1293                                 bKeep = true;
1294                             }
1295                         }
1296                     }
1297                     else
1298                     {
1299                         bKeep = true;
1300                     }
1301                 }
1302             }
1303         }
1304
1305         /* keep all angles with no virtual sites in them or
1306            with virtual sites with more than 3 constr. atoms */
1307         if (nvsite == 0 && vsnral > 3)
1308         {
1309             bKeep = true;
1310         }
1311
1312         /* check if all non-vsite atoms are used in construction: */
1313         bool bFirstTwo = true;
1314         for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1315         {
1316             atom = atoms[k];
1317             if (vsite_type[atom] == NOTSET)
1318             {
1319                 bool bUsed = false;
1320                 for (int m = 0; (m < vsnral) && !bUsed; m++)
1321                 {
1322                     GMX_ASSERT(
1323                             first_atoms != nullptr,
1324                             "If we've seen a vsite before, we know what its first atom index was");
1325
1326                     if (atom == first_atoms[m])
1327                     {
1328                         bUsed     = true;
1329                         bFirstTwo = bFirstTwo && m < 2;
1330                     }
1331                 }
1332                 if (!bUsed)
1333                 {
1334                     bKeep = true;
1335                 }
1336             }
1337         }
1338
1339         if (!(bAll3FAD && bFirstTwo))
1340         {
1341             /* check if all constructing atoms are constrained together */
1342             for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1343             {
1344                 at1           = first_atoms[m];
1345                 at2           = first_atoms[(m + 1) % vsnral];
1346                 bool bPresent = false;
1347                 auto found    = std::find(at2vc[at1].begin(), at2vc[at1].end(), at2);
1348                 if (found != at2vc[at1].end())
1349                 {
1350                     bPresent = true;
1351                 }
1352                 if (!bPresent)
1353                 {
1354                     bKeep = true;
1355                 }
1356             }
1357         }
1358
1359         if (bKeep)
1360         {
1361             ++angle;
1362         }
1363         else
1364         {
1365             angle = ps->interactionTypes.erase(angle);
1366         }
1367     }
1368
1369     if (oldSize != gmx::ssize(*ps))
1370     {
1371         GMX_LOG(logger.info)
1372                 .asParagraph()
1373                 .appendTextFormatted("Removed   %4zu %15ss with virtual sites, %zu left",
1374                                      oldSize - ps->size(), interaction_function[cftype].longname,
1375                                      ps->size());
1376     }
1377 }
1378
1379 static void clean_vsite_dihs(gmx::ArrayRef<InteractionsOfType>     plist,
1380                              gmx::ArrayRef<const VsiteAtomMapping> pindex,
1381                              int                                   cftype,
1382                              const int                             vsite_type[],
1383                              const gmx::MDLogger&                  logger)
1384 {
1385     InteractionsOfType* ps;
1386
1387     ps = &(plist[cftype]);
1388
1389     int oldSize = ps->size();
1390     for (auto dih = ps->interactionTypes.begin(); dih != ps->interactionTypes.end();)
1391     {
1392         int        vsnral      = 0;
1393         const int* first_atoms = nullptr;
1394         int        atom;
1395
1396         gmx::ArrayRef<const int> atoms = dih->atoms();
1397         bool                     bKeep = false;
1398         /* check if all virtual sites are constructed from the same atoms */
1399         int nvsite = 0;
1400         for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1401         {
1402             atom = atoms[k];
1403             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1404             {
1405                 if (nvsite == 0)
1406                 {
1407                     /* store construction atoms of first vsite */
1408                     vsnral      = NRAL(pindex[atom].functionType()) - 1;
1409                     first_atoms = plist[pindex[atom].functionType()]
1410                                           .interactionTypes[pindex[atom].interactionIndex()]
1411                                           .atoms()
1412                                           .data()
1413                                   + 1;
1414                 }
1415                 else
1416                 {
1417                     GMX_ASSERT(vsnral != 0,
1418                                "If we've seen a vsite before, we know how many constructing atoms "
1419                                "it had");
1420                     GMX_ASSERT(
1421                             first_atoms != nullptr,
1422                             "If we've seen a vsite before, we know what its first atom index was");
1423                     /* check if this vsite is constructed from the same atoms */
1424                     if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1425                     {
1426                         for (int m = 0; (m < vsnral) && !bKeep; m++)
1427                         {
1428                             const int* subAtoms;
1429
1430                             bool bPresent = false;
1431                             subAtoms      = plist[pindex[atom].functionType()]
1432                                                .interactionTypes[pindex[atom].interactionIndex()]
1433                                                .atoms()
1434                                                .data()
1435                                        + 1;
1436                             for (int n = 0; (n < vsnral) && !bPresent; n++)
1437                             {
1438                                 if (subAtoms[m] == first_atoms[n])
1439                                 {
1440                                     bPresent = true;
1441                                 }
1442                             }
1443                             if (!bPresent)
1444                             {
1445                                 bKeep = true;
1446                             }
1447                         }
1448                     }
1449                 }
1450                 /* TODO clean_site_bonds and _angles do this increment
1451                    at the top of the loop. Refactor this for
1452                    consistency */
1453                 nvsite++;
1454             }
1455         }
1456
1457         /* keep all dihedrals with no virtual sites in them */
1458         if (nvsite == 0)
1459         {
1460             bKeep = true;
1461         }
1462
1463         /* check if all atoms in dihedral are either virtual sites, or used in
1464            construction of virtual sites. If so, keep it, if not throw away: */
1465         for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1466         {
1467             GMX_ASSERT(vsnral != 0,
1468                        "If we've seen a vsite before, we know how many constructing atoms it had");
1469             GMX_ASSERT(first_atoms != nullptr,
1470                        "If we've seen a vsite before, we know what its first atom index was");
1471             atom = atoms[k];
1472             if (vsite_type[atom] == NOTSET)
1473             {
1474                 /* vsnral will be set here, we don't get here with nvsite==0 */
1475                 bool bUsed = false;
1476                 for (int m = 0; (m < vsnral) && !bUsed; m++)
1477                 {
1478                     if (atom == first_atoms[m])
1479                     {
1480                         bUsed = true;
1481                     }
1482                 }
1483                 if (!bUsed)
1484                 {
1485                     bKeep = true;
1486                 }
1487             }
1488         }
1489
1490         if (bKeep)
1491         {
1492             ++dih;
1493         }
1494         else
1495         {
1496             dih = ps->interactionTypes.erase(dih);
1497         }
1498     }
1499
1500     if (oldSize != gmx::ssize(*ps))
1501     {
1502         GMX_LOG(logger.info)
1503                 .asParagraph()
1504                 .appendTextFormatted("Removed   %4zu %15ss with virtual sites, %zu left",
1505                                      oldSize - ps->size(), interaction_function[cftype].longname,
1506                                      ps->size());
1507     }
1508 }
1509
1510 // TODO use gmx::compat::optional for pindex.
1511 void clean_vsite_bondeds(gmx::ArrayRef<InteractionsOfType> plist,
1512                          int                               natoms,
1513                          bool                              bRmVSiteBds,
1514                          const gmx::MDLogger&              logger)
1515 {
1516     int                               nvsite, vsite;
1517     int*                              vsite_type;
1518     std::vector<VsiteAtomMapping>     pindex;
1519     std::vector<Atom2VsiteConnection> at2vc;
1520
1521     /* make vsite_type array */
1522     snew(vsite_type, natoms);
1523     for (int i = 0; i < natoms; i++)
1524     {
1525         vsite_type[i] = NOTSET;
1526     }
1527     nvsite = 0;
1528     for (int ftype = 0; ftype < F_NRE; ftype++)
1529     {
1530         if (interaction_function[ftype].flags & IF_VSITE)
1531         {
1532             nvsite += plist[ftype].size();
1533             int i = 0;
1534             while (i < gmx::ssize(plist[ftype]))
1535             {
1536                 vsite = plist[ftype].interactionTypes[i].ai();
1537                 if (vsite_type[vsite] == NOTSET)
1538                 {
1539                     vsite_type[vsite] = ftype;
1540                 }
1541                 else
1542                 {
1543                     gmx_fatal(FARGS, "multiple vsite constructions for atom %d", vsite + 1);
1544                 }
1545                 if (ftype == F_VSITEN)
1546                 {
1547                     while (i < gmx::ssize(plist[ftype]) && plist[ftype].interactionTypes[i].ai() == vsite)
1548                     {
1549                         i++;
1550                     }
1551                 }
1552                 else
1553                 {
1554                     i++;
1555                 }
1556             }
1557         }
1558     }
1559
1560     /* the rest only if we have virtual sites: */
1561     if (nvsite)
1562     {
1563         GMX_LOG(logger.info)
1564                 .asParagraph()
1565                 .appendTextFormatted("Cleaning up constraints %swith virtual sites",
1566                                      bRmVSiteBds ? "and constant bonded interactions " : "");
1567
1568         /* Make a reverse list to avoid ninteractions^2 operations */
1569         at2vc = make_at2vsitecon(natoms, plist);
1570
1571         pindex.resize(natoms);
1572         for (int ftype = 0; ftype < F_NRE; ftype++)
1573         {
1574             /* Here we skip VSITEN. In neary all practical use cases this
1575              * is not an issue, since VSITEN is intended for constructing
1576              * additional interaction sites, not for replacing normal atoms
1577              * with bonded interactions. Thus we do not expect constant
1578              * bonded interactions. If VSITEN does get used with constant
1579              * bonded interactions, these are not removed which only leads
1580              * to very minor extra computation and constant energy.
1581              * The only problematic case is onstraints between VSITEN
1582              * constructions with fixed distance (which is anyhow useless).
1583              * This will generate a fatal error in check_vsite_constraints.
1584              */
1585             if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
1586             {
1587                 for (gmx::index interactionIndex = 0; interactionIndex < gmx::ssize(plist[ftype]);
1588                      interactionIndex++)
1589                 {
1590                     int k     = plist[ftype].interactionTypes[interactionIndex].ai();
1591                     pindex[k] = VsiteAtomMapping(ftype, interactionIndex);
1592                 }
1593             }
1594         }
1595
1596         /* remove interactions that include virtual sites */
1597         for (int ftype = 0; ftype < F_NRE; ftype++)
1598         {
1599             if (((interaction_function[ftype].flags & IF_BOND) && bRmVSiteBds)
1600                 || (interaction_function[ftype].flags & IF_CONSTRAINT))
1601             {
1602                 if (interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT))
1603                 {
1604                     clean_vsite_bonds(plist, pindex, ftype, vsite_type, logger);
1605                 }
1606                 else if (interaction_function[ftype].flags & IF_ATYPE)
1607                 {
1608                     clean_vsite_angles(plist, pindex, ftype, vsite_type, at2vc, logger);
1609                 }
1610                 else if ((ftype == F_PDIHS) || (ftype == F_IDIHS))
1611                 {
1612                     clean_vsite_dihs(plist, pindex, ftype, vsite_type, logger);
1613                 }
1614             }
1615         }
1616         /* check that no remaining constraints include virtual sites */
1617         for (int ftype = 0; ftype < F_NRE; ftype++)
1618         {
1619             if (interaction_function[ftype].flags & IF_CONSTRAINT)
1620             {
1621                 check_vsite_constraints(plist, ftype, vsite_type, logger);
1622             }
1623         }
1624     }
1625     sfree(vsite_type);
1626 }