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,2021, 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, angle.parameterValue());
 309         }
 310         fprintf(fp, "\n");
 311     }
 312     if (!idihs.empty())
 313     {
 314         fprintf(fp, "idihs:");
 315         for (const auto& idih : idihs)
 316         {
 317             fprintf(fp,
 318                     " %d-%d-%d-%d (%g)",
 319                     idih.ai() + 1,
 320                     idih.aj() + 1,
 321                     idih.ak() + 1,
 322                     idih.al() + 1,
 323                     idih.parameterValue());
 324         }
 325         fprintf(fp, "\n");
 326     }
 327 }
 328
 329 static void printInteractionOfType(FILE* fp, int ftype, int i, const InteractionOfType& type)
 330 {
 331     static int pass       = 0;
 332     static int prev_ftype = NOTSET;
 333     static int prev_i     = NOTSET;
 334
 335     if ((ftype != prev_ftype) || (i != prev_i))
 336     {
 337         pass       = 0;
 338         prev_ftype = ftype;
 339         prev_i     = i;
 340     }
 341     fprintf(fp, "(%d) plist[%s].param[%d]", pass, interaction_function[ftype].name, i);
 342     gmx::ArrayRef<const real> forceParam = type.forceParam();
 343     for (int j = 0; j < NRFP(ftype); j++)
 344     {
 345         fprintf(fp, ".c[%d]=%g ", j, forceParam[j]);
 346     }
 347     fprintf(fp, "\n");
 348     pass++;
 349 }
 350
 351 static real get_bond_length(gmx::ArrayRef<const VsiteBondedInteraction> bonds, t_iatom ai, t_iatom aj)
 352 {
 353     real bondlen;
 354
 355     bondlen = NOTSET;
 356     for (const auto& bond : bonds)
 357     {
 358         /* check both ways */
 359         if (((ai == bond.ai()) && (aj == bond.aj())) || ((ai == bond.aj()) && (aj == bond.ai())))
 360         {
 361             bondlen = bond.parameterValue(); /* note: parameterValue might be NOTSET */
 362             break;
 363         }
 364     }
 365     return bondlen;
 366 }
 367
 368 static real get_angle(gmx::ArrayRef<const VsiteBondedInteraction> angles, t_iatom ai, t_iatom aj, t_iatom ak)
 369 {
 370     real angle;
 371
 372     angle = NOTSET;
 373     for (const auto& ang : angles)
 374     {
 375         /* check both ways */
 376         if (((ai == ang.ai()) && (aj == ang.aj()) && (ak == ang.ak()))
 377             || ((ai == ang.ak()) && (aj == ang.aj()) && (ak == ang.ai())))
 378         {
 379             angle = DEG2RAD * ang.parameterValue();
 380             break;
 381         }
 382     }
 383     return angle;
 384 }
 385
 386 static const char* get_atomtype_name_AB(t_atom* atom, PreprocessingAtomTypes* atypes)
 387 {
 388     const char* name = *atypes->atomNameFromAtomType(atom->type);
 389
 390     /* When using the decoupling option, atom types are changed
 391      * to decoupled for the non-bonded interactions, but the virtual
 392      * sites constructions should be based on the "normal" interactions.
 393      * So we return the state B atom type name if the state A atom
 394      * type is the decoupled one. We should actually check for the atom
 395      * type number, but that's not passed here. So we check for
 396      * the decoupled atom type name. This should not cause trouble
 397      * as this code is only used for topologies with v-sites without
 398      * parameters generated by pdb2gmx.
 399      */
 400     if (strcmp(name, "decoupled") == 0)
 401     {
 402         name = *atypes->atomNameFromAtomType(atom->typeB);
 403     }
 404
 405     return name;
 406 }
 407
 408 static bool calc_vsite3_param(PreprocessingAtomTypes*                     atypes,
 409                               InteractionOfType*                          vsite,
 410                               t_atoms*                                    at,
 411                               gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 412                               gmx::ArrayRef<const VsiteBondedInteraction> angles)
 413 {
 414     /* i = virtual site          |    ,k
 415      * j = 1st bonded heavy atom | i-j
 416      * k,l = 2nd bonded atoms    |    `l
 417      */
 418
 419     bool bXH3, bError;
 420     real bjk, bjl, a = -1, b = -1;
 421     /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
 422      * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 423     bXH3 = ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MNH", 3))
 424             && (gmx::equalCaseInsensitive(
 425                        get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MNH", 3)))
 426            || ((gmx::equalCaseInsensitive(
 427                        get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MCH3", 4))
 428                && (gmx::equalCaseInsensitive(
 429                           get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MCH3", 4)));
 430
 431     bjk    = get_bond_length(bonds, vsite->aj(), vsite->ak());
 432     bjl    = get_bond_length(bonds, vsite->aj(), vsite->al());
 433     bError = (bjk == NOTSET) || (bjl == NOTSET);
 434     if (bXH3)
 435     {
 436         /* now we get some XH2/XH3 group specific construction */
 437         /* note: we call the heavy atom 'C' and the X atom 'N' */
 438         real bMM, bCM, bCN, bNH, aCNH, dH, rH, dM, rM;
 439         int  aN;
 440
 441         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 442         bError = bError || (bjk != bjl);
 443
 444         /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
 445         aN = std::max(vsite->ak(), vsite->al()) + 1;
 446
 447         /* get common bonds */
 448         bMM    = get_bond_length(bonds, vsite->ak(), vsite->al());
 449         bCM    = bjk;
 450         bCN    = get_bond_length(bonds, vsite->aj(), aN);
 451         bError = bError || (bMM == NOTSET) || (bCN == NOTSET);
 452
 453         /* calculate common things */
 454         rM = 0.5 * bMM;
 455         dM = std::sqrt(gmx::square(bCM) - gmx::square(rM));
 456
 457         /* are we dealing with the X atom? */
 458         if (vsite->ai() == aN)
 459         {
 460             /* this is trivial */
 461             a = b = 0.5 * bCN / dM;
 462         }
 463         else
 464         {
 465             /* get other bondlengths and angles: */
 466             bNH    = get_bond_length(bonds, aN, vsite->ai());
 467             aCNH   = get_angle(angles, vsite->aj(), aN, vsite->ai());
 468             bError = bError || (bNH == NOTSET) || (aCNH == NOTSET);
 469
 470             /* calculate */
 471             dH = bCN - bNH * std::cos(aCNH);
 472             rH = bNH * std::sin(aCNH);
 473
 474             a = 0.5 * (dH / dM + rH / rM);
 475             b = 0.5 * (dH / dM - rH / rM);
 476         }
 477     }
 478     else
 479     {
 480         gmx_fatal(FARGS,
 481                   "calc_vsite3_param not implemented for the general case "
 482                   "(atom %d)",
 483                   vsite->ai() + 1);
 484     }
 485     vsite->setForceParameter(0, a);
 486     vsite->setForceParameter(1, b);
 487
 488     return bError;
 489 }
 490
 491 static bool calc_vsite3fd_param(InteractionOfType*                          vsite,
 492                                 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 493                                 gmx::ArrayRef<const VsiteBondedInteraction> angles)
 494 {
 495     /* i = virtual site          |    ,k
 496      * j = 1st bonded heavy atom | i-j
 497      * k,l = 2nd bonded atoms    |    `l
 498      */
 499
 500     bool bError;
 501     real bij, bjk, bjl, aijk, aijl, rk, rl;
 502
 503     bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 504     bjk  = get_bond_length(bonds, vsite->aj(), vsite->ak());
 505     bjl  = get_bond_length(bonds, vsite->aj(), vsite->al());
 506     aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 507     aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 508     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET);
 509
 510     rk = bjk * std::sin(aijk);
 511     rl = bjl * std::sin(aijl);
 512     vsite->setForceParameter(0, rk / (rk + rl));
 513     vsite->setForceParameter(1, -bij);
 514
 515     return bError;
 516 }
 517
 518 static bool calc_vsite3fad_param(InteractionOfType*                          vsite,
 519                                  gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 520                                  gmx::ArrayRef<const VsiteBondedInteraction> angles)
 521 {
 522     /* i = virtual site          |
 523      * j = 1st bonded heavy atom | i-j
 524      * k = 2nd bonded heavy atom |    `k-l
 525      * l = 3d bonded heavy atom  |
 526      */
 527
 528     bool bSwapParity, bError;
 529     real bij, aijk;
 530
 531     bSwapParity = (vsite->c1() == -1);
 532
 533     bij    = get_bond_length(bonds, vsite->ai(), vsite->aj());
 534     aijk   = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 535     bError = (bij == NOTSET) || (aijk == NOTSET);
 536
 537     vsite->setForceParameter(1, bij);
 538     vsite->setForceParameter(0, RAD2DEG * aijk);
 539
 540     if (bSwapParity)
 541     {
 542         vsite->setForceParameter(0, 360 - vsite->c0());
 543     }
 544
 545     return bError;
 546 }
 547
 548 static bool calc_vsite3out_param(PreprocessingAtomTypes*                     atypes,
 549                                  InteractionOfType*                          vsite,
 550                                  t_atoms*                                    at,
 551                                  gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 552                                  gmx::ArrayRef<const VsiteBondedInteraction> angles)
 553 {
 554     /* i = virtual site          |    ,k
 555      * j = 1st bonded heavy atom | i-j
 556      * k,l = 2nd bonded atoms    |    `l
 557      * NOTE: i is out of the j-k-l plane!
 558      */
 559
 560     bool bXH3, bError, bSwapParity;
 561     real bij, bjk, bjl, aijk, aijl, akjl, pijk, pijl, a, b, c;
 562
 563     /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
 564      * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 565     bXH3 = ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MNH", 3))
 566             && (gmx::equalCaseInsensitive(
 567                        get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MNH", 3)))
 568            || ((gmx::equalCaseInsensitive(
 569                        get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MCH3", 4))
 570                && (gmx::equalCaseInsensitive(
 571                           get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MCH3", 4)));
 572
 573     /* check if construction parity must be swapped */
 574     bSwapParity = (vsite->c1() == -1);
 575
 576     bjk    = get_bond_length(bonds, vsite->aj(), vsite->ak());
 577     bjl    = get_bond_length(bonds, vsite->aj(), vsite->al());
 578     bError = (bjk == NOTSET) || (bjl == NOTSET);
 579     if (bXH3)
 580     {
 581         /* now we get some XH3 group specific construction */
 582         /* note: we call the heavy atom 'C' and the X atom 'N' */
 583         real bMM, bCM, bCN, bNH, aCNH, dH, rH, rHx, rHy, dM, rM;
 584         int  aN;
 585
 586         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 587         bError = bError || (bjk != bjl);
 588
 589         /* the X atom (C or N) in the XH3 group is the first after the masses: */
 590         aN = std::max(vsite->ak(), vsite->al()) + 1;
 591
 592         /* get all bondlengths and angles: */
 593         bMM    = get_bond_length(bonds, vsite->ak(), vsite->al());
 594         bCM    = bjk;
 595         bCN    = get_bond_length(bonds, vsite->aj(), aN);
 596         bNH    = get_bond_length(bonds, aN, vsite->ai());
 597         aCNH   = get_angle(angles, vsite->aj(), aN, vsite->ai());
 598         bError = bError || (bMM == NOTSET) || (bCN == NOTSET) || (bNH == NOTSET) || (aCNH == NOTSET);
 599
 600         /* calculate */
 601         dH = bCN - bNH * std::cos(aCNH);
 602         rH = bNH * std::sin(aCNH);
 603         /* we assume the H's are symmetrically distributed */
 604         rHx = rH * std::cos(DEG2RAD * 30);
 605         rHy = rH * std::sin(DEG2RAD * 30);
 606         rM  = 0.5 * bMM;
 607         dM  = std::sqrt(gmx::square(bCM) - gmx::square(rM));
 608         a   = 0.5 * ((dH / dM) - (rHy / rM));
 609         b   = 0.5 * ((dH / dM) + (rHy / rM));
 610         c   = rHx / (2 * dM * rM);
 611     }
 612     else
 613     {
 614         /* this is the general construction */
 615
 616         bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 617         aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 618         aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 619         akjl = get_angle(angles, vsite->ak(), vsite->aj(), vsite->al());
 620         bError = bError || (bij == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET) || (akjl == NOTSET);
 621
 622         pijk = std::cos(aijk) * bij;
 623         pijl = std::cos(aijl) * bij;
 624         a = (pijk + (pijk * std::cos(akjl) - pijl) * std::cos(akjl) / gmx::square(std::sin(akjl))) / bjk;
 625         b = (pijl + (pijl * std::cos(akjl) - pijk) * std::cos(akjl) / gmx::square(std::sin(akjl))) / bjl;
 626         c = -std::sqrt(gmx::square(bij)
 627                        - (gmx::square(pijk) - 2 * pijk * pijl * std::cos(akjl) + gmx::square(pijl))
 628                                  / gmx::square(std::sin(akjl)))
 629             / (bjk * bjl * std::sin(akjl));
 630     }
 631
 632     vsite->setForceParameter(0, a);
 633     vsite->setForceParameter(1, b);
 634     if (bSwapParity)
 635     {
 636         vsite->setForceParameter(2, -c);
 637     }
 638     else
 639     {
 640         vsite->setForceParameter(2, c);
 641     }
 642     return bError;
 643 }
 644
 645 static bool calc_vsite4fd_param(InteractionOfType*                          vsite,
 646                                 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 647                                 gmx::ArrayRef<const VsiteBondedInteraction> angles,
 648                                 const gmx::MDLogger&                        logger)
 649 {
 650     /* i = virtual site          |    ,k
 651      * j = 1st bonded heavy atom | i-j-m
 652      * k,l,m = 2nd bonded atoms  |    `l
 653      */
 654
 655     bool bError;
 656     real bij, bjk, bjl, bjm, aijk, aijl, aijm, akjm, akjl;
 657     real pk, pl, pm, cosakl, cosakm, sinakl, sinakm, cl, cm;
 658
 659     bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 660     bjk  = get_bond_length(bonds, vsite->aj(), vsite->ak());
 661     bjl  = get_bond_length(bonds, vsite->aj(), vsite->al());
 662     bjm  = get_bond_length(bonds, vsite->aj(), vsite->am());
 663     aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 664     aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 665     aijm = get_angle(angles, vsite->ai(), vsite->aj(), vsite->am());
 666     akjm = get_angle(angles, vsite->ak(), vsite->aj(), vsite->am());
 667     akjl = get_angle(angles, vsite->ak(), vsite->aj(), vsite->al());
 668     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET) || (aijk == NOTSET)
 669              || (aijl == NOTSET) || (aijm == NOTSET) || (akjm == NOTSET) || (akjl == NOTSET);
 670
 671     if (!bError)
 672     {
 673         pk     = bjk * std::sin(aijk);
 674         pl     = bjl * std::sin(aijl);
 675         pm     = bjm * std::sin(aijm);
 676         cosakl = (std::cos(akjl) - std::cos(aijk) * std::cos(aijl)) / (std::sin(aijk) * std::sin(aijl));
 677         cosakm = (std::cos(akjm) - std::cos(aijk) * std::cos(aijm)) / (std::sin(aijk) * std::sin(aijm));
 678         if (cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1)
 679         {
 680             GMX_LOG(logger.warning)
 681                     .asParagraph()
 682                     .appendTextFormatted(
 683                             "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f",
 684                             vsite->ai() + 1,
 685                             RAD2DEG * aijk,
 686                             RAD2DEG * aijl,
 687                             RAD2DEG * aijm);
 688             gmx_fatal(FARGS,
 689                       "invalid construction in calc_vsite4fd for atom %d: "
 690                       "cosakl=%f, cosakm=%f\n",
 691                       vsite->ai() + 1,
 692                       cosakl,
 693                       cosakm);
 694         }
 695         sinakl = std::sqrt(1 - gmx::square(cosakl));
 696         sinakm = std::sqrt(1 - gmx::square(cosakm));
 697
 698         /* note: there is a '+' because of the way the sines are calculated */
 699         cl = -pk / (pl * cosakl - pk + pl * sinakl * (pm * cosakm - pk) / (pm * sinakm));
 700         cm = -pk / (pm * cosakm - pk + pm * sinakm * (pl * cosakl - pk) / (pl * sinakl));
 701
 702         vsite->setForceParameter(0, cl);
 703         vsite->setForceParameter(1, cm);
 704         vsite->setForceParameter(2, -bij);
 705     }
 706
 707     return bError;
 708 }
 709
 710
 711 static bool calc_vsite4fdn_param(InteractionOfType*                          vsite,
 712                                  gmx::ArrayRef<const VsiteBondedInteraction> bonds,
 713                                  gmx::ArrayRef<const VsiteBondedInteraction> angles,
 714                                  const gmx::MDLogger&                        logger)
 715 {
 716     /* i = virtual site          |    ,k
 717      * j = 1st bonded heavy atom | i-j-m
 718      * k,l,m = 2nd bonded atoms  |    `l
 719      */
 720
 721     bool bError;
 722     real bij, bjk, bjl, bjm, aijk, aijl, aijm;
 723     real pk, pl, pm, a, b;
 724
 725     bij  = get_bond_length(bonds, vsite->ai(), vsite->aj());
 726     bjk  = get_bond_length(bonds, vsite->aj(), vsite->ak());
 727     bjl  = get_bond_length(bonds, vsite->aj(), vsite->al());
 728     bjm  = get_bond_length(bonds, vsite->aj(), vsite->am());
 729     aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
 730     aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
 731     aijm = get_angle(angles, vsite->ai(), vsite->aj(), vsite->am());
 732
 733     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET)
 734              || (aijk == NOTSET) || (aijl == NOTSET) || (aijm == NOTSET);
 735
 736     if (!bError)
 737     {
 738
 739         /* Calculate component of bond j-k along the direction i-j */
 740         pk = -bjk * std::cos(aijk);
 741
 742         /* Calculate component of bond j-l along the direction i-j */
 743         pl = -bjl * std::cos(aijl);
 744
 745         /* Calculate component of bond j-m along the direction i-j */
 746         pm = -bjm * std::cos(aijm);
 747
 748         if (fabs(pl) < 1000 * GMX_REAL_MIN || fabs(pm) < 1000 * GMX_REAL_MIN)
 749         {
 750             GMX_LOG(logger.warning)
 751                     .asParagraph()
 752                     .appendTextFormatted(
 753                             "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f",
 754                             vsite->ai() + 1,
 755                             RAD2DEG * aijk,
 756                             RAD2DEG * aijl,
 757                             RAD2DEG * aijm);
 758             gmx_fatal(FARGS,
 759                       "invalid construction in calc_vsite4fdn for atom %d: "
 760                       "pl=%f, pm=%f\n",
 761                       vsite->ai() + 1,
 762                       pl,
 763                       pm);
 764         }
 765
 766         a = pk / pl;
 767         b = pk / pm;
 768
 769         vsite->setForceParameter(0, a);
 770         vsite->setForceParameter(1, b);
 771         vsite->setForceParameter(2, bij);
 772     }
 773
 774     return bError;
 775 }
 776
 777
 778 int set_vsites(bool                              bVerbose,
 779                t_atoms*                          atoms,
 780                PreprocessingAtomTypes*           atypes,
 781                gmx::ArrayRef<InteractionsOfType> plist,
 782                const gmx::MDLogger&              logger)
 783 {
 784     int  ftype;
 785     int  nvsite, nrset;
 786     bool bFirst, bERROR;
 787
 788     bFirst = TRUE;
 789     nvsite = 0;
 790
 791     /* Make a reverse list to avoid ninteractions^2 operations */
 792     std::vector<Atom2VsiteBond> at2vb = make_at2vsitebond(atoms->nr, plist);
 793
 794     for (ftype = 0; (ftype < F_NRE); ftype++)
 795     {
 796         if (interaction_function[ftype].flags & IF_VSITE)
 797         {
 798             nvsite += plist[ftype].size();
 799
 800             if (ftype == F_VSITEN)
 801             {
 802                 /* We don't calculate parameters for VSITEN */
 803                 continue;
 804             }
 805
 806             nrset = 0;
 807             int i = 0;
 808             for (auto& param : plist[ftype].interactionTypes)
 809             {
 810                 /* check if all parameters are set */
 811                 bool                      bSet       = true;
 812                 gmx::ArrayRef<const real> forceParam = param.forceParam();
 813                 for (int j = 0; (j < NRFP(ftype)) && bSet; j++)
 814                 {
 815                     bSet = forceParam[j] != NOTSET;
 816                 }
 817
 818                 if (debug)
 819                 {
 820                     fprintf(debug, "bSet=%s ", gmx::boolToString(bSet));
 821                     printInteractionOfType(debug, ftype, i, plist[ftype].interactionTypes[i]);
 822                 }
 823                 if (!bSet)
 824                 {
 825                     if (bVerbose && bFirst)
 826                     {
 827                         GMX_LOG(logger.info)
 828                                 .asParagraph()
 829                                 .appendTextFormatted("Calculating parameters for virtual sites");
 830                         bFirst = FALSE;
 831                     }
 832
 833                     nrset++;
 834                     /* now set the vsite parameters: */
 835                     AllVsiteBondedInteractions allVsiteBondeds =
 836                             createVsiteBondedInformation(NRAL(ftype), param.atoms(), at2vb);
 837                     if (debug)
 838                     {
 839                         fprintf(debug,
 840                                 "Found %zu bonds, %zu angles and %zu idihs "
 841                                 "for virtual site %d (%s)\n",
 842                                 allVsiteBondeds.bonds.size(),
 843                                 allVsiteBondeds.angles.size(),
 844                                 allVsiteBondeds.dihedrals.size(),
 845                                 param.ai() + 1,
 846                                 interaction_function[ftype].longname);
 847                         print_bad(debug,
 848                                   allVsiteBondeds.bonds,
 849                                   allVsiteBondeds.angles,
 850                                   allVsiteBondeds.dihedrals);
 851                     } /* debug */
 852                     switch (ftype)
 853                     {
 854                         case F_VSITE3:
 855                             bERROR = calc_vsite3_param(
 856                                     atypes, &param, atoms, allVsiteBondeds.bonds, allVsiteBondeds.angles);
 857                             break;
 858                         case F_VSITE3FD:
 859                             bERROR = calc_vsite3fd_param(
 860                                     &param, allVsiteBondeds.bonds, allVsiteBondeds.angles);
 861                             break;
 862                         case F_VSITE3FAD:
 863                             bERROR = calc_vsite3fad_param(
 864                                     &param, allVsiteBondeds.bonds, allVsiteBondeds.angles);
 865                             break;
 866                         case F_VSITE3OUT:
 867                             bERROR = calc_vsite3out_param(
 868                                     atypes, &param, atoms, allVsiteBondeds.bonds, allVsiteBondeds.angles);
 869                             break;
 870                         case F_VSITE4FD:
 871                             bERROR = calc_vsite4fd_param(
 872                                     &param, allVsiteBondeds.bonds, allVsiteBondeds.angles, logger);
 873                             break;
 874                         case F_VSITE4FDN:
 875                             bERROR = calc_vsite4fdn_param(
 876                                     &param, allVsiteBondeds.bonds, allVsiteBondeds.angles, logger);
 877                             break;
 878                         default:
 879                             gmx_fatal(FARGS,
 880                                       "Automatic parameter generation not supported "
 881                                       "for %s atom %d",
 882                                       interaction_function[ftype].longname,
 883                                       param.ai() + 1);
 884                             bERROR = TRUE;
 885                     } /* switch */
 886                     if (bERROR)
 887                     {
 888                         gmx_fatal(FARGS,
 889                                   "Automatic parameter generation not supported "
 890                                   "for %s atom %d for this bonding configuration",
 891                                   interaction_function[ftype].longname,
 892                                   param.ai() + 1);
 893                     }
 894                 } /* if bSet */
 895                 i++;
 896             }
 897         } /* if IF_VSITE */
 898     }
 899     return nvsite;
 900 }
 901
 902 void set_vsites_ptype(bool bVerbose, gmx_moltype_t* molt, const gmx::MDLogger& logger)
 903 {
 904     int ftype, i;
 905
 906     if (bVerbose)
 907     {
 908         GMX_LOG(logger.info)
 909                 .asParagraph()
 910                 .appendTextFormatted("Setting particle type to V for virtual sites");
 911     }
 912     for (ftype = 0; ftype < F_NRE; ftype++)
 913     {
 914         InteractionList* il = &molt->ilist[ftype];
 915         if (interaction_function[ftype].flags & IF_VSITE)
 916         {
 917             const int                nra = interaction_function[ftype].nratoms;
 918             const int                nrd = il->size();
 919             gmx::ArrayRef<const int> ia  = il->iatoms;
 920
 921             if (debug && nrd)
 922             {
 923                 GMX_LOG(logger.info)
 924                         .asParagraph()
 925                         .appendTextFormatted("doing %d %s virtual sites",
 926                                              (nrd / (nra + 1)),
 927                                              interaction_function[ftype].longname);
 928             }
 929
 930             for (i = 0; (i < nrd);)
 931             {
 932                 /* The virtual site */
 933                 int avsite                     = ia[i + 1];
 934                 molt->atoms.atom[avsite].ptype = eptVSite;
 935
 936                 i += nra + 1;
 937             }
 938         }
 939     }
 940 }
 941
 942 /*! \brief
 943  *  Convenience typedef for linking function type to parameter numbers.
 944  *
 945  *  The entries in this datastructure are valid if the particle participates in
 946  *  a virtual site interaction and has a valid vsite function type other than VSITEN.
 947  *  \todo Change to remove empty constructor when gmx::compat::optional is available.
 948  */
 949 class VsiteAtomMapping
 950 {
 951 public:
 952     //! Only construct with all information in place or nothing
 953     VsiteAtomMapping(int functionType, int interactionIndex) :
 954         functionType_(functionType),
 955         interactionIndex_(interactionIndex)
 956     {
 957     }
 958     VsiteAtomMapping() : functionType_(-1), interactionIndex_(-1) {}
 959     //! Get function type.
 960     const int& functionType() const { return functionType_; }
 961     //! Get parameter number.
 962     const int& interactionIndex() const { return interactionIndex_; }
 963
 964 private:
 965     //! Function type for the linked parameter.
 966     int functionType_;
 967     //! The linked parameter.
 968     int interactionIndex_;
 969 };
 970
 971 static void check_vsite_constraints(gmx::ArrayRef<InteractionsOfType> plist,
 972                                     int                               cftype,
 973                                     const int                         vsite_type[],
 974                                     const gmx::MDLogger&              logger)
 975 {
 976     int n = 0;
 977     for (const auto& param : plist[cftype].interactionTypes)
 978     {
 979         gmx::ArrayRef<const int> atoms = param.atoms();
 980         for (int k = 0; k < 2; k++)
 981         {
 982             int atom = atoms[k];
 983             if (vsite_type[atom] != NOTSET)
 984             {
 985                 GMX_LOG(logger.info)
 986                         .asParagraph()
 987                         .appendTextFormatted(
 988                                 "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)",
 989                                 param.ai() + 1,
 990                                 param.aj() + 1,
 991                                 atom + 1);
 992                 n++;
 993             }
 994         }
 995     }
 996     if (n)
 997     {
 998         gmx_fatal(FARGS, "There were %d virtual sites involved in constraints", n);
 999     }
1000 }
1001
1002 static void clean_vsite_bonds(gmx::ArrayRef<InteractionsOfType>     plist,
1003                               gmx::ArrayRef<const VsiteAtomMapping> pindex,
1004                               int                                   cftype,
1005                               const int                             vsite_type[],
1006                               const gmx::MDLogger&                  logger)
1007 {
1008     int                 ftype, nOut;
1009     int                 nconverted, nremoved;
1010     int                 oatom, at1, at2;
1011     bool                bKeep, bRemove, bAllFD;
1012     InteractionsOfType* ps;
1013
1014     if (cftype == F_CONNBONDS)
1015     {
1016         return;
1017     }
1018
1019     ps         = &(plist[cftype]);
1020     nconverted = 0;
1021     nremoved   = 0;
1022     nOut       = 0;
1023     for (auto bond = ps->interactionTypes.begin(); bond != ps->interactionTypes.end();)
1024     {
1025         int        vsnral      = 0;
1026         const int* first_atoms = nullptr;
1027
1028         bKeep   = false;
1029         bRemove = false;
1030         bAllFD  = true;
1031         /* check if all virtual sites are constructed from the same atoms */
1032         int                      nvsite = 0;
1033         gmx::ArrayRef<const int> atoms  = bond->atoms();
1034         for (int k = 0; (k < 2) && !bKeep && !bRemove; k++)
1035         {
1036             /* for all atoms in the bond */
1037             int atom = atoms[k];
1038             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1039             {
1040                 nvsite++;
1041                 bool bThisFD  = ((pindex[atom].functionType() == F_VSITE3FD)
1042                                 || (pindex[atom].functionType() == F_VSITE3FAD)
1043                                 || (pindex[atom].functionType() == F_VSITE4FD)
1044                                 || (pindex[atom].functionType() == F_VSITE4FDN));
1045                 bool bThisOUT = ((pindex[atom].functionType() == F_VSITE3OUT)
1046                                  && ((interaction_function[cftype].flags & IF_CONSTRAINT) != 0U));
1047                 bAllFD        = bAllFD && bThisFD;
1048                 if (bThisFD || bThisOUT)
1049                 {
1050                     oatom = atoms[1 - k]; /* the other atom */
1051                     if (vsite_type[oatom] == NOTSET
1052                         && oatom
1053                                    == plist[pindex[atom].functionType()]
1054                                               .interactionTypes[pindex[atom].interactionIndex()]
1055                                               .aj())
1056                     {
1057                         /* if the other atom isn't a vsite, and it is AI */
1058                         bRemove = true;
1059                         if (bThisOUT)
1060                         {
1061                             nOut++;
1062                         }
1063                     }
1064                 }
1065                 if (!bRemove)
1066                 {
1067                     /* TODO This fragment, and corresponding logic in
1068                        clean_vsite_angles and clean_vsite_dihs should
1069                        be refactored into a common function */
1070                     if (nvsite == 1)
1071                     {
1072                         /* if this is the first vsite we encounter then
1073                            store construction atoms */
1074                         /* TODO This would be nicer to implement with
1075                            a C++ "vector view" class" with an
1076                            STL-container-like interface. */
1077                         vsnral      = NRAL(pindex[atom].functionType()) - 1;
1078                         first_atoms = plist[pindex[atom].functionType()]
1079                                               .interactionTypes[pindex[atom].interactionIndex()]
1080                                               .atoms()
1081                                               .data()
1082                                       + 1;
1083                     }
1084                     else
1085                     {
1086                         GMX_ASSERT(vsnral != 0, "nvsite > 1 must have vsnral != 0");
1087                         GMX_ASSERT(first_atoms != nullptr,
1088                                    "nvsite > 1 must have first_atoms != NULL");
1089                         /* if it is not the first then
1090                            check if this vsite is constructed from the same atoms */
1091                         if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1092                         {
1093                             for (int m = 0; (m < vsnral) && !bKeep; m++)
1094                             {
1095                                 const int* atoms;
1096
1097                                 bool bPresent = false;
1098                                 atoms         = plist[pindex[atom].functionType()]
1099                                                 .interactionTypes[pindex[atom].interactionIndex()]
1100                                                 .atoms()
1101                                                 .data()
1102                                         + 1;
1103                                 for (int n = 0; (n < vsnral) && !bPresent; n++)
1104                                 {
1105                                     if (atoms[m] == first_atoms[n])
1106                                     {
1107                                         bPresent = true;
1108                                     }
1109                                 }
1110                                 if (!bPresent)
1111                                 {
1112                                     bKeep = true;
1113                                 }
1114                             }
1115                         }
1116                         else
1117                         {
1118                             bKeep = true;
1119                         }
1120                     }
1121                 }
1122             }
1123         }
1124
1125         if (bRemove)
1126         {
1127             bKeep = false;
1128         }
1129         else
1130         {
1131             /* if we have no virtual sites in this bond, keep it */
1132             if (nvsite == 0)
1133             {
1134                 bKeep = true;
1135             }
1136
1137             /* TODO This loop and the corresponding loop in
1138                check_vsite_angles should be refactored into a common
1139                function */
1140             /* check if all non-vsite atoms are used in construction: */
1141             bool bFirstTwo = true;
1142             for (int k = 0; (k < 2) && !bKeep; k++) /* for all atoms in the bond */
1143             {
1144                 int atom = atoms[k];
1145                 if (vsite_type[atom] == NOTSET)
1146                 {
1147                     bool bUsed = false;
1148                     for (int m = 0; (m < vsnral) && !bUsed; m++)
1149                     {
1150                         GMX_ASSERT(first_atoms != nullptr,
1151                                    "If we've seen a vsite before, we know what its first atom "
1152                                    "index was");
1153
1154                         if (atom == first_atoms[m])
1155                         {
1156                             bUsed     = true;
1157                             bFirstTwo = bFirstTwo && m < 2;
1158                         }
1159                     }
1160                     if (!bUsed)
1161                     {
1162                         bKeep = true;
1163                     }
1164                 }
1165             }
1166
1167             if (!(bAllFD && bFirstTwo))
1168             {
1169                 /* Two atom bonded interactions include constraints.
1170                  * We need to remove constraints between vsite pairs that have
1171                  * a fixed distance due to being constructed from the same
1172                  * atoms, since this can be numerically unstable.
1173                  */
1174                 for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1175                 {
1176                     at1           = first_atoms[m];
1177                     at2           = first_atoms[(m + 1) % vsnral];
1178                     bool bPresent = false;
1179                     for (ftype = 0; ftype < F_NRE; ftype++)
1180                     {
1181                         if (interaction_function[ftype].flags & IF_CONSTRAINT)
1182                         {
1183                             for (auto entry = plist[ftype].interactionTypes.begin();
1184                                  (entry != plist[ftype].interactionTypes.end()) && !bPresent;
1185                                  entry++)
1186                             {
1187                                 /* all constraints until one matches */
1188                                 bPresent = (((entry->ai() == at1) && (entry->aj() == at2))
1189                                             || ((entry->ai() == at2) && (entry->aj() == at1)));
1190                             }
1191                         }
1192                     }
1193                     if (!bPresent)
1194                     {
1195                         bKeep = true;
1196                     }
1197                 }
1198             }
1199         }
1200
1201         if (bKeep)
1202         {
1203             ++bond;
1204         }
1205         else if (IS_CHEMBOND(cftype))
1206         {
1207             plist[F_CONNBONDS].interactionTypes.emplace_back(*bond);
1208             bond = ps->interactionTypes.erase(bond);
1209             nconverted++;
1210         }
1211         else
1212         {
1213             bond = ps->interactionTypes.erase(bond);
1214             nremoved++;
1215         }
1216     }
1217
1218     if (nremoved)
1219     {
1220         GMX_LOG(logger.info)
1221                 .asParagraph()
1222                 .appendTextFormatted("Removed   %4d %15ss with virtual sites, %zu left",
1223                                      nremoved,
1224                                      interaction_function[cftype].longname,
1225                                      ps->size());
1226     }
1227     if (nconverted)
1228     {
1229         GMX_LOG(logger.info)
1230                 .asParagraph()
1231                 .appendTextFormatted(
1232                         "Converted %4d %15ss with virtual sites to connections, %zu left",
1233                         nconverted,
1234                         interaction_function[cftype].longname,
1235                         ps->size());
1236     }
1237     if (nOut)
1238     {
1239         GMX_LOG(logger.info)
1240                 .asParagraph()
1241                 .appendTextFormatted(
1242                         "Warning: removed %d %ss with vsite with %s construction\n"
1243                         "         This vsite construction does not guarantee constant "
1244                         "bond-length\n"
1245                         "         If the constructions were generated by pdb2gmx ignore "
1246                         "this warning",
1247                         nOut,
1248                         interaction_function[cftype].longname,
1249                         interaction_function[F_VSITE3OUT].longname);
1250     }
1251 }
1252
1253 static void clean_vsite_angles(gmx::ArrayRef<InteractionsOfType>         plist,
1254                                gmx::ArrayRef<VsiteAtomMapping>           pindex,
1255                                int                                       cftype,
1256                                const int                                 vsite_type[],
1257                                gmx::ArrayRef<const Atom2VsiteConnection> at2vc,
1258                                const gmx::MDLogger&                      logger)
1259 {
1260     int                 atom, at1, at2;
1261     InteractionsOfType* ps;
1262
1263     ps          = &(plist[cftype]);
1264     int oldSize = ps->size();
1265     for (auto angle = ps->interactionTypes.begin(); angle != ps->interactionTypes.end();)
1266     {
1267         int        vsnral      = 0;
1268         const int* first_atoms = nullptr;
1269
1270         bool bKeep    = false;
1271         bool bAll3FAD = true;
1272         /* check if all virtual sites are constructed from the same atoms */
1273         int                      nvsite = 0;
1274         gmx::ArrayRef<const int> atoms  = angle->atoms();
1275         for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1276         {
1277             int atom = atoms[k];
1278             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1279             {
1280                 nvsite++;
1281                 bAll3FAD = bAll3FAD && (pindex[atom].functionType() == F_VSITE3FAD);
1282                 if (nvsite == 1)
1283                 {
1284                     /* store construction atoms of first vsite */
1285                     vsnral      = NRAL(pindex[atom].functionType()) - 1;
1286                     first_atoms = plist[pindex[atom].functionType()]
1287                                           .interactionTypes[pindex[atom].interactionIndex()]
1288                                           .atoms()
1289                                           .data()
1290                                   + 1;
1291                 }
1292                 else
1293                 {
1294                     GMX_ASSERT(vsnral != 0,
1295                                "If we've seen a vsite before, we know how many constructing atoms "
1296                                "it had");
1297                     GMX_ASSERT(
1298                             first_atoms != nullptr,
1299                             "If we've seen a vsite before, we know what its first atom index was");
1300                     /* check if this vsite is constructed from the same atoms */
1301                     if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1302                     {
1303                         for (int m = 0; (m < vsnral) && !bKeep; m++)
1304                         {
1305                             const int* subAtoms;
1306
1307                             bool bPresent = false;
1308                             subAtoms      = plist[pindex[atom].functionType()]
1309                                                .interactionTypes[pindex[atom].interactionIndex()]
1310                                                .atoms()
1311                                                .data()
1312                                        + 1;
1313                             for (int n = 0; (n < vsnral) && !bPresent; n++)
1314                             {
1315                                 if (subAtoms[m] == first_atoms[n])
1316                                 {
1317                                     bPresent = true;
1318                                 }
1319                             }
1320                             if (!bPresent)
1321                             {
1322                                 bKeep = true;
1323                             }
1324                         }
1325                     }
1326                     else
1327                     {
1328                         bKeep = true;
1329                     }
1330                 }
1331             }
1332         }
1333
1334         /* keep all angles with no virtual sites in them or
1335            with virtual sites with more than 3 constr. atoms */
1336         if (nvsite == 0 && vsnral > 3)
1337         {
1338             bKeep = true;
1339         }
1340
1341         /* check if all non-vsite atoms are used in construction: */
1342         bool bFirstTwo = true;
1343         for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1344         {
1345             atom = atoms[k];
1346             if (vsite_type[atom] == NOTSET)
1347             {
1348                 bool bUsed = false;
1349                 for (int m = 0; (m < vsnral) && !bUsed; m++)
1350                 {
1351                     GMX_ASSERT(
1352                             first_atoms != nullptr,
1353                             "If we've seen a vsite before, we know what its first atom index was");
1354
1355                     if (atom == first_atoms[m])
1356                     {
1357                         bUsed     = true;
1358                         bFirstTwo = bFirstTwo && m < 2;
1359                     }
1360                 }
1361                 if (!bUsed)
1362                 {
1363                     bKeep = true;
1364                 }
1365             }
1366         }
1367
1368         if (!(bAll3FAD && bFirstTwo))
1369         {
1370             /* check if all constructing atoms are constrained together */
1371             for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1372             {
1373                 at1           = first_atoms[m];
1374                 at2           = first_atoms[(m + 1) % vsnral];
1375                 bool bPresent = false;
1376                 auto found    = std::find(at2vc[at1].begin(), at2vc[at1].end(), at2);
1377                 if (found != at2vc[at1].end())
1378                 {
1379                     bPresent = true;
1380                 }
1381                 if (!bPresent)
1382                 {
1383                     bKeep = true;
1384                 }
1385             }
1386         }
1387
1388         if (bKeep)
1389         {
1390             ++angle;
1391         }
1392         else
1393         {
1394             angle = ps->interactionTypes.erase(angle);
1395         }
1396     }
1397
1398     if (oldSize != gmx::ssize(*ps))
1399     {
1400         GMX_LOG(logger.info)
1401                 .asParagraph()
1402                 .appendTextFormatted("Removed   %4zu %15ss with virtual sites, %zu left",
1403                                      oldSize - ps->size(),
1404                                      interaction_function[cftype].longname,
1405                                      ps->size());
1406     }
1407 }
1408
1409 static void clean_vsite_dihs(gmx::ArrayRef<InteractionsOfType>     plist,
1410                              gmx::ArrayRef<const VsiteAtomMapping> pindex,
1411                              int                                   cftype,
1412                              const int                             vsite_type[],
1413                              const gmx::MDLogger&                  logger)
1414 {
1415     InteractionsOfType* ps;
1416
1417     ps = &(plist[cftype]);
1418
1419     int oldSize = ps->size();
1420     for (auto dih = ps->interactionTypes.begin(); dih != ps->interactionTypes.end();)
1421     {
1422         int        vsnral      = 0;
1423         const int* first_atoms = nullptr;
1424         int        atom;
1425
1426         gmx::ArrayRef<const int> atoms = dih->atoms();
1427         bool                     bKeep = false;
1428         /* check if all virtual sites are constructed from the same atoms */
1429         int nvsite = 0;
1430         for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1431         {
1432             atom = atoms[k];
1433             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1434             {
1435                 if (nvsite == 0)
1436                 {
1437                     /* store construction atoms of first vsite */
1438                     vsnral      = NRAL(pindex[atom].functionType()) - 1;
1439                     first_atoms = plist[pindex[atom].functionType()]
1440                                           .interactionTypes[pindex[atom].interactionIndex()]
1441                                           .atoms()
1442                                           .data()
1443                                   + 1;
1444                 }
1445                 else
1446                 {
1447                     GMX_ASSERT(vsnral != 0,
1448                                "If we've seen a vsite before, we know how many constructing atoms "
1449                                "it had");
1450                     GMX_ASSERT(
1451                             first_atoms != nullptr,
1452                             "If we've seen a vsite before, we know what its first atom index was");
1453                     /* check if this vsite is constructed from the same atoms */
1454                     if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1455                     {
1456                         for (int m = 0; (m < vsnral) && !bKeep; m++)
1457                         {
1458                             const int* subAtoms;
1459
1460                             bool bPresent = false;
1461                             subAtoms      = plist[pindex[atom].functionType()]
1462                                                .interactionTypes[pindex[atom].interactionIndex()]
1463                                                .atoms()
1464                                                .data()
1465                                        + 1;
1466                             for (int n = 0; (n < vsnral) && !bPresent; n++)
1467                             {
1468                                 if (subAtoms[m] == first_atoms[n])
1469                                 {
1470                                     bPresent = true;
1471                                 }
1472                             }
1473                             if (!bPresent)
1474                             {
1475                                 bKeep = true;
1476                             }
1477                         }
1478                     }
1479                 }
1480                 /* TODO clean_site_bonds and _angles do this increment
1481                    at the top of the loop. Refactor this for
1482                    consistency */
1483                 nvsite++;
1484             }
1485         }
1486
1487         /* keep all dihedrals with no virtual sites in them */
1488         if (nvsite == 0)
1489         {
1490             bKeep = true;
1491         }
1492
1493         /* check if all atoms in dihedral are either virtual sites, or used in
1494            construction of virtual sites. If so, keep it, if not throw away: */
1495         for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1496         {
1497             GMX_ASSERT(vsnral != 0,
1498                        "If we've seen a vsite before, we know how many constructing atoms it had");
1499             GMX_ASSERT(first_atoms != nullptr,
1500                        "If we've seen a vsite before, we know what its first atom index was");
1501             atom = atoms[k];
1502             if (vsite_type[atom] == NOTSET)
1503             {
1504                 /* vsnral will be set here, we don't get here with nvsite==0 */
1505                 bool bUsed = false;
1506                 for (int m = 0; (m < vsnral) && !bUsed; m++)
1507                 {
1508                     if (atom == first_atoms[m])
1509                     {
1510                         bUsed = true;
1511                     }
1512                 }
1513                 if (!bUsed)
1514                 {
1515                     bKeep = true;
1516                 }
1517             }
1518         }
1519
1520         if (bKeep)
1521         {
1522             ++dih;
1523         }
1524         else
1525         {
1526             dih = ps->interactionTypes.erase(dih);
1527         }
1528     }
1529
1530     if (oldSize != gmx::ssize(*ps))
1531     {
1532         GMX_LOG(logger.info)
1533                 .asParagraph()
1534                 .appendTextFormatted("Removed   %4zu %15ss with virtual sites, %zu left",
1535                                      oldSize - ps->size(),
1536                                      interaction_function[cftype].longname,
1537                                      ps->size());
1538     }
1539 }
1540
1541 // TODO use gmx::compat::optional for pindex.
1542 void clean_vsite_bondeds(gmx::ArrayRef<InteractionsOfType> plist,
1543                          int                               natoms,
1544                          bool                              bRmVSiteBds,
1545                          const gmx::MDLogger&              logger)
1546 {
1547     int                               nvsite, vsite;
1548     int*                              vsite_type;
1549     std::vector<VsiteAtomMapping>     pindex;
1550     std::vector<Atom2VsiteConnection> at2vc;
1551
1552     /* make vsite_type array */
1553     snew(vsite_type, natoms);
1554     for (int i = 0; i < natoms; i++)
1555     {
1556         vsite_type[i] = NOTSET;
1557     }
1558     nvsite = 0;
1559     for (int ftype = 0; ftype < F_NRE; ftype++)
1560     {
1561         if (interaction_function[ftype].flags & IF_VSITE)
1562         {
1563             nvsite += plist[ftype].size();
1564             int i = 0;
1565             while (i < gmx::ssize(plist[ftype]))
1566             {
1567                 vsite = plist[ftype].interactionTypes[i].ai();
1568                 if (vsite_type[vsite] == NOTSET)
1569                 {
1570                     vsite_type[vsite] = ftype;
1571                 }
1572                 else
1573                 {
1574                     gmx_fatal(FARGS, "multiple vsite constructions for atom %d", vsite + 1);
1575                 }
1576                 if (ftype == F_VSITEN)
1577                 {
1578                     while (i < gmx::ssize(plist[ftype]) && plist[ftype].interactionTypes[i].ai() == vsite)
1579                     {
1580                         i++;
1581                     }
1582                 }
1583                 else
1584                 {
1585                     i++;
1586                 }
1587             }
1588         }
1589     }
1590
1591     /* the rest only if we have virtual sites: */
1592     if (nvsite)
1593     {
1594         GMX_LOG(logger.info)
1595                 .asParagraph()
1596                 .appendTextFormatted("Cleaning up constraints %swith virtual sites",
1597                                      bRmVSiteBds ? "and constant bonded interactions " : "");
1598
1599         /* Make a reverse list to avoid ninteractions^2 operations */
1600         at2vc = make_at2vsitecon(natoms, plist);
1601
1602         pindex.resize(natoms);
1603         for (int ftype = 0; ftype < F_NRE; ftype++)
1604         {
1605             /* Here we skip VSITEN. In neary all practical use cases this
1606              * is not an issue, since VSITEN is intended for constructing
1607              * additional interaction sites, not for replacing normal atoms
1608              * with bonded interactions. Thus we do not expect constant
1609              * bonded interactions. If VSITEN does get used with constant
1610              * bonded interactions, these are not removed which only leads
1611              * to very minor extra computation and constant energy.
1612              * The only problematic case is onstraints between VSITEN
1613              * constructions with fixed distance (which is anyhow useless).
1614              * This will generate a fatal error in check_vsite_constraints.
1615              */
1616             if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
1617             {
1618                 for (gmx::index interactionIndex = 0; interactionIndex < gmx::ssize(plist[ftype]);
1619                      interactionIndex++)
1620                 {
1621                     int k     = plist[ftype].interactionTypes[interactionIndex].ai();
1622                     pindex[k] = VsiteAtomMapping(ftype, interactionIndex);
1623                 }
1624             }
1625         }
1626
1627         /* remove interactions that include virtual sites */
1628         for (int ftype = 0; ftype < F_NRE; ftype++)
1629         {
1630             if (((interaction_function[ftype].flags & IF_BOND) && bRmVSiteBds)
1631                 || (interaction_function[ftype].flags & IF_CONSTRAINT))
1632             {
1633                 if (interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT))
1634                 {
1635                     clean_vsite_bonds(plist, pindex, ftype, vsite_type, logger);
1636                 }
1637                 else if (interaction_function[ftype].flags & IF_ATYPE)
1638                 {
1639                     clean_vsite_angles(plist, pindex, ftype, vsite_type, at2vc, logger);
1640                 }
1641                 else if ((ftype == F_PDIHS) || (ftype == F_IDIHS))
1642                 {
1643                     clean_vsite_dihs(plist, pindex, ftype, vsite_type, logger);
1644                 }
1645             }
1646         }
1647         /* check that no remaining constraints include virtual sites */
1648         for (int ftype = 0; ftype < F_NRE; ftype++)
1649         {
1650             if (interaction_function[ftype].flags & IF_CONSTRAINT)
1651             {
1652                 check_vsite_constraints(plist, ftype, vsite_type, logger);
1653             }
1654         }
1655     }
1656     sfree(vsite_type);
1657 }