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