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