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