src/gromacs/gmxpreprocess/vsite_parm.cpp

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