src/gromacs/gmxpreprocess/vsite_parm.cpp

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