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