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