src/gromacs/gmxpreprocess/vsite_parm.c

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