src/gromacs/gmxpreprocess/vsite_parm.cpp

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