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