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