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 <stdio.h>
  42 #include <string.h>
  43
  44 #include <cmath>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/gmxpreprocess/add_par.h"
  49 #include "gromacs/gmxpreprocess/notset.h"
  50 #include "gromacs/gmxpreprocess/resall.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/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, 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     gmx_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     gmx_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 gmx_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     gmx_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     if (debug)
 432     {
 433         int i;
 434         for (i = 0; i < 4; i++)
 435         {
 436             fprintf(debug, "atom %d type %s ",
 437                     param->a[i]+1,
 438                     get_atomtype_name_AB(&at->atom[param->a[i]], atype));
 439         }
 440         fprintf(debug, "\n");
 441     }
 442     bXH3 =
 443         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MNH", 3) == 0) &&
 444           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MNH", 3) == 0) ) ||
 445         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MCH3", 4) == 0) &&
 446           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MCH3", 4) == 0) );
 447
 448     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 449     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 450     bError = (bjk == NOTSET) || (bjl == NOTSET);
 451     if (bXH3)
 452     {
 453         /* now we get some XH2/XH3 group specific construction */
 454         /* note: we call the heavy atom 'C' and the X atom 'N' */
 455         real bMM, bCM, bCN, bNH, aCNH, dH, rH, dM, rM;
 456         int  aN;
 457
 458         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 459         bError = bError || (bjk != bjl);
 460
 461         /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
 462         aN = std::max(param->ak(), param->al())+1;
 463
 464         /* get common bonds */
 465         bMM    = get_bond_length(nrbond, bonds, param->ak(), param->al());
 466         bCM    = bjk;
 467         bCN    = get_bond_length(nrbond, bonds, param->aj(), aN);
 468         bError = bError || (bMM == NOTSET) || (bCN == NOTSET);
 469
 470         /* calculate common things */
 471         rM  = 0.5*bMM;
 472         dM  = std::sqrt( gmx::square(bCM) - gmx::square(rM) );
 473
 474         /* are we dealing with the X atom? */
 475         if (param->ai() == aN)
 476         {
 477             /* this is trivial */
 478             a = b = 0.5 * bCN/dM;
 479
 480         }
 481         else
 482         {
 483             /* get other bondlengths and angles: */
 484             bNH    = get_bond_length(nrbond, bonds, aN, param->ai());
 485             aCNH   = get_angle      (nrang, angles, param->aj(), aN, param->ai());
 486             bError = bError || (bNH == NOTSET) || (aCNH == NOTSET);
 487
 488             /* calculate */
 489             dH  = bCN - bNH * std::cos(aCNH);
 490             rH  = bNH * std::sin(aCNH);
 491
 492             a = 0.5 * ( dH/dM + rH/rM );
 493             b = 0.5 * ( dH/dM - rH/rM );
 494         }
 495     }
 496     else
 497     {
 498         gmx_fatal(FARGS, "calc_vsite3_param not implemented for the general case "
 499                   "(atom %d)", param->ai()+1);
 500     }
 501
 502     param->c0() = a;
 503     param->c1() = b;
 504
 505     if (debug)
 506     {
 507         fprintf(debug, "params for vsite3 %d: %g %g\n",
 508                 param->ai()+1, param->c0(), param->c1());
 509     }
 510
 511     return bError;
 512 }
 513
 514 static gmx_bool calc_vsite3fd_param(t_param *param,
 515                                     int nrbond, t_mybonded *bonds,
 516                                     int nrang,  t_mybonded *angles)
 517 {
 518     /* i = virtual site          |    ,k
 519      * j = 1st bonded heavy atom | i-j
 520      * k,l = 2nd bonded atoms    |    `l
 521      */
 522
 523     gmx_bool bError;
 524     real     bij, bjk, bjl, aijk, aijl, rk, rl;
 525
 526     bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 527     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 528     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 529     aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 530     aijl   = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 531     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) ||
 532         (aijk == NOTSET) || (aijl == NOTSET);
 533
 534     rk          = bjk * std::sin(aijk);
 535     rl          = bjl * std::sin(aijl);
 536     param->c0() = rk / (rk + rl);
 537     param->c1() = -bij; /* 'bond'-length for fixed distance vsite */
 538
 539     if (debug)
 540     {
 541         fprintf(debug, "params for vsite3fd %d: %g %g\n",
 542                 param->ai()+1, param->c0(), param->c1());
 543     }
 544     return bError;
 545 }
 546
 547 static gmx_bool calc_vsite3fad_param(t_param *param,
 548                                      int nrbond, t_mybonded *bonds,
 549                                      int nrang,  t_mybonded *angles)
 550 {
 551     /* i = virtual site          |
 552      * j = 1st bonded heavy atom | i-j
 553      * k = 2nd bonded heavy atom |    `k-l
 554      * l = 3d bonded heavy atom  |
 555      */
 556
 557     gmx_bool bSwapParity, bError;
 558     real     bij, aijk;
 559
 560     bSwapParity = ( param->c1() == -1 );
 561
 562     bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 563     aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 564     bError = (bij == NOTSET) || (aijk == NOTSET);
 565
 566     param->c1() = bij;          /* 'bond'-length for fixed distance vsite */
 567     param->c0() = RAD2DEG*aijk; /* 'bond'-angle for fixed angle vsite */
 568
 569     if (bSwapParity)
 570     {
 571         param->c0() = 360 - param->c0();
 572     }
 573
 574     if (debug)
 575     {
 576         fprintf(debug, "params for vsite3fad %d: %g %g\n",
 577                 param->ai()+1, param->c0(), param->c1());
 578     }
 579     return bError;
 580 }
 581
 582 static gmx_bool calc_vsite3out_param(gpp_atomtype_t atype,
 583                                      t_param *param, t_atoms *at,
 584                                      int nrbond, t_mybonded *bonds,
 585                                      int nrang,  t_mybonded *angles)
 586 {
 587     /* i = virtual site          |    ,k
 588      * j = 1st bonded heavy atom | i-j
 589      * k,l = 2nd bonded atoms    |    `l
 590      * NOTE: i is out of the j-k-l plane!
 591      */
 592
 593     gmx_bool bXH3, bError, bSwapParity;
 594     real     bij, bjk, bjl, aijk, aijl, akjl, pijk, pijl, a, b, c;
 595
 596     /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
 597      * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 598     if (debug)
 599     {
 600         int i;
 601         for (i = 0; i < 4; i++)
 602         {
 603             fprintf(debug, "atom %d type %s ",
 604                     param->a[i]+1, get_atomtype_name_AB(&at->atom[param->a[i]], atype));
 605         }
 606         fprintf(debug, "\n");
 607     }
 608     bXH3 =
 609         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MNH", 3) == 0) &&
 610           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MNH", 3) == 0) ) ||
 611         ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MCH3", 4) == 0) &&
 612           (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MCH3", 4) == 0) );
 613
 614     /* check if construction parity must be swapped */
 615     bSwapParity = ( param->c1() == -1 );
 616
 617     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 618     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 619     bError = (bjk == NOTSET) || (bjl == NOTSET);
 620     if (bXH3)
 621     {
 622         /* now we get some XH3 group specific construction */
 623         /* note: we call the heavy atom 'C' and the X atom 'N' */
 624         real bMM, bCM, bCN, bNH, aCNH, dH, rH, rHx, rHy, dM, rM;
 625         int  aN;
 626
 627         /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 628         bError = bError || (bjk != bjl);
 629
 630         /* the X atom (C or N) in the XH3 group is the first after the masses: */
 631         aN = std::max(param->ak(), param->al())+1;
 632
 633         /* get all bondlengths and angles: */
 634         bMM    = get_bond_length(nrbond, bonds, param->ak(), param->al());
 635         bCM    = bjk;
 636         bCN    = get_bond_length(nrbond, bonds, param->aj(), aN);
 637         bNH    = get_bond_length(nrbond, bonds, aN, param->ai());
 638         aCNH   = get_angle      (nrang, angles, param->aj(), aN, param->ai());
 639         bError = bError ||
 640             (bMM == NOTSET) || (bCN == NOTSET) || (bNH == NOTSET) || (aCNH == NOTSET);
 641
 642         /* calculate */
 643         dH  = bCN - bNH * std::cos(aCNH);
 644         rH  = bNH * std::sin(aCNH);
 645         /* we assume the H's are symmetrically distributed */
 646         rHx = rH*std::cos(DEG2RAD*30);
 647         rHy = rH*std::sin(DEG2RAD*30);
 648         rM  = 0.5*bMM;
 649         dM  = std::sqrt( gmx::square(bCM) - gmx::square(rM) );
 650         a   = 0.5*( (dH/dM) - (rHy/rM) );
 651         b   = 0.5*( (dH/dM) + (rHy/rM) );
 652         c   = rHx / (2*dM*rM);
 653
 654     }
 655     else
 656     {
 657         /* this is the general construction */
 658
 659         bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 660         aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 661         aijl   = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 662         akjl   = get_angle      (nrang, angles, param->ak(), param->aj(), param->al());
 663         bError = bError ||
 664             (bij == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET) || (akjl == NOTSET);
 665
 666         pijk = std::cos(aijk)*bij;
 667         pijl = std::cos(aijl)*bij;
 668         a    = ( pijk + (pijk*std::cos(akjl)-pijl) * std::cos(akjl) / gmx::square(std::sin(akjl)) ) / bjk;
 669         b    = ( pijl + (pijl*std::cos(akjl)-pijk) * std::cos(akjl) / gmx::square(std::sin(akjl)) ) / bjl;
 670         c    = -std::sqrt( gmx::square(bij) -
 671                            ( gmx::square(pijk) - 2*pijk*pijl*std::cos(akjl) + gmx::square(pijl) )
 672                            / gmx::square(std::sin(akjl)) )
 673             / ( bjk*bjl*std::sin(akjl) );
 674     }
 675
 676     param->c0() = a;
 677     param->c1() = b;
 678     if (bSwapParity)
 679     {
 680         param->c2() = -c;
 681     }
 682     else
 683     {
 684         param->c2() =  c;
 685     }
 686     if (debug)
 687     {
 688         fprintf(debug, "params for vsite3out %d: %g %g %g\n",
 689                 param->ai()+1, param->c0(), param->c1(), param->c2());
 690     }
 691     return bError;
 692 }
 693
 694 static gmx_bool calc_vsite4fd_param(t_param *param,
 695                                     int nrbond, t_mybonded *bonds,
 696                                     int nrang,  t_mybonded *angles)
 697 {
 698     /* i = virtual site          |    ,k
 699      * j = 1st bonded heavy atom | i-j-m
 700      * k,l,m = 2nd bonded atoms  |    `l
 701      */
 702
 703     gmx_bool bError;
 704     real     bij, bjk, bjl, bjm, aijk, aijl, aijm, akjm, akjl;
 705     real     pk, pl, pm, cosakl, cosakm, sinakl, sinakm, cl, cm;
 706
 707     bij    = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 708     bjk    = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 709     bjl    = get_bond_length(nrbond, bonds, param->aj(), param->al());
 710     bjm    = get_bond_length(nrbond, bonds, param->aj(), param->am());
 711     aijk   = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 712     aijl   = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 713     aijm   = get_angle      (nrang, angles, param->ai(), param->aj(), param->am());
 714     akjm   = get_angle      (nrang, angles, param->ak(), param->aj(), param->am());
 715     akjl   = get_angle      (nrang, angles, param->ak(), param->aj(), param->al());
 716     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET) ||
 717         (aijk == NOTSET) || (aijl == NOTSET) || (aijm == NOTSET) || (akjm == NOTSET) ||
 718         (akjl == NOTSET);
 719
 720     if (!bError)
 721     {
 722         pk     = bjk*std::sin(aijk);
 723         pl     = bjl*std::sin(aijl);
 724         pm     = bjm*std::sin(aijm);
 725         cosakl = (std::cos(akjl) - std::cos(aijk)*std::cos(aijl)) / (std::sin(aijk)*std::sin(aijl));
 726         cosakm = (std::cos(akjm) - std::cos(aijk)*std::cos(aijm)) / (std::sin(aijk)*std::sin(aijm));
 727         if (cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1)
 728         {
 729             fprintf(stderr, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
 730                     param->ai()+1, RAD2DEG*aijk, RAD2DEG*aijl, RAD2DEG*aijm);
 731             gmx_fatal(FARGS, "invalid construction in calc_vsite4fd for atom %d: "
 732                       "cosakl=%f, cosakm=%f\n", param->ai()+1, cosakl, cosakm);
 733         }
 734         sinakl = std::sqrt(1-gmx::square(cosakl));
 735         sinakm = std::sqrt(1-gmx::square(cosakm));
 736
 737         /* note: there is a '+' because of the way the sines are calculated */
 738         cl = -pk / ( pl*cosakl - pk + pl*sinakl*(pm*cosakm-pk)/(pm*sinakm) );
 739         cm = -pk / ( pm*cosakm - pk + pm*sinakm*(pl*cosakl-pk)/(pl*sinakl) );
 740
 741         param->c0() = cl;
 742         param->c1() = cm;
 743         param->c2() = -bij;
 744         if (debug)
 745         {
 746             fprintf(debug, "params for vsite4fd %d: %g %g %g\n",
 747                     param->ai()+1, param->c0(), param->c1(), param->c2());
 748         }
 749     }
 750
 751     return bError;
 752 }
 753
 754
 755 static gmx_bool
 756 calc_vsite4fdn_param(t_param *param,
 757                      int nrbond, t_mybonded *bonds,
 758                      int nrang,  t_mybonded *angles)
 759 {
 760     /* i = virtual site          |    ,k
 761      * j = 1st bonded heavy atom | i-j-m
 762      * k,l,m = 2nd bonded atoms  |    `l
 763      */
 764
 765     gmx_bool bError;
 766     real     bij, bjk, bjl, bjm, aijk, aijl, aijm;
 767     real     pk, pl, pm, a, b;
 768
 769     bij  = get_bond_length(nrbond, bonds, param->ai(), param->aj());
 770     bjk  = get_bond_length(nrbond, bonds, param->aj(), param->ak());
 771     bjl  = get_bond_length(nrbond, bonds, param->aj(), param->al());
 772     bjm  = get_bond_length(nrbond, bonds, param->aj(), param->am());
 773     aijk = get_angle      (nrang, angles, param->ai(), param->aj(), param->ak());
 774     aijl = get_angle      (nrang, angles, param->ai(), param->aj(), param->al());
 775     aijm = get_angle      (nrang, angles, param->ai(), param->aj(), param->am());
 776
 777     bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET) ||
 778         (aijk == NOTSET) || (aijl == NOTSET) || (aijm == NOTSET);
 779
 780     if (!bError)
 781     {
 782
 783         /* Calculate component of bond j-k along the direction i-j */
 784         pk = -bjk*std::cos(aijk);
 785
 786         /* Calculate component of bond j-l along the direction i-j */
 787         pl = -bjl*std::cos(aijl);
 788
 789         /* Calculate component of bond j-m along the direction i-j */
 790         pm = -bjm*std::cos(aijm);
 791
 792         if (fabs(pl) < 1000*GMX_REAL_MIN || fabs(pm) < 1000*GMX_REAL_MIN)
 793         {
 794             fprintf(stderr, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
 795                     param->ai()+1, RAD2DEG*aijk, RAD2DEG*aijl, RAD2DEG*aijm);
 796             gmx_fatal(FARGS, "invalid construction in calc_vsite4fdn for atom %d: "
 797                       "pl=%f, pm=%f\n", param->ai()+1, pl, pm);
 798         }
 799
 800         a = pk/pl;
 801         b = pk/pm;
 802
 803         param->c0() = a;
 804         param->c1() = b;
 805         param->c2() = bij;
 806
 807         if (debug)
 808         {
 809             fprintf(debug, "params for vsite4fdn %d: %g %g %g\n",
 810                     param->ai()+1, param->c0(), param->c1(), param->c2());
 811         }
 812     }
 813
 814     return bError;
 815 }
 816
 817
 818
 819 int set_vsites(gmx_bool bVerbose, t_atoms *atoms, gpp_atomtype_t atype,
 820                t_params plist[])
 821 {
 822     int             i, j, ftype;
 823     int             nvsite, nrbond, nrang, nridih, nrset;
 824     gmx_bool        bFirst, bSet, bERROR;
 825     at2vsitebond_t *at2vb;
 826     t_mybonded     *bonds;
 827     t_mybonded     *angles;
 828     t_mybonded     *idihs;
 829
 830     bFirst = TRUE;
 831     nvsite = 0;
 832     if (debug)
 833     {
 834         fprintf(debug, "\nCalculating parameters for virtual sites\n");
 835     }
 836
 837     /* Make a reverse list to avoid ninteractions^2 operations */
 838     at2vb = make_at2vsitebond(atoms->nr, plist);
 839
 840     for (ftype = 0; (ftype < F_NRE); ftype++)
 841     {
 842         if (interaction_function[ftype].flags & IF_VSITE)
 843         {
 844             nvsite += plist[ftype].nr;
 845
 846             if (ftype == F_VSITEN)
 847             {
 848                 /* We don't calculate parameters for VSITEN */
 849                 continue;
 850             }
 851
 852             nrset = 0;
 853             for (i = 0; (i < plist[ftype].nr); i++)
 854             {
 855                 /* check if all parameters are set */
 856                 bSet = TRUE;
 857                 for (j = 0; j < NRFP(ftype) && bSet; j++)
 858                 {
 859                     bSet = plist[ftype].param[i].c[j] != NOTSET;
 860                 }
 861
 862                 if (debug)
 863                 {
 864                     fprintf(debug, "bSet=%s ", gmx::boolToString(bSet));
 865                     print_param(debug, ftype, i, &plist[ftype].param[i]);
 866                 }
 867                 if (!bSet)
 868                 {
 869                     if (bVerbose && bFirst)
 870                     {
 871                         fprintf(stderr, "Calculating parameters for virtual sites\n");
 872                         bFirst = FALSE;
 873                     }
 874
 875                     nrbond = nrang = nridih = 0;
 876                     bonds  = nullptr;
 877                     angles = nullptr;
 878                     idihs  = nullptr;
 879                     nrset++;
 880                     /* now set the vsite parameters: */
 881                     get_bondeds(NRAL(ftype), plist[ftype].param[i].a, at2vb,
 882                                 &nrbond, &bonds, &nrang,  &angles, &nridih, &idihs);
 883                     if (debug)
 884                     {
 885                         fprintf(debug, "Found %d bonds, %d angles and %d idihs "
 886                                 "for virtual site %d (%s)\n", nrbond, nrang, nridih,
 887                                 plist[ftype].param[i].ai()+1,
 888                                 interaction_function[ftype].longname);
 889                         print_bad(debug, nrbond, bonds, nrang, angles, nridih, idihs);
 890                     } /* debug */
 891                     switch (ftype)
 892                     {
 893                         case F_VSITE3:
 894                             bERROR =
 895                                 calc_vsite3_param(atype, &(plist[ftype].param[i]), atoms,
 896                                                   nrbond, bonds, nrang, angles);
 897                             break;
 898                         case F_VSITE3FD:
 899                             bERROR =
 900                                 calc_vsite3fd_param(&(plist[ftype].param[i]),
 901                                                     nrbond, bonds, nrang, angles);
 902                             break;
 903                         case F_VSITE3FAD:
 904                             bERROR =
 905                                 calc_vsite3fad_param(&(plist[ftype].param[i]),
 906                                                      nrbond, bonds, nrang, angles);
 907                             break;
 908                         case F_VSITE3OUT:
 909                             bERROR =
 910                                 calc_vsite3out_param(atype, &(plist[ftype].param[i]), atoms,
 911                                                      nrbond, bonds, nrang, angles);
 912                             break;
 913                         case F_VSITE4FD:
 914                             bERROR =
 915                                 calc_vsite4fd_param(&(plist[ftype].param[i]),
 916                                                     nrbond, bonds, nrang, angles);
 917                             break;
 918                         case F_VSITE4FDN:
 919                             bERROR =
 920                                 calc_vsite4fdn_param(&(plist[ftype].param[i]),
 921                                                      nrbond, bonds, nrang, angles);
 922                             break;
 923                         default:
 924                             gmx_fatal(FARGS, "Automatic parameter generation not supported "
 925                                       "for %s atom %d",
 926                                       interaction_function[ftype].longname,
 927                                       plist[ftype].param[i].ai()+1);
 928                             bERROR = TRUE;
 929                     } /* switch */
 930                     if (bERROR)
 931                     {
 932                         gmx_fatal(FARGS, "Automatic parameter generation not supported "
 933                                   "for %s atom %d for this bonding configuration",
 934                                   interaction_function[ftype].longname,
 935                                   plist[ftype].param[i].ai()+1);
 936                     }
 937                     sfree(bonds);
 938                     sfree(angles);
 939                     sfree(idihs);
 940                 } /* if bSet */
 941             }     /* for i */
 942             if (debug && plist[ftype].nr)
 943             {
 944                 fprintf(stderr, "Calculated parameters for %d out of %d %s atoms\n",
 945                         nrset, plist[ftype].nr, interaction_function[ftype].longname);
 946             }
 947         } /* if IF_VSITE */
 948
 949     }
 950     done_at2vsitebond(atoms->nr, at2vb);
 951
 952     return nvsite;
 953 }
 954
 955 void set_vsites_ptype(gmx_bool bVerbose, gmx_moltype_t *molt)
 956 {
 957     int      ftype, i;
 958     int      nra, nrd;
 959     t_ilist *il;
 960     t_iatom *ia, avsite;
 961
 962     if (bVerbose)
 963     {
 964         fprintf(stderr, "Setting particle type to V for virtual sites\n");
 965     }
 966     if (debug)
 967     {
 968         fprintf(stderr, "checking %d functypes\n", F_NRE);
 969     }
 970     for (ftype = 0; ftype < F_NRE; ftype++)
 971     {
 972         il = &molt->ilist[ftype];
 973         if (interaction_function[ftype].flags & IF_VSITE)
 974         {
 975             nra    = interaction_function[ftype].nratoms;
 976             nrd    = il->nr;
 977             ia     = il->iatoms;
 978
 979             if (debug && nrd)
 980             {
 981                 fprintf(stderr, "doing %d %s virtual sites\n",
 982                         (nrd / (nra+1)), interaction_function[ftype].longname);
 983             }
 984
 985             for (i = 0; (i < nrd); )
 986             {
 987                 /* The virtual site */
 988                 avsite = ia[1];
 989                 molt->atoms.atom[avsite].ptype = eptVSite;
 990
 991                 i  += nra+1;
 992                 ia += nra+1;
 993             }
 994         }
 995     }
 996
 997 }
 998
 999 typedef struct {
1000     int ftype, parnr;
1001 } t_pindex;
1002
1003 static void check_vsite_constraints(t_params *plist,
1004                                     int cftype, int vsite_type[])
1005 {
1006     int       i, k, n;
1007     int       atom;
1008     t_params *ps;
1009
1010     n  = 0;
1011     ps = &(plist[cftype]);
1012     for (i = 0; (i < ps->nr); i++)
1013     {
1014         for (k = 0; k < 2; k++)
1015         {
1016             atom = ps->param[i].a[k];
1017             if (vsite_type[atom] != NOTSET)
1018             {
1019                 fprintf(stderr, "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)\n",
1020                         ps->param[i].ai()+1, ps->param[i].aj()+1, atom+1);
1021                 n++;
1022             }
1023         }
1024     }
1025     if (n)
1026     {
1027         gmx_fatal(FARGS, "There were %d virtual sites involved in constraints", n);
1028     }
1029 }
1030
1031 static void clean_vsite_bonds(t_params *plist, t_pindex pindex[],
1032                               int cftype, int vsite_type[])
1033 {
1034     int          ftype, i, j, k, m, n, nvsite, nOut, kept_i;
1035     int          nconverted, nremoved;
1036     int          atom, oatom, at1, at2;
1037     gmx_bool     bKeep, bRemove, bUsed, bPresent, bThisFD, bThisOUT, bAllFD, bFirstTwo;
1038     t_params    *ps;
1039
1040     if (cftype == F_CONNBONDS)
1041     {
1042         return;
1043     }
1044
1045     ps         = &(plist[cftype]);
1046     kept_i     = 0;
1047     nconverted = 0;
1048     nremoved   = 0;
1049     nOut       = 0;
1050     for (i = 0; (i < ps->nr); i++) /* for all bonds in the plist */
1051     {
1052         int            vsnral      = 0;
1053         const int     *first_atoms = nullptr;
1054
1055         bKeep   = FALSE;
1056         bRemove = FALSE;
1057         bAllFD  = TRUE;
1058         /* check if all virtual sites are constructed from the same atoms */
1059         nvsite = 0;
1060         if (debug)
1061         {
1062             fprintf(debug, "constr %d %d:", ps->param[i].ai()+1, ps->param[i].aj()+1);
1063         }
1064         for (k = 0; (k < 2) && !bKeep && !bRemove; k++)
1065         {
1066             /* for all atoms in the bond */
1067             atom = ps->param[i].a[k];
1068             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1069             {
1070                 nvsite++;
1071                 bThisFD = ( (pindex[atom].ftype == F_VSITE3FD ) ||
1072                             (pindex[atom].ftype == F_VSITE3FAD) ||
1073                             (pindex[atom].ftype == F_VSITE4FD ) ||
1074                             (pindex[atom].ftype == F_VSITE4FDN ) );
1075                 bThisOUT = ( (pindex[atom].ftype == F_VSITE3OUT) &&
1076                              (interaction_function[cftype].flags & IF_CONSTRAINT) );
1077                 bAllFD = bAllFD && bThisFD;
1078                 if (bThisFD || bThisOUT)
1079                 {
1080                     if (debug)
1081                     {
1082                         fprintf(debug, " %s", bThisOUT ? "out" : "fd");
1083                     }
1084                     oatom = ps->param[i].a[1-k]; /* the other atom */
1085                     if (vsite_type[oatom] == NOTSET &&
1086                         oatom == plist[pindex[atom].ftype].param[pindex[atom].parnr].aj())
1087                     {
1088                         /* if the other atom isn't a vsite, and it is AI */
1089                         bRemove = TRUE;
1090                         if (bThisOUT)
1091                         {
1092                             nOut++;
1093                         }
1094                         if (debug)
1095                         {
1096                             fprintf(debug, " D-AI");
1097                         }
1098                     }
1099                 }
1100                 if (!bRemove)
1101                 {
1102                     /* TODO This fragment, and corresponding logic in
1103                        clean_vsite_angles and clean_vsite_dihs should
1104                        be refactored into a common function */
1105                     if (nvsite == 1)
1106                     {
1107                         /* if this is the first vsite we encounter then
1108                            store construction atoms */
1109                         /* TODO This would be nicer to implement with
1110                            a C++ "vector view" class" with an
1111                            STL-container-like interface. */
1112                         vsnral      = NRAL(pindex[atom].ftype) - 1;
1113                         first_atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1;
1114                     }
1115                     else
1116                     {
1117                         GMX_ASSERT(vsnral != 0, "nvsite > 1 must have vsnral != 0");
1118                         GMX_ASSERT(first_atoms != NULL, "nvsite > 1 must have first_atoms != NULL");
1119                         /* if it is not the first then
1120                            check if this vsite is constructed from the same atoms */
1121                         if (vsnral == NRAL(pindex[atom].ftype)-1)
1122                         {
1123                             for (m = 0; (m < vsnral) && !bKeep; m++)
1124                             {
1125                                 const int *atoms;
1126
1127                                 bPresent = FALSE;
1128                                 atoms    = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1;
1129                                 for (n = 0; (n < vsnral) && !bPresent; n++)
1130                                 {
1131                                     if (atoms[m] == first_atoms[n])
1132                                     {
1133                                         bPresent = TRUE;
1134                                     }
1135                                 }
1136                                 if (!bPresent)
1137                                 {
1138                                     bKeep = TRUE;
1139                                     if (debug)
1140                                     {
1141                                         fprintf(debug, " !present");
1142                                     }
1143                                 }
1144                             }
1145                         }
1146                         else
1147                         {
1148                             bKeep = TRUE;
1149                             if (debug)
1150                             {
1151                                 fprintf(debug, " !same#at");
1152                             }
1153                         }
1154                     }
1155                 }
1156             }
1157         }
1158
1159         if (bRemove)
1160         {
1161             bKeep = FALSE;
1162         }
1163         else
1164         {
1165             /* if we have no virtual sites in this bond, keep it */
1166             if (nvsite == 0)
1167             {
1168                 if (debug)
1169                 {
1170                     fprintf(debug, " no vsite");
1171                 }
1172                 bKeep = TRUE;
1173             }
1174
1175             /* TODO This loop and the corresponding loop in
1176                check_vsite_angles should be refactored into a common
1177                function */
1178             /* check if all non-vsite atoms are used in construction: */
1179             bFirstTwo = TRUE;
1180             for (k = 0; (k < 2) && !bKeep; k++) /* for all atoms in the bond */
1181             {
1182                 atom = ps->param[i].a[k];
1183                 if (vsite_type[atom] == NOTSET)
1184                 {
1185                     bUsed = FALSE;
1186                     for (m = 0; (m < vsnral) && !bUsed; m++)
1187                     {
1188                         GMX_ASSERT(first_atoms != NULL, "If we've seen a vsite before, we know what its first atom index was");
1189
1190                         if (atom == first_atoms[m])
1191                         {
1192                             bUsed     = TRUE;
1193                             bFirstTwo = bFirstTwo && m < 2;
1194                         }
1195                     }
1196                     if (!bUsed)
1197                     {
1198                         bKeep = TRUE;
1199                         if (debug)
1200                         {
1201                             fprintf(debug, " !used");
1202                         }
1203                     }
1204                 }
1205             }
1206
1207             if (!( bAllFD && bFirstTwo ) )
1208             {
1209                 /* Two atom bonded interactions include constraints.
1210                  * We need to remove constraints between vsite pairs that have
1211                  * a fixed distance due to being constructed from the same
1212                  * atoms, since this can be numerically unstable.
1213                  */
1214                 for (m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1215                 {
1216                     at1      = first_atoms[m];
1217                     at2      = first_atoms[(m+1) % vsnral];
1218                     bPresent = FALSE;
1219                     for (ftype = 0; ftype < F_NRE; ftype++)
1220                     {
1221                         if (interaction_function[ftype].flags & IF_CONSTRAINT)
1222                         {
1223                             for (j = 0; (j < plist[ftype].nr) && !bPresent; j++)
1224                             {
1225                                 /* all constraints until one matches */
1226                                 bPresent = ( ( (plist[ftype].param[j].ai() == at1) &&
1227                                                (plist[ftype].param[j].aj() == at2) ) ||
1228                                              ( (plist[ftype].param[j].ai() == at2) &&
1229                                                (plist[ftype].param[j].aj() == at1) ) );
1230                             }
1231                         }
1232                     }
1233                     if (!bPresent)
1234                     {
1235                         bKeep = TRUE;
1236                         if (debug)
1237                         {
1238                             fprintf(debug, " !bonded");
1239                         }
1240                     }
1241                 }
1242             }
1243         }
1244
1245         if (bKeep)
1246         {
1247             if (debug)
1248             {
1249                 fprintf(debug, " keeping");
1250             }
1251             /* now copy the bond to the new array */
1252             ps->param[kept_i] = ps->param[i];
1253             kept_i++;
1254         }
1255         else if (IS_CHEMBOND(cftype))
1256         {
1257             srenew(plist[F_CONNBONDS].param, plist[F_CONNBONDS].nr+1);
1258             plist[F_CONNBONDS].param[plist[F_CONNBONDS].nr] = ps->param[i];
1259             plist[F_CONNBONDS].nr++;
1260             nconverted++;
1261         }
1262         else
1263         {
1264             nremoved++;
1265         }
1266         if (debug)
1267         {
1268             fprintf(debug, "\n");
1269         }
1270     }
1271
1272     if (nremoved)
1273     {
1274         fprintf(stderr, "Removed   %4d %15ss with virtual sites, %5d left\n",
1275                 nremoved, interaction_function[cftype].longname, kept_i);
1276     }
1277     if (nconverted)
1278     {
1279         fprintf(stderr, "Converted %4d %15ss with virtual sites to connections, %5d left\n",
1280                 nconverted, interaction_function[cftype].longname, kept_i);
1281     }
1282     if (nOut)
1283     {
1284         fprintf(stderr, "Warning: removed %d %ss with vsite with %s construction\n"
1285                 "         This vsite construction does not guarantee constant "
1286                 "bond-length\n"
1287                 "         If the constructions were generated by pdb2gmx ignore "
1288                 "this warning\n",
1289                 nOut, interaction_function[cftype].longname,
1290                 interaction_function[F_VSITE3OUT].longname );
1291     }
1292     ps->nr = kept_i;
1293 }
1294
1295 static void clean_vsite_angles(t_params *plist, t_pindex pindex[],
1296                                int cftype, int vsite_type[],
1297                                at2vsitecon_t *at2vc)
1298 {
1299     int          i, j, k, m, n, nvsite, kept_i;
1300     int          atom, at1, at2;
1301     gmx_bool     bKeep, bUsed, bPresent, bAll3FAD, bFirstTwo;
1302     t_params    *ps;
1303
1304     ps     = &(plist[cftype]);
1305     kept_i = 0;
1306     for (i = 0; (i < ps->nr); i++) /* for all angles in the plist */
1307     {
1308         int            vsnral      = 0;
1309         const int     *first_atoms = nullptr;
1310
1311         bKeep    = FALSE;
1312         bAll3FAD = TRUE;
1313         /* check if all virtual sites are constructed from the same atoms */
1314         nvsite = 0;
1315         for (k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1316         {
1317             atom = ps->param[i].a[k];
1318             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1319             {
1320                 nvsite++;
1321                 bAll3FAD = bAll3FAD && (pindex[atom].ftype == F_VSITE3FAD);
1322                 if (nvsite == 1)
1323                 {
1324                     /* store construction atoms of first vsite */
1325                     vsnral      = NRAL(pindex[atom].ftype) - 1;
1326                     first_atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1;
1327                 }
1328                 else
1329                 {
1330                     GMX_ASSERT(vsnral != 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1331                     GMX_ASSERT(first_atoms != NULL, "If we've seen a vsite before, we know what its first atom index was");
1332                     /* check if this vsite is constructed from the same atoms */
1333                     if (vsnral == NRAL(pindex[atom].ftype)-1)
1334                     {
1335                         for (m = 0; (m < vsnral) && !bKeep; m++)
1336                         {
1337                             const int *atoms;
1338
1339                             bPresent = FALSE;
1340                             atoms    = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1;
1341                             for (n = 0; (n < vsnral) && !bPresent; n++)
1342                             {
1343                                 if (atoms[m] == first_atoms[n])
1344                                 {
1345                                     bPresent = TRUE;
1346                                 }
1347                             }
1348                             if (!bPresent)
1349                             {
1350                                 bKeep = TRUE;
1351                             }
1352                         }
1353                     }
1354                     else
1355                     {
1356                         bKeep = TRUE;
1357                     }
1358                 }
1359             }
1360         }
1361
1362         /* keep all angles with no virtual sites in them or
1363            with virtual sites with more than 3 constr. atoms */
1364         if (nvsite == 0 && vsnral > 3)
1365         {
1366             bKeep = TRUE;
1367         }
1368
1369         /* check if all non-vsite atoms are used in construction: */
1370         bFirstTwo = TRUE;
1371         for (k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1372         {
1373             atom = ps->param[i].a[k];
1374             if (vsite_type[atom] == NOTSET)
1375             {
1376                 bUsed = FALSE;
1377                 for (m = 0; (m < vsnral) && !bUsed; m++)
1378                 {
1379                     GMX_ASSERT(first_atoms != NULL, "If we've seen a vsite before, we know what its first atom index was");
1380
1381                     if (atom == first_atoms[m])
1382                     {
1383                         bUsed     = TRUE;
1384                         bFirstTwo = bFirstTwo && m < 2;
1385                     }
1386                 }
1387                 if (!bUsed)
1388                 {
1389                     bKeep = TRUE;
1390                 }
1391             }
1392         }
1393
1394         if (!( bAll3FAD && bFirstTwo ) )
1395         {
1396             /* check if all constructing atoms are constrained together */
1397             for (m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1398             {
1399                 at1      = first_atoms[m];
1400                 at2      = first_atoms[(m+1) % vsnral];
1401                 bPresent = FALSE;
1402                 for (j = 0; j < at2vc[at1].nr; j++)
1403                 {
1404                     if (at2vc[at1].aj[j] == at2)
1405                     {
1406                         bPresent = TRUE;
1407                     }
1408                 }
1409                 if (!bPresent)
1410                 {
1411                     bKeep = TRUE;
1412                 }
1413             }
1414         }
1415
1416         if (bKeep)
1417         {
1418             /* now copy the angle to the new array */
1419             ps->param[kept_i] = ps->param[i];
1420             kept_i++;
1421         }
1422     }
1423
1424     if (ps->nr != kept_i)
1425     {
1426         fprintf(stderr, "Removed   %4d %15ss with virtual sites, %5d left\n",
1427                 ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1428     }
1429     ps->nr = kept_i;
1430 }
1431
1432 static void clean_vsite_dihs(t_params *plist, t_pindex pindex[],
1433                              int cftype, int vsite_type[])
1434 {
1435     int       i, kept_i;
1436     t_params *ps;
1437
1438     ps = &(plist[cftype]);
1439
1440     kept_i = 0;
1441     for (i = 0; (i < ps->nr); i++) /* for all dihedrals in the plist */
1442     {
1443         int            k, m, n, nvsite;
1444         int            vsnral      = 0;
1445         const int     *first_atoms = nullptr;
1446         int            atom;
1447         gmx_bool       bKeep, bUsed, bPresent;
1448
1449
1450         bKeep = FALSE;
1451         /* check if all virtual sites are constructed from the same atoms */
1452         nvsite = 0;
1453         for (k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1454         {
1455             atom = ps->param[i].a[k];
1456             if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1457             {
1458                 if (nvsite == 0)
1459                 {
1460                     /* store construction atoms of first vsite */
1461                     vsnral      = NRAL(pindex[atom].ftype) - 1;
1462                     first_atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1;
1463                     if (debug)
1464                     {
1465                         fprintf(debug, "dih w. vsite: %d %d %d %d\n",
1466                                 ps->param[i].ai()+1, ps->param[i].aj()+1,
1467                                 ps->param[i].ak()+1, ps->param[i].al()+1);
1468                         fprintf(debug, "vsite %d from: %d %d %d\n",
1469                                 atom+1, first_atoms[0]+1, first_atoms[1]+1, first_atoms[2]+1);
1470                     }
1471                 }
1472                 else
1473                 {
1474                     GMX_ASSERT(vsnral != 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1475                     GMX_ASSERT(first_atoms != NULL, "If we've seen a vsite before, we know what its first atom index was");
1476                     /* check if this vsite is constructed from the same atoms */
1477                     if (vsnral == NRAL(pindex[atom].ftype)-1)
1478                     {
1479                         for (m = 0; (m < vsnral) && !bKeep; m++)
1480                         {
1481                             const int *atoms;
1482
1483                             bPresent = FALSE;
1484                             atoms    = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1;
1485                             for (n = 0; (n < vsnral) && !bPresent; n++)
1486                             {
1487                                 if (atoms[m] == first_atoms[n])
1488                                 {
1489                                     bPresent = TRUE;
1490                                 }
1491                             }
1492                             if (!bPresent)
1493                             {
1494                                 bKeep = TRUE;
1495                             }
1496                         }
1497                     }
1498                 }
1499                 /* TODO clean_site_bonds and _angles do this increment
1500                    at the top of the loop. Refactor this for
1501                    consistency */
1502                 nvsite++;
1503             }
1504         }
1505
1506         /* keep all dihedrals with no virtual sites in them */
1507         if (nvsite == 0)
1508         {
1509             bKeep = TRUE;
1510         }
1511
1512         /* check if all atoms in dihedral are either virtual sites, or used in
1513            construction of virtual sites. If so, keep it, if not throw away: */
1514         for (k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1515         {
1516             GMX_ASSERT(vsnral != 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1517             GMX_ASSERT(first_atoms != NULL, "If we've seen a vsite before, we know what its first atom index was");
1518             atom = ps->param[i].a[k];
1519             if (vsite_type[atom] == NOTSET)
1520             {
1521                 /* vsnral will be set here, we don't get here with nvsite==0 */
1522                 bUsed = FALSE;
1523                 for (m = 0; (m < vsnral) && !bUsed; m++)
1524                 {
1525                     if (atom == first_atoms[m])
1526                     {
1527                         bUsed = TRUE;
1528                     }
1529                 }
1530                 if (!bUsed)
1531                 {
1532                     bKeep = TRUE;
1533                     if (debug)
1534                     {
1535                         fprintf(debug, "unused atom in dih: %d\n", atom+1);
1536                     }
1537                 }
1538             }
1539         }
1540
1541         if (bKeep)
1542         {
1543             ps->param[kept_i] = ps->param[i];
1544             kept_i++;
1545         }
1546     }
1547
1548     if (ps->nr != kept_i)
1549     {
1550         fprintf(stderr, "Removed   %4d %15ss with virtual sites, %5d left\n",
1551                 ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1552     }
1553     ps->nr = kept_i;
1554 }
1555
1556 void clean_vsite_bondeds(t_params *plist, int natoms, gmx_bool bRmVSiteBds)
1557 {
1558     int            i, k, nvsite, ftype, vsite, parnr;
1559     int           *vsite_type;
1560     t_pindex      *pindex;
1561     at2vsitecon_t *at2vc;
1562
1563     pindex = nullptr; /* avoid warnings */
1564     /* make vsite_type array */
1565     snew(vsite_type, natoms);
1566     for (i = 0; i < natoms; i++)
1567     {
1568         vsite_type[i] = NOTSET;
1569     }
1570     nvsite = 0;
1571     for (ftype = 0; ftype < F_NRE; ftype++)
1572     {
1573         if (interaction_function[ftype].flags & IF_VSITE)
1574         {
1575             nvsite += plist[ftype].nr;
1576             i       = 0;
1577             while (i < plist[ftype].nr)
1578             {
1579                 vsite = plist[ftype].param[i].ai();
1580                 if (vsite_type[vsite] == NOTSET)
1581                 {
1582                     vsite_type[vsite] = ftype;
1583                 }
1584                 else
1585                 {
1586                     gmx_fatal(FARGS, "multiple vsite constructions for atom %d", vsite+1);
1587                 }
1588                 if (ftype == F_VSITEN)
1589                 {
1590                     while (i < plist[ftype].nr && plist[ftype].param[i].ai() == vsite)
1591                     {
1592                         i++;
1593                     }
1594                 }
1595                 else
1596                 {
1597                     i++;
1598                 }
1599             }
1600         }
1601     }
1602
1603     /* the rest only if we have virtual sites: */
1604     if (nvsite)
1605     {
1606         fprintf(stderr, "Cleaning up constraints %swith virtual sites\n",
1607                 bRmVSiteBds ? "and constant bonded interactions " : "");
1608
1609         /* Make a reverse list to avoid ninteractions^2 operations */
1610         at2vc = make_at2vsitecon(natoms, plist);
1611
1612         snew(pindex, natoms);
1613         for (ftype = 0; ftype < F_NRE; ftype++)
1614         {
1615             /* Here we skip VSITEN. In neary all practical use cases this
1616              * is not an issue, since VSITEN is intended for constructing
1617              * additional interaction sites, not for replacing normal atoms
1618              * with bonded interactions. Thus we do not expect constant
1619              * bonded interactions. If VSITEN does get used with constant
1620              * bonded interactions, these are not removed which only leads
1621              * to very minor extra computation and constant energy.
1622              * The only problematic case is onstraints between VSITEN
1623              * constructions with fixed distance (which is anyhow useless).
1624              * This will generate a fatal error in check_vsite_constraints.
1625              */
1626             if ((interaction_function[ftype].flags & IF_VSITE) &&
1627                 ftype != F_VSITEN)
1628             {
1629                 for (parnr = 0; (parnr < plist[ftype].nr); parnr++)
1630                 {
1631                     k               = plist[ftype].param[parnr].ai();
1632                     pindex[k].ftype = ftype;
1633                     pindex[k].parnr = parnr;
1634                 }
1635             }
1636         }
1637
1638         if (debug)
1639         {
1640             for (i = 0; i < natoms; i++)
1641             {
1642                 fprintf(debug, "atom %d vsite_type %s\n", i,
1643                         vsite_type[i] == NOTSET ? "NOTSET" :
1644                         interaction_function[vsite_type[i]].name);
1645             }
1646         }
1647
1648         /* remove interactions that include virtual sites */
1649         for (ftype = 0; ftype < F_NRE; ftype++)
1650         {
1651             if ( ( ( interaction_function[ftype].flags & IF_BOND ) && bRmVSiteBds ) ||
1652                  ( interaction_function[ftype].flags & IF_CONSTRAINT ) )
1653             {
1654                 if (interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT) )
1655                 {
1656                     clean_vsite_bonds (plist, pindex, ftype, vsite_type);
1657                 }
1658                 else if (interaction_function[ftype].flags & IF_ATYPE)
1659                 {
1660                     clean_vsite_angles(plist, pindex, ftype, vsite_type, at2vc);
1661                 }
1662                 else if ( (ftype == F_PDIHS) || (ftype == F_IDIHS) )
1663                 {
1664                     clean_vsite_dihs  (plist, pindex, ftype, vsite_type);
1665                 }
1666             }
1667         }
1668         /* check that no remaining constraints include virtual sites */
1669         for (ftype = 0; ftype < F_NRE; ftype++)
1670         {
1671             if (interaction_function[ftype].flags & IF_CONSTRAINT)
1672             {
1673                 check_vsite_constraints(plist, ftype, vsite_type);
1674             }
1675         }
1676
1677         done_at2vsitecon(natoms, at2vc);
1678     }
1679     sfree(pindex);
1680     sfree(vsite_type);
1681 }