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