src/gromacs/gmxpreprocess/vsite_parm.cpp

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