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