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