src/gromacs/gmxpreprocess/vsite_parm.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "vsite_parm.h"
  40
  41 #include <cmath>
  42 #include <cstdio>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/gmxpreprocess/add_par.h"
  48 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  49 #include "gromacs/gmxpreprocess/grompp_impl.h"
  50 #include "gromacs/gmxpreprocess/notset.h"
  51 #include "gromacs/gmxpreprocess/toputil.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/mdtypes/md_enums.h"
  56 #include "gromacs/topology/ifunc.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/utility/arrayref.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/gmxassert.h"
  62 #include "gromacs/utility/smalloc.h"
  63 #include "gromacs/utility/strconvert.h"
  64
  65 #include "resall.h"
  66
  67 typedef struct {
  68     t_iatom  a[4];
  69     real     c;
  70     t_iatom &ai() { return a[0]; }
  71     t_iatom &aj() { return a[1]; }
  72     t_iatom &ak() { return a[2]; }
  73     t_iatom &al() { return a[3]; }
  74 } t_mybonded;
  75
  76 typedef struct {
  77     int      ftype;
  78     t_param *param;
  79 } vsitebondparam_t;
  80
  81 typedef struct {
  82     int               nr;
  83     int               ftype;
  84     vsitebondparam_t *vsbp;
  85 } at2vsitebond_t;
  86
  87 typedef struct {
  88     int  nr;
  89     int *aj;
  90 } at2vsitecon_t;
  91
  92 static int vsite_bond_nrcheck(int ftype)
  93 {
  94     int nrcheck;
  95
  96     if ((interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT | IF_ATYPE)) || (ftype == F_IDIHS))
  97     {
  98         nrcheck = NRAL(ftype);
  99     }
 100     else
 101     {
 102         nrcheck = 0;
 103     }
 104
 105     return nrcheck;
 106 }
 107
 108 static void enter_bonded(int nratoms, int *nrbonded, t_mybonded **bondeds,
 109                          t_param *param)
 110 {
 111     int j;
 112
 113     srenew(*bondeds, *nrbonded+1);
 114
 115     /* copy atom numbers */
 116     for (j = 0; j < nratoms; j++)
 117     {
 118         (*bondeds)[*nrbonded].a[j] = param->a[j];
 119     }
 120     /* copy parameter */
 121     (*bondeds)[*nrbonded].c = param->c0();
 122
 123     (*nrbonded)++;
 124 }
 125
 126 static void get_bondeds(int nrat, const t_iatom atoms[],
 127                         at2vsitebond_t *at2vb,
 128                         int *nrbond, t_mybonded **bonds,
 129                         int *nrang,  t_mybonded **angles,
 130                         int *nridih, t_mybonded **idihs )
 131 {
 132     int      k, i, ftype, nrcheck;
 133     t_param *param;
 134
 135     for (k = 0; k < nrat; k++)
 136     {
 137         for (i = 0; i < at2vb[atoms[k]].nr; i++)
 138         {
 139             ftype   = at2vb[atoms[k]].vsbp[i].ftype;
 140             param   = at2vb[atoms[k]].vsbp[i].param;
 141             nrcheck = vsite_bond_nrcheck(ftype);
 142             /* abuse nrcheck to see if we're adding bond, angle or idih */
 143             switch (nrcheck)
 144             {
 145                 case 2: enter_bonded(nrcheck, nrbond, bonds, param); break;
 146                 case 3: enter_bonded(nrcheck, nrang, angles, param); break;
 147                 case 4: enter_bonded(nrcheck, nridih, idihs, param); break;
 148             }
 149         }
 150     }
 151 }
 152
 153 static at2vsitebond_t *make_at2vsitebond(int natoms, gmx::ArrayRef<InteractionTypeParameters> plist)
 154 {
 155     bool           *bVSI;
 156     int             ftype, i, j, nrcheck, nr;
 157     t_iatom        *aa;
 158     at2vsitebond_t *at2vb;
 159
 160     snew(at2vb, natoms);
 161
 162     snew(bVSI, natoms);
 163     for (ftype = 0; (ftype < F_NRE); ftype++)
 164     {
 165         if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
 166         {
 167             for (i = 0; (i < plist[ftype].nr); i++)
 168             {
 169                 for (j = 0; j < NRAL(ftype); j++)
 170                 {
 171                     bVSI[plist[ftype].param[i].a[j]] = TRUE;
 172                 }
 173             }
 174         }
 175     }
 176
 177     for (ftype = 0; (ftype < F_NRE); ftype++)
 178     {
 179         nrcheck = vsite_bond_nrcheck(ftype);
 180         if (nrcheck > 0)
 181         {
 182             for (i = 0; (i < plist[ftype].nr); i++)
 183             {
 184                 aa = plist[ftype].param[i].a;
 185                 for (j = 0; j < nrcheck; j++)
 186                 {
 187                     if (bVSI[aa[j]])
 188                     {
 189                         nr = at2vb[aa[j]].nr;
 190                         if (nr % 10 == 0)
 191                         {
 192                             srenew(at2vb[aa[j]].vsbp, nr+10);
 193                         }
 194                         at2vb[aa[j]].vsbp[nr].ftype = ftype;
 195                         at2vb[aa[j]].vsbp[nr].param = &plist[ftype].param[i];
 196                         at2vb[aa[j]].nr++;
 197                     }
 198                 }
 199             }
 200         }
 201     }
 202     sfree(bVSI);
 203
 204     return at2vb;
 205 }
 206
 207 static void done_at2vsitebond(int natoms, at2vsitebond_t *at2vb)
 208 {
 209     int i;
 210
 211     for (i = 0; i < natoms; i++)
 212     {
 213         if (at2vb[i].nr)
 214         {
 215             sfree(at2vb[i].vsbp);
 216         }
 217     }
 218     sfree(at2vb);
 219 }
 220
 221 static at2vsitecon_t *make_at2vsitecon(int natoms, gmx::ArrayRef<InteractionTypeParameters> plist)
 222 {
 223     bool          *bVSI;
 224     int            ftype, i, j, ai, aj, nr;
 225     at2vsitecon_t *at2vc;
 226
 227     snew(at2vc, natoms);
 228
 229     snew(bVSI, natoms);
 230     for (ftype = 0; (ftype < F_NRE); ftype++)
 231     {
 232         if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
 233         {
 234             for (i = 0; (i < plist[ftype].nr); i++)
 235             {
 236                 for (j = 0; j < NRAL(ftype); j++)
 237                 {
 238                     bVSI[plist[ftype].param[i].a[j]] = TRUE;
 239                 }
 240             }
 241         }
 242     }
 243
 244     for (ftype = 0; (ftype < F_NRE); ftype++)
 245     {
 246         if (interaction_function[ftype].flags & IF_CONSTRAINT)
 247         {
 248             for (i = 0; (i < plist[ftype].nr); i++)
 249             {
 250                 ai = plist[ftype].param[i].ai();
 251                 aj = plist[ftype].param[i].aj();
 252                 if (bVSI[ai] && bVSI[aj])
 253                 {
 254                     /* Store forward direction */
 255                     nr = at2vc[ai].nr;
 256                     if (nr % 10 == 0)
 257                     {
 258                         srenew(at2vc[ai].aj, nr+10);
 259                     }
 260                     at2vc[ai].aj[nr] = aj;
 261                     at2vc[ai].nr++;
 262                     /* Store backward direction */
 263                     nr = at2vc[aj].nr;
 264                     if (nr % 10 == 0)
 265                     {
 266                         srenew(at2vc[aj].aj, nr+10);
 267                     }
 268                     at2vc[aj].aj[nr] = ai;
 269                     at2vc[aj].nr++;
 270                 }
 271             }
 272         }
 273     }
 274     sfree(bVSI);
 275
 276     return at2vc;
 277 }
 278
 279 static void done_at2vsitecon(int natoms, at2vsitecon_t *at2vc)
 280 {
 281     int i;
 282
 283     for (i = 0; i < natoms; i++)
 284     {
 285         if (at2vc[i].nr)
 286         {
 287             sfree(at2vc[i].aj);
 288         }
 289     }
 290     sfree(at2vc);
 291 }
 292
 293 /* for debug */
 294 static void print_bad(FILE *fp,
 295                       int nrbond, t_mybonded *bonds,
 296                       int nrang,  t_mybonded *angles,
 297                       int nridih, t_mybonded *idihs )
 298 {
 299     int i;
 300
 301     if (nrbond)
 302     {
 303         fprintf(fp, "bonds:");
 304         for (i = 0; i < nrbond; i++)
 305         {
 306             fprintf(fp, " %d-%d (%g)",
 307                     bonds[i].ai()+1, bonds[i].aj()+1, bonds[i].c);
 308         }
 309         fprintf(fp, "\n");
 310     }
 311     if (nrang)
 312     {
 313         fprintf(fp, "angles:");
 314         for (i = 0; i < nrang; i++)
 315         {
 316             fprintf(fp, " %d-%d-%d (%g)",
 317                     angles[i].ai()+1, angles[i].aj()+1,
 318                     angles[i].ak()+1, angles[i].c);
 319         }
 320         fprintf(fp, "\n");
 321     }
 322     if (nridih)
 323     {
 324         fprintf(fp, "idihs:");
 325         for (i = 0; i < nridih; i++)
 326         {
 327             fprintf(fp, " %d-%d-%d-%d (%g)",
 328                     idihs[i].ai()+1, idihs[i].aj()+1,
 329                     idihs[i].ak()+1, idihs[i].al()+1, idihs[i].c);
 330         }
 331         fprintf(fp, "\n");
 332     }
 333 }
 334
 335 static void print_param(FILE *fp, int ftype, int i, t_param *param)
 336 {
 337     static int pass       = 0;
 338     static int prev_ftype = NOTSET;
 339     static int prev_i     = NOTSET;
 340     int        j;
 341
 342     if ( (ftype != prev_ftype) || (i != prev_i) )
 343     {
 344         pass       = 0;
 345         prev_ftype = ftype;
 346         prev_i     = i;
 347     }
 348     fprintf(fp, "(%d) plist[%s].param[%d]",
 349             pass, interaction_function[ftype].name, i);
 350     for (j = 0; j < NRFP(ftype); j++)
 351     {
 352         fprintf(fp, ".c[%d]=%g ", j, param->c[j]);
 353     }
 354     fprintf(fp, "\n");
 355     pass++;
 356 }
 357
 358 static real get_bond_length(int nrbond, t_mybonded bonds[],
 359                             t_iatom ai, t_iatom aj)
 360 {
 361     int  i;
 362     real bondlen;
 363
 364     bondlen = NOTSET;
 365     for (i = 0; i < nrbond && (bondlen == NOTSET); i++)
 366     {
 367         /* check both ways */
 368         if ( ( (ai == bonds[i].ai()) && (aj == bonds[i].aj()) ) ||
 369              ( (ai == bonds[i].aj()) && (aj == bonds[i].ai()) ) )
 370         {
 371             bondlen = bonds[i].c; /* note: bonds[i].c might be NOTSET */
 372         }
 373     }
 374     return bondlen;
 375 }
 376
 377 static real get_angle(int nrang, t_mybonded angles[],
 378                       t_iatom ai, t_iatom aj, t_iatom ak)
 379 {
 380     int  i;
 381     real angle;
 382
 383     angle = NOTSET;
 384     for (i = 0; i < nrang && (angle == NOTSET); i++)
 385     {
 386         /* check both ways */
 387         if ( ( (ai == angles[i].ai()) && (aj == angles[i].aj()) && (ak == angles[i].ak()) ) ||
 388              ( (ai == angles[i].ak()) && (aj == angles[i].aj()) && (ak == angles[i].ai()) ) )
 389         {
 390             angle = DEG2RAD*angles[i].c;
 391         }
 392     }
 393     return angle;
 394 }
 395
 396 static char *get_atomtype_name_AB(t_atom *atom, gpp_atomtype *atype)
 397 {
 398     char *name;
 399
 400     name = get_atomtype_name(atom->type, atype);
 401
 402     /* When using the decoupling option, atom types are changed
 403      * to decoupled for the non-bonded interactions, but the virtual
 404      * sites constructions should be based on the "normal" interactions.
 405      * So we return the state B atom type name if the state A atom
 406      * type is the decoupled one. We should actually check for the atom
 407      * type number, but that's not passed here. So we check for
 408      * the decoupled atom type name. This should not cause trouble
 409      * as this code is only used for topologies with v-sites without
 410      * parameters generated by pdb2gmx.
 411      */
 412     if (strcmp(name, "decoupled") == 0)
 413     {
 414         name = get_atomtype_name(atom->typeB, atype);
 415     }
 416
 417     return name;
 418 }
 419
 420 static bool calc_vsite3_param(gpp_atomtype *atype,
 421                               t_param *param, t_atoms *at,
 422                               int nrbond, t_mybonded *bonds,
 423                               int nrang,  t_mybonded *angles )
 424 {
 425     /* i = virtual site          |    ,k
 426      * j = 1st bonded heavy atom | i-j
 427      * k,l = 2nd bonded atoms    |    `l
 428      */
 429
 430     bool     bXH3, bError;
 431     real     bjk, bjl, a = -1, b = -1;
 432     /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
 433      * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 434     bXH3 =
 435         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MNH", 3) == 0) &&
 436           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MNH", 3) == 0) ) ||
 437         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MCH3", 4) == 0) &&
 438           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MCH3", 4) == 0) );
 439
 440     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 441     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 442     bError = (bjk == NOTSET) || (bjl == NOTSET);
 443     if (bXH3)
 444     {
 445         /* now we get some XH2/XH3 group specific construction */
 446         /* note: we call the heavy atom 'C' and the X atom 'N' */
 447         real bMM, bCM, bCN, bNH, aCNH, dH, rH, dM, rM;
 448         int  aN;
 449
 450         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 451         bError = bError || (bjk != bjl);
 452
 453         /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
 454         aN = std::max(param->ak(), param->al())+1;
 455
 456         /* get common bonds */
 457         bMM    = get_bond_length(nrbond, bonds, param->ak(), param->al());
 458         bCM    = bjk;
 459         bCN    = get_bond_length(nrbond, bonds, param->aj(), aN);
 460         bError = bError || (bMM == NOTSET) || (bCN == NOTSET);
 461
 462         /* calculate common things */
 463         rM  = 0.5*bMM;
 464         dM  = std::sqrt( gmx::square(bCM) - gmx::square(rM) );
 465
 466         /* are we dealing with the X atom? */
 467         if (param->ai() == aN)
 468         {
 469             /* this is trivial */
 470             a = b = 0.5 * bCN/dM;
 471
 472         }
 473         else
 474         {
 475             /* get other bondlengths and angles: */
 476             bNH    = get_bond_length(nrbond, bonds, aN, param->ai());
 477             aCNH   = get_angle      (nrang, angles, param->aj(), aN, param->ai());
 478             bError = bError || (bNH == NOTSET) || (aCNH == NOTSET);
 479
 480             /* calculate */
 481             dH  = bCN - bNH * std::cos(aCNH);
 482             rH  = bNH * std::sin(aCNH);
 483
 484             a = 0.5 * ( dH/dM + rH/rM );
 485             b = 0.5 * ( dH/dM - rH/rM );
 486         }
 487     }
 488     else
 489     {
 490         gmx_fatal(FARGS, "calc_vsite3_param not implemented for the general case "
 491                   "(atom %d)", param->ai()+1);
 492     }
 493
 494     param->c0() = a;
 495     param->c1() = b;
 496
 497     return bError;
 498 }
 499
 500 static bool calc_vsite3fd_param(t_param *param,
 501                                 int nrbond, t_mybonded *bonds,
 502                                 int nrang,  t_mybonded *angles)
 503 {
 504     /* i = virtual site          |    ,k
 505      * j = 1st bonded heavy atom | i-j
 506      * k,l = 2nd bonded atoms    |    `l
 507      */
 508
 509     bool     bError;
 510     real     bij, bjk, bjl, aijk, aijl, rk, rl;
 511
 512     bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 513     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 514     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 515     aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 516     aijl   = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 517     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) ||
 518         (aijk == NOTSET) || (aijl == NOTSET);
 519
 520     rk          = bjk * std::sin(aijk);
 521     rl          = bjl * std::sin(aijl);
 522     param->c0() = rk / (rk + rl);
 523     param->c1() = -bij; /* 'bond'-length for fixed distance vsite */
 524
 525     return bError;
 526 }
 527
 528 static bool calc_vsite3fad_param(t_param *param,
 529                                  int nrbond, t_mybonded *bonds,
 530                                  int nrang,  t_mybonded *angles)
 531 {
 532     /* i = virtual site          |
 533      * j = 1st bonded heavy atom | i-j
 534      * k = 2nd bonded heavy atom |    `k-l
 535      * l = 3d bonded heavy atom  |
 536      */
 537
 538     bool     bSwapParity, bError;
 539     real     bij, aijk;
 540
 541     bSwapParity = ( param->c1() == -1 );
 542
 543     bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 544     aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 545     bError = (bij == NOTSET) || (aijk == NOTSET);
 546
 547     param->c1() = bij;          /* 'bond'-length for fixed distance vsite */
 548     param->c0() = RAD2DEG*aijk; /* 'bond'-angle for fixed angle vsite */
 549
 550     if (bSwapParity)
 551     {
 552         param->c0() = 360 - param->c0();
 553     }
 554
 555     return bError;
 556 }
 557
 558 static bool calc_vsite3out_param(gpp_atomtype *atype,
 559                                  t_param *param, t_atoms *at,
 560                                  int nrbond, t_mybonded *bonds,
 561                                  int nrang,  t_mybonded *angles)
 562 {
 563     /* i = virtual site          |    ,k
 564      * j = 1st bonded heavy atom | i-j
 565      * k,l = 2nd bonded atoms    |    `l
 566      * NOTE: i is out of the j-k-l plane!
 567      */
 568
 569     bool     bXH3, bError, bSwapParity;
 570     real     bij, bjk, bjl, aijk, aijl, akjl, pijk, pijl, a, b, c;
 571
 572     /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
 573      * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 574     bXH3 =
 575         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MNH", 3) == 0) &&
 576           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MNH", 3) == 0) ) ||
 577         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MCH3", 4) == 0) &&
 578           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MCH3", 4) == 0) );
 579
 580     /* check if construction parity must be swapped */
 581     bSwapParity = ( param->c1() == -1 );
 582
 583     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 584     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 585     bError = (bjk == NOTSET) || (bjl == NOTSET);
 586     if (bXH3)
 587     {
 588         /* now we get some XH3 group specific construction */
 589         /* note: we call the heavy atom 'C' and the X atom 'N' */
 590         real bMM, bCM, bCN, bNH, aCNH, dH, rH, rHx, rHy, dM, rM;
 591         int  aN;
 592
 593         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 594         bError = bError || (bjk != bjl);
 595
 596         /* the X atom (C or N) in the XH3 group is the first after the masses: */
 597         aN = std::max(param->ak(), param->al())+1;
 598
 599         /* get all bondlengths and angles: */
 600         bMM    = get_bond_length(nrbond, bonds, param->ak(), param->al());
 601         bCM    = bjk;
 602         bCN    = get_bond_length(nrbond, bonds, param->aj(), aN);
 603         bNH    = get_bond_length(nrbond, bonds, aN, param->ai());
 604         aCNH   = get_angle      (nrang, angles, param->aj(), aN, param->ai());
 605         bError = bError ||
 606             (bMM == NOTSET) || (bCN == NOTSET) || (bNH == NOTSET) || (aCNH == NOTSET);
 607
 608         /* calculate */
 609         dH  = bCN - bNH * std::cos(aCNH);
 610         rH  = bNH * std::sin(aCNH);
 611         /* we assume the H's are symmetrically distributed */
 612         rHx = rH*std::cos(DEG2RAD*30);
 613         rHy = rH*std::sin(DEG2RAD*30);
 614         rM  = 0.5*bMM;
 615         dM  = std::sqrt( gmx::square(bCM) - gmx::square(rM) );
 616         a   = 0.5*( (dH/dM) - (rHy/rM) );
 617         b   = 0.5*( (dH/dM) + (rHy/rM) );
 618         c   = rHx / (2*dM*rM);
 619
 620     }
 621     else
 622     {
 623         /* this is the general construction */
 624
 625         bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 626         aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 627         aijl   = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 628         akjl   = get_angle      (nrang, angles, param->ak(), param->aj(), param->al());
 629         bError = bError ||
 630             (bij == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET) || (akjl == NOTSET);
 631
 632         pijk = std::cos(aijk)*bij;
 633         pijl = std::cos(aijl)*bij;
 634         a    = ( pijk + (pijk*std::cos(akjl)-pijl) * std::cos(akjl) / gmx::square(std::sin(akjl)) ) / bjk;
 635         b    = ( pijl + (pijl*std::cos(akjl)-pijk) * std::cos(akjl) / gmx::square(std::sin(akjl)) ) / bjl;
 636         c    = -std::sqrt( gmx::square(bij) -
 637                            ( gmx::square(pijk) - 2*pijk*pijl*std::cos(akjl) + gmx::square(pijl) )
 638                            / gmx::square(std::sin(akjl)) )
 639             / ( bjk*bjl*std::sin(akjl) );
 640     }
 641
 642     param->c0() = a;
 643     param->c1() = b;
 644     if (bSwapParity)
 645     {
 646         param->c2() = -c;
 647     }
 648     else
 649     {
 650         param->c2() =  c;
 651     }
 652     return bError;
 653 }
 654
 655 static bool calc_vsite4fd_param(t_param *param,
 656                                 int nrbond, t_mybonded *bonds,
 657                                 int nrang,  t_mybonded *angles)
 658 {
 659     /* i = virtual site          |    ,k
 660      * j = 1st bonded heavy atom | i-j-m
 661      * k,l,m = 2nd bonded atoms  |    `l
 662      */
 663
 664     bool     bError;
 665     real     bij, bjk, bjl, bjm, aijk, aijl, aijm, akjm, akjl;
 666     real     pk, pl, pm, cosakl, cosakm, sinakl, sinakm, cl, cm;
 667
 668     bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 669     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 670     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 671     bjm    = get_bond_length(nrbond, bonds, param->aj(), param->am());
 672     aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 673     aijl   = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 674     aijm   = get_angle      (nrang, angles, param->ai(), param->aj(), param->am());
 675     akjm   = get_angle      (nrang, angles, param->ak(), param->aj(), param->am());
 676     akjl   = get_angle      (nrang, angles, param->ak(), param->aj(), param->al());
 677     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET) ||
 678         (aijk == NOTSET) || (aijl == NOTSET) || (aijm == NOTSET) || (akjm == NOTSET) ||
 679         (akjl == NOTSET);
 680
 681     if (!bError)
 682     {
 683         pk     = bjk*std::sin(aijk);
 684         pl     = bjl*std::sin(aijl);
 685         pm     = bjm*std::sin(aijm);
 686         cosakl = (std::cos(akjl) - std::cos(aijk)*std::cos(aijl)) / (std::sin(aijk)*std::sin(aijl));
 687         cosakm = (std::cos(akjm) - std::cos(aijk)*std::cos(aijm)) / (std::sin(aijk)*std::sin(aijm));
 688         if (cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1)
 689         {
 690             fprintf(stderr, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
 691                     param->ai()+1, RAD2DEG*aijk, RAD2DEG*aijl, RAD2DEG*aijm);
 692             gmx_fatal(FARGS, "invalid construction in calc_vsite4fd for atom %d: "
 693                       "cosakl=%f, cosakm=%f\n", param->ai()+1, cosakl, cosakm);
 694         }
 695         sinakl = std::sqrt(1-gmx::square(cosakl));
 696         sinakm = std::sqrt(1-gmx::square(cosakm));
 697
 698         /* note: there is a '+' because of the way the sines are calculated */
 699         cl = -pk / ( pl*cosakl - pk + pl*sinakl*(pm*cosakm-pk)/(pm*sinakm) );
 700         cm = -pk / ( pm*cosakm - pk + pm*sinakm*(pl*cosakl-pk)/(pl*sinakl) );
 701
 702         param->c0() = cl;
 703         param->c1() = cm;
 704         param->c2() = -bij;
 705     }
 706
 707     return bError;
 708 }
 709
 710
 711 static bool
 712 calc_vsite4fdn_param(t_param *param,
 713                      int nrbond, t_mybonded *bonds,
 714                      int nrang,  t_mybonded *angles)
 715 {
 716     /* i = virtual site          |    ,k
 717      * j = 1st bonded heavy atom | i-j-m
 718      * k,l,m = 2nd bonded atoms  |    `l
 719      */
 720
 721     bool     bError;
 722     real     bij, bjk, bjl, bjm, aijk, aijl, aijm;
 723     real     pk, pl, pm, a, b;
 724
 725     bij  = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 726     bjk  = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 727     bjl  = get_bond_length(nrbond, bonds, param->aj(), param->al());
 728     bjm  = get_bond_length(nrbond, bonds, param->aj(), param->am());
 729     aijk = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 730     aijl = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 731     aijm = get_angle      (nrang, angles, param->ai(), param->aj(), param->am());
 732
 733     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET) ||
 734         (aijk == NOTSET) || (aijl == NOTSET) || (aijm == NOTSET);
 735
 736     if (!bError)
 737     {
 738
 739         /* Calculate component of bond j-k along the direction i-j */
 740         pk = -bjk*std::cos(aijk);
 741
 742         /* Calculate component of bond j-l along the direction i-j */
 743         pl = -bjl*std::cos(aijl);
 744
 745         /* Calculate component of bond j-m along the direction i-j */
 746         pm = -bjm*std::cos(aijm);
 747
 748         if (fabs(pl) < 1000*GMX_REAL_MIN || fabs(pm) < 1000*GMX_REAL_MIN)
 749         {
 750             fprintf(stderr, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
 751                     param->ai()+1, RAD2DEG*aijk, RAD2DEG*aijl, RAD2DEG*aijm);
 752             gmx_fatal(FARGS, "invalid construction in calc_vsite4fdn for atom %d: "
 753                       "pl=%f, pm=%f\n", param->ai()+1, pl, pm);
 754         }
 755
 756         a = pk/pl;
 757         b = pk/pm;
 758
 759         param->c0() = a;
 760         param->c1() = b;
 761         param->c2() = bij;
 762
 763     }
 764
 765     return bError;
 766 }
 767
 768
 769
 770 int set_vsites(bool bVerbose, t_atoms *atoms, gpp_atomtype *atype,
 771                gmx::ArrayRef<InteractionTypeParameters> plist)
 772 {
 773     int             i, j, ftype;
 774     int             nvsite, nrbond, nrang, nridih, nrset;
 775     bool            bFirst, bSet, bERROR;
 776     at2vsitebond_t *at2vb;
 777     t_mybonded     *bonds;
 778     t_mybonded     *angles;
 779     t_mybonded     *idihs;
 780
 781     bFirst = TRUE;
 782     nvsite = 0;
 783
 784     /* Make a reverse list to avoid ninteractions^2 operations */
 785     at2vb = make_at2vsitebond(atoms->nr, plist);
 786
 787     for (ftype = 0; (ftype < F_NRE); ftype++)
 788     {
 789         if (interaction_function[ftype].flags & IF_VSITE)
 790         {
 791             nvsite += plist[ftype].nr;
 792
 793             if (ftype == F_VSITEN)
 794             {
 795                 /* We don't calculate parameters for VSITEN */
 796                 continue;
 797             }
 798
 799             nrset = 0;
 800             for (i = 0; (i < plist[ftype].nr); i++)
 801             {
 802                 /* check if all parameters are set */
 803                 bSet = TRUE;
 804                 for (j = 0; j < NRFP(ftype) && bSet; j++)
 805                 {
 806                     bSet = plist[ftype].param[i].c[j] != NOTSET;
 807                 }
 808
 809                 if (debug)
 810                 {
 811                     fprintf(debug, "bSet=%s ", gmx::boolToString(bSet));
 812                     print_param(debug, ftype, i, &plist[ftype].param[i]);
 813                 }
 814                 if (!bSet)
 815                 {
 816                     if (bVerbose && bFirst)
 817                     {
 818                         fprintf(stderr, "Calculating parameters for virtual sites\n");
 819                         bFirst = FALSE;
 820                     }
 821
 822                     nrbond = nrang = nridih = 0;
 823                     bonds  = nullptr;
 824                     angles = nullptr;
 825                     idihs  = nullptr;
 826                     nrset++;
 827                     /* now set the vsite parameters: */
 828                     get_bondeds(NRAL(ftype), plist[ftype].param[i].a, at2vb,
 829                                 &nrbond, &bonds, &nrang,  &angles, &nridih, &idihs);
 830                     if (debug)
 831                     {
 832                         fprintf(debug, "Found %d bonds, %d angles and %d idihs "
 833                                 "for virtual site %d (%s)\n", nrbond, nrang, nridih,
 834                                 plist[ftype].param[i].ai()+1,
 835                                 interaction_function[ftype].longname);
 836                         print_bad(debug, nrbond, bonds, nrang, angles, nridih, idihs);
 837                     } /* debug */
 838                     switch (ftype)
 839                     {
 840                         case F_VSITE3:
 841                             bERROR =
 842                                 calc_vsite3_param(atype, &(plist[ftype].param[i]), atoms,
 843                                                   nrbond, bonds, nrang, angles);
 844                             break;
 845                         case F_VSITE3FD:
 846                             bERROR =
 847                                 calc_vsite3fd_param(&(plist[ftype].param[i]),
 848                                                     nrbond, bonds, nrang, angles);
 849                             break;
 850                         case F_VSITE3FAD:
 851                             bERROR =
 852                                 calc_vsite3fad_param(&(plist[ftype].param[i]),
 853                                                      nrbond, bonds, nrang, angles);
 854                             break;
 855                         case F_VSITE3OUT:
 856                             bERROR =
 857                                 calc_vsite3out_param(atype, &(plist[ftype].param[i]), atoms,
 858                                                      nrbond, bonds, nrang, angles);
 859                             break;
 860                         case F_VSITE4FD:
 861                             bERROR =
 862                                 calc_vsite4fd_param(&(plist[ftype].param[i]),
 863                                                     nrbond, bonds, nrang, angles);
 864                             break;
 865                         case F_VSITE4FDN:
 866                             bERROR =
 867                                 calc_vsite4fdn_param(&(plist[ftype].param[i]),
 868                                                      nrbond, bonds, nrang, angles);
 869                             break;
 870                         default:
 871                             gmx_fatal(FARGS, "Automatic parameter generation not supported "
 872                                       "for %s atom %d",
 873                                       interaction_function[ftype].longname,
 874                                       plist[ftype].param[i].ai()+1);
 875                             bERROR = TRUE;
 876                     } /* switch */
 877                     if (bERROR)
 878                     {
 879                         gmx_fatal(FARGS, "Automatic parameter generation not supported "
 880                                   "for %s atom %d for this bonding configuration",
 881                                   interaction_function[ftype].longname,
 882                                   plist[ftype].param[i].ai()+1);
 883                     }
 884                     sfree(bonds);
 885                     sfree(angles);
 886                     sfree(idihs);
 887                 } /* if bSet */
 888             }     /* for i */
 889         }         /* if IF_VSITE */
 890
 891     }
 892     done_at2vsitebond(atoms->nr, at2vb);
 893
 894     return nvsite;
 895 }
 896
 897 void set_vsites_ptype(bool bVerbose, gmx_moltype_t *molt)
 898 {
 899     int      ftype, i;
 900
 901     if (bVerbose)
 902     {
 903         fprintf(stderr, "Setting particle type to V for virtual sites\n");
 904     }
 905     for (ftype = 0; ftype < F_NRE; ftype++)
 906     {
 907         InteractionList *il = &molt->ilist[ftype];
 908         if (interaction_function[ftype].flags & IF_VSITE)
 909         {
 910             const int                nra = interaction_function[ftype].nratoms;
 911             const int                nrd = il->size();
 912             gmx::ArrayRef<const int> ia  = il->iatoms;
 913
 914             if (debug && nrd)
 915             {
 916                 fprintf(stderr, "doing %d %s virtual sites\n",
 917                         (nrd / (nra+1)), interaction_function[ftype].longname);
 918             }
 919
 920             for (i = 0; (i < nrd); )
 921             {
 922                 /* The virtual site */
 923                 int avsite = ia[i + 1];
 924                 molt->atoms.atom[avsite].ptype = eptVSite;
 925
 926                 i  += 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(gmx::ArrayRef<InteractionTypeParameters> plist,
 938                                     int cftype, const int vsite_type[])
 939 {
 940     int                        i, k, n;
 941     int                        atom;
 942     InteractionTypeParameters *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(gmx::ArrayRef<InteractionTypeParameters> 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     InteractionTypeParameters    *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(gmx::ArrayRef<InteractionTypeParameters> 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     InteractionTypeParameters    *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(gmx::ArrayRef<InteractionTypeParameters> plist, t_pindex pindex[],
1327                              int cftype, const int vsite_type[])
1328 {
1329     int                        i, kept_i;
1330     InteractionTypeParameters *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(gmx::ArrayRef<InteractionTypeParameters> 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 }