src/kernel/vsite_parm.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <math.h>
  41 #include <string.h>
  42 #include "vsite_parm.h"
  43 #include "smalloc.h"
  44 #include "resall.h"
  45 #include "add_par.h"
  46 #include "vec.h"
  47 #include "toputil.h"
  48 #include "physics.h"
  49 #include "index.h"
  50 #include "names.h"
  51 #include "gmx_fatal.h"
  52 #include "string2.h"
  53 #include "physics.h"
  54
  55 typedef struct {
  56   t_iatom a[4];
  57   real    c;
  58 } t_mybonded;
  59
  60 typedef struct {
  61   int     ftype;
  62   t_param *param;
  63 } vsitebondparam_t;
  64
  65 typedef struct {
  66   int              nr;
  67   int              ftype;
  68   vsitebondparam_t *vsbp;
  69 } at2vsitebond_t;
  70
  71 typedef struct {
  72   int nr;
  73   int *aj;
  74 } at2vsitecon_t;
  75
  76 static int vsite_bond_nrcheck(int ftype)
  77 {
  78   int nrcheck;
  79
  80   if ((interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT | IF_ATYPE)) || (ftype == F_IDIHS))
  81     nrcheck = NRAL(ftype);
  82   else
  83     nrcheck = 0;
  84
  85   return nrcheck;
  86 }
  87
  88 static void enter_bonded(int nratoms, int *nrbonded, t_mybonded **bondeds,
  89                          t_param *param)
  90 {
  91   int j;
  92
  93   srenew(*bondeds, *nrbonded+1);
  94
  95   /* copy atom numbers */
  96   for(j=0; j<nratoms; j++)
  97     (*bondeds)[*nrbonded].a[j] = param->a[j];
  98   /* copy parameter */
  99   (*bondeds)[*nrbonded].c = param->C0;
 100
 101   (*nrbonded)++;
 102 }
 103
 104 static void get_bondeds(int nrat, t_iatom atoms[],
 105                         at2vsitebond_t *at2vb, t_params plist[],
 106                         int *nrbond, t_mybonded **bonds,
 107                         int *nrang,  t_mybonded **angles,
 108                         int *nridih, t_mybonded **idihs )
 109 {
 110   int     k,i,ftype,nrcheck;
 111   t_param *param;
 112
 113   for(k=0; k<nrat; k++) {
 114     for(i=0; i<at2vb[atoms[k]].nr; i++) {
 115       ftype = at2vb[atoms[k]].vsbp[i].ftype;
 116       param = at2vb[atoms[k]].vsbp[i].param;
 117       nrcheck = vsite_bond_nrcheck(ftype);
 118       /* abuse nrcheck to see if we're adding bond, angle or idih */
 119       switch (nrcheck) {
 120       case 2: enter_bonded(nrcheck,nrbond,bonds, param); break;
 121       case 3: enter_bonded(nrcheck,nrang, angles,param); break;
 122       case 4: enter_bonded(nrcheck,nridih,idihs, param); break;
 123       }
 124     }
 125   }
 126 }
 127
 128 static at2vsitebond_t *make_at2vsitebond(int natoms,t_params plist[])
 129 {
 130   bool *bVSI;
 131   int  ftype,i,j,nrcheck,nr;
 132   t_iatom *aa;
 133   at2vsitebond_t *at2vb;
 134
 135   snew(at2vb,natoms);
 136
 137   snew(bVSI,natoms);
 138   for(ftype=0; (ftype<F_NRE); ftype++) {
 139     if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
 140       for(i=0; (i<plist[ftype].nr); i++) {
 141         for(j=0; j<NRAL(ftype); j++)
 142           bVSI[plist[ftype].param[i].a[j]] = TRUE;
 143       }
 144     }
 145   }
 146
 147   for(ftype=0; (ftype<F_NRE); ftype++) {
 148     nrcheck = vsite_bond_nrcheck(ftype);
 149     if (nrcheck > 0) {
 150       for(i=0; (i<plist[ftype].nr); i++) {
 151         aa = plist[ftype].param[i].a;
 152         for(j=0; j<nrcheck; j++) {
 153           if (bVSI[aa[j]]) {
 154             nr = at2vb[aa[j]].nr;
 155             if (nr % 10 == 0)
 156               srenew(at2vb[aa[j]].vsbp,nr+10);
 157             at2vb[aa[j]].vsbp[nr].ftype = ftype;
 158             at2vb[aa[j]].vsbp[nr].param = &plist[ftype].param[i];
 159             at2vb[aa[j]].nr++;
 160           }
 161         }
 162       }
 163     }
 164   }
 165   sfree(bVSI);
 166
 167   return at2vb;
 168 }
 169
 170 static void done_at2vsitebond(int natoms,at2vsitebond_t *at2vb)
 171 {
 172   int i;
 173
 174   for(i=0; i<natoms; i++)
 175     if (at2vb[i].nr)
 176       sfree(at2vb[i].vsbp);
 177   sfree(at2vb);
 178 }
 179
 180 static at2vsitecon_t *make_at2vsitecon(int natoms,t_params plist[])
 181 {
 182   bool *bVSI;
 183   int  ftype,i,j,ai,aj,nr;
 184   at2vsitecon_t *at2vc;
 185
 186   snew(at2vc,natoms);
 187
 188   snew(bVSI,natoms);
 189   for(ftype=0; (ftype<F_NRE); ftype++) {
 190     if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
 191       for(i=0; (i<plist[ftype].nr); i++) {
 192         for(j=0; j<NRAL(ftype); j++)
 193           bVSI[plist[ftype].param[i].a[j]] = TRUE;
 194       }
 195     }
 196   }
 197
 198   for(ftype=0; (ftype<F_NRE); ftype++) {
 199     if (interaction_function[ftype].flags & IF_CONSTRAINT) {
 200       for(i=0; (i<plist[ftype].nr); i++) {
 201         ai = plist[ftype].param[i].AI;
 202         aj = plist[ftype].param[i].AJ;
 203         if (bVSI[ai] && bVSI[aj]) {
 204           /* Store forward direction */
 205           nr = at2vc[ai].nr;
 206           if (nr % 10 == 0)
 207             srenew(at2vc[ai].aj,nr+10);
 208           at2vc[ai].aj[nr] = aj;
 209           at2vc[ai].nr++;
 210           /* Store backward direction */
 211           nr = at2vc[aj].nr;
 212           if (nr % 10 == 0)
 213             srenew(at2vc[aj].aj,nr+10);
 214           at2vc[aj].aj[nr] = ai;
 215           at2vc[aj].nr++;
 216         }
 217       }
 218     }
 219   }
 220   sfree(bVSI);
 221
 222   return at2vc;
 223 }
 224
 225 static void done_at2vsitecon(int natoms,at2vsitecon_t *at2vc)
 226 {
 227   int i;
 228
 229   for(i=0; i<natoms; i++)
 230     if (at2vc[i].nr)
 231       sfree(at2vc[i].aj);
 232   sfree(at2vc);
 233 }
 234
 235 /* for debug */
 236 static void print_bad(FILE *fp,
 237                       int nrbond, t_mybonded *bonds,
 238                       int nrang,  t_mybonded *angles,
 239                       int nridih, t_mybonded *idihs )
 240 {
 241   int i;
 242
 243   if (nrbond) {
 244     fprintf(fp,"bonds:");
 245     for(i=0; i<nrbond; i++)
 246       fprintf(fp," %u-%u (%g)",
 247               bonds[i].AI+1, bonds[i].AJ+1, bonds[i].c);
 248     fprintf(fp,"\n");
 249   }
 250   if (nrang) {
 251     fprintf(fp,"angles:");
 252     for(i=0; i<nrang; i++)
 253       fprintf(fp," %u-%u-%u (%g)",
 254               angles[i].AI+1, angles[i].AJ+1,
 255               angles[i].AK+1, angles[i].c);
 256     fprintf(fp,"\n");
 257   }
 258   if (nridih) {
 259     fprintf(fp,"idihs:");
 260     for(i=0; i<nridih; i++)
 261       fprintf(fp," %u-%u-%u-%u (%g)",
 262               idihs[i].AI+1, idihs[i].AJ+1,
 263               idihs[i].AK+1, idihs[i].AL+1, idihs[i].c);
 264     fprintf(fp,"\n");
 265   }
 266 }
 267
 268 static void print_param(FILE *fp, int ftype, int i, t_param *param)
 269 {
 270   static int pass = 0;
 271   static int prev_ftype= NOTSET;
 272   static int prev_i    = NOTSET;
 273   int j;
 274
 275   if ( (ftype!=prev_ftype) || (i!=prev_i) ) {
 276     pass = 0;
 277     prev_ftype= ftype;
 278     prev_i    = i;
 279   }
 280   fprintf(fp,"(%d) plist[%s].param[%d]",
 281           pass,interaction_function[ftype].name,i);
 282   for(j=0; j<NRFP(ftype); j++)
 283     fprintf(fp,".c[%d]=%g ",j,param->c[j]);
 284   fprintf(fp,"\n");
 285   pass++;
 286 }
 287
 288 static real get_bond_length(int nrbond, t_mybonded bonds[],
 289                             t_iatom ai, t_iatom aj)
 290 {
 291   int  i;
 292   real bondlen;
 293
 294   bondlen=NOTSET;
 295   for (i=0; i < nrbond && (bondlen==NOTSET); i++) {
 296     /* check both ways */
 297     if ( ( (ai == bonds[i].AI) && (aj == bonds[i].AJ) ) ||
 298          ( (ai == bonds[i].AJ) && (aj == bonds[i].AI) ) )
 299       bondlen = bonds[i].c; /* note: bonds[i].c might be NOTSET */
 300   }
 301   return bondlen;
 302 }
 303
 304 static real get_angle(int nrang, t_mybonded angles[],
 305                       t_iatom ai, t_iatom aj, t_iatom ak)
 306 {
 307   int  i;
 308   real angle;
 309
 310   angle=NOTSET;
 311   for (i=0; i < nrang && (angle==NOTSET); i++) {
 312     /* check both ways */
 313     if ( ( (ai==angles[i].AI) && (aj==angles[i].AJ) && (ak==angles[i].AK) ) ||
 314          ( (ai==angles[i].AK) && (aj==angles[i].AJ) && (ak==angles[i].AI) ) )
 315       angle = DEG2RAD*angles[i].c;
 316   }
 317   return angle;
 318 }
 319
 320 static bool calc_vsite3_param(gpp_atomtype_t atype,
 321                               t_param *param, t_atoms *at,
 322                               int nrbond, t_mybonded *bonds,
 323                               int nrang,  t_mybonded *angles )
 324 {
 325   /* i = virtual site          |    ,k
 326    * j = 1st bonded heavy atom | i-j
 327    * k,l = 2nd bonded atoms    |    `l
 328    */
 329
 330   bool bXH3,bError;
 331   real bjk,bjl,a=-1,b=-1;
 332   /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
 333    * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 334   if (debug) {
 335     int i;
 336     for (i=0; i<4; i++)
 337       fprintf(debug,"atom %u type %s ",
 338               param->a[i]+1,
 339               get_atomtype_name(at->atom[param->a[i]].type,atype));
 340     fprintf(debug,"\n");
 341   }
 342   bXH3 =
 343     ( (gmx_strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MNH",3)==0) &&
 344       (gmx_strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MNH",3)==0) ) ||
 345     ( (gmx_strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MCH3",4)==0) &&
 346       (gmx_strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MCH3",4)==0) );
 347
 348   bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
 349   bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
 350   bError = (bjk==NOTSET) || (bjl==NOTSET);
 351   if (bXH3) {
 352     /* now we get some XH2/XH3 group specific construction */
 353     /* note: we call the heavy atom 'C' and the X atom 'N' */
 354     real bMM,bCM,bCN,bNH,aCNH,dH,rH,dM,rM;
 355     int aN;
 356
 357     /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 358     bError = bError || (bjk!=bjl);
 359
 360     /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
 361     aN = max(param->AK,param->AL)+1;
 362
 363     /* get common bonds */
 364     bMM = get_bond_length(nrbond, bonds, param->AK, param->AL);
 365     bCM = bjk;
 366     bCN = get_bond_length(nrbond, bonds, param->AJ, aN);
 367     bError = bError || (bMM==NOTSET) || (bCN==NOTSET);
 368
 369     /* calculate common things */
 370     rM  = 0.5*bMM;
 371     dM  = sqrt( sqr(bCM) - sqr(rM) );
 372
 373     /* are we dealing with the X atom? */
 374     if ( param->AI == aN ) {
 375       /* this is trivial */
 376       a = b = 0.5 * bCN/dM;
 377
 378     } else {
 379       /* get other bondlengths and angles: */
 380       bNH = get_bond_length(nrbond, bonds, aN, param->AI);
 381       aCNH= get_angle      (nrang, angles, param->AJ, aN, param->AI);
 382       bError = bError || (bNH==NOTSET) || (aCNH==NOTSET);
 383
 384       /* calculate */
 385       dH  = bCN - bNH * cos(aCNH);
 386       rH  = bNH * sin(aCNH);
 387
 388       a = 0.5 * ( dH/dM + rH/rM );
 389       b = 0.5 * ( dH/dM - rH/rM );
 390     }
 391   } else
 392     gmx_fatal(FARGS,"calc_vsite3_param not implemented for the general case "
 393                 "(atom %d)",param->AI+1);
 394
 395   param->C0 = a;
 396   param->C1 = b;
 397
 398   if (debug)
 399     fprintf(debug,"params for vsite3 %u: %g %g\n",
 400             param->AI+1,param->C0,param->C1);
 401
 402   return bError;
 403 }
 404
 405 static bool calc_vsite3fd_param(t_param *param,
 406                                 int nrbond, t_mybonded *bonds,
 407                                 int nrang,  t_mybonded *angles)
 408 {
 409   /* i = virtual site          |    ,k
 410    * j = 1st bonded heavy atom | i-j
 411    * k,l = 2nd bonded atoms    |    `l
 412    */
 413
 414   bool bError;
 415   real bij,bjk,bjl,aijk,aijl,rk,rl;
 416
 417   bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
 418   bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
 419   bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
 420   aijk= get_angle      (nrang, angles, param->AI, param->AJ, param->AK);
 421   aijl= get_angle      (nrang, angles, param->AI, param->AJ, param->AL);
 422   bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) ||
 423     (aijk==NOTSET) || (aijl==NOTSET);
 424
 425   rk = bjk * sin(aijk);
 426   rl = bjl * sin(aijl);
 427   param->C0 = rk / (rk + rl);
 428   param->C1 = -bij; /* 'bond'-length for fixed distance vsite */
 429
 430   if (debug)
 431     fprintf(debug,"params for vsite3fd %u: %g %g\n",
 432             param->AI+1,param->C0,param->C1);
 433   return bError;
 434 }
 435
 436 static bool calc_vsite3fad_param(t_param *param,
 437                                  int nrbond, t_mybonded *bonds,
 438                                  int nrang,  t_mybonded *angles)
 439 {
 440   /* i = virtual site          |
 441    * j = 1st bonded heavy atom | i-j
 442    * k = 2nd bonded heavy atom |    `k-l
 443    * l = 3d bonded heavy atom  |
 444    */
 445
 446   bool bSwapParity,bError;
 447   real bij,aijk;
 448
 449   bSwapParity = ( param->C1 == -1 );
 450
 451   bij  = get_bond_length(nrbond, bonds, param->AI, param->AJ);
 452   aijk = get_angle      (nrang, angles, param->AI, param->AJ, param->AK);
 453   bError = (bij==NOTSET) || (aijk==NOTSET);
 454
 455   param->C1 = bij;          /* 'bond'-length for fixed distance vsite */
 456   param->C0 = RAD2DEG*aijk; /* 'bond'-angle for fixed angle vsite */
 457
 458   if (bSwapParity)
 459     param->C0 = 360 - param->C0;
 460
 461   if (debug)
 462     fprintf(debug,"params for vsite3fad %u: %g %g\n",
 463             param->AI+1,param->C0,param->C1);
 464   return bError;
 465 }
 466
 467 static bool calc_vsite3out_param(gpp_atomtype_t atype,
 468                                  t_param *param, t_atoms *at,
 469                                  int nrbond, t_mybonded *bonds,
 470                                  int nrang,  t_mybonded *angles)
 471 {
 472   /* i = virtual site          |    ,k
 473    * j = 1st bonded heavy atom | i-j
 474    * k,l = 2nd bonded atoms    |    `l
 475    * NOTE: i is out of the j-k-l plane!
 476    */
 477
 478   bool bXH3,bError,bSwapParity;
 479   real bij,bjk,bjl,aijk,aijl,akjl,pijk,pijl,a,b,c;
 480
 481   /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
 482    * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
 483   if (debug) {
 484     int i;
 485     for (i=0; i<4; i++)
 486       fprintf(debug,"atom %u type %s ",
 487               param->a[i]+1,get_atomtype_name(at->atom[param->a[i]].type,atype));
 488     fprintf(debug,"\n");
 489   }
 490   bXH3 =
 491     ( (gmx_strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MNH",3)==0) &&
 492       (gmx_strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MNH",3)==0) ) ||
 493     ( (gmx_strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MCH3",4)==0) &&
 494       (gmx_strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MCH3",4)==0) );
 495
 496   /* check if construction parity must be swapped */
 497   bSwapParity = ( param->C1 == -1 );
 498
 499   bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
 500   bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
 501   bError = (bjk==NOTSET) || (bjl==NOTSET);
 502   if (bXH3) {
 503     /* now we get some XH3 group specific construction */
 504     /* note: we call the heavy atom 'C' and the X atom 'N' */
 505     real bMM,bCM,bCN,bNH,aCNH,dH,rH,rHx,rHy,dM,rM;
 506     int aN;
 507
 508     /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
 509     bError = bError || (bjk!=bjl);
 510
 511     /* the X atom (C or N) in the XH3 group is the first after the masses: */
 512     aN = max(param->AK,param->AL)+1;
 513
 514     /* get all bondlengths and angles: */
 515     bMM = get_bond_length(nrbond, bonds, param->AK, param->AL);
 516     bCM = bjk;
 517     bCN = get_bond_length(nrbond, bonds, param->AJ, aN);
 518     bNH = get_bond_length(nrbond, bonds, aN, param->AI);
 519     aCNH= get_angle      (nrang, angles, param->AJ, aN, param->AI);
 520     bError = bError ||
 521       (bMM==NOTSET) || (bCN==NOTSET) || (bNH==NOTSET) || (aCNH==NOTSET);
 522
 523     /* calculate */
 524     dH  = bCN - bNH * cos(aCNH);
 525     rH  = bNH * sin(aCNH);
 526     /* we assume the H's are symmetrically distributed */
 527     rHx = rH*cos(DEG2RAD*30);
 528     rHy = rH*sin(DEG2RAD*30);
 529     rM  = 0.5*bMM;
 530     dM  = sqrt( sqr(bCM) - sqr(rM) );
 531     a   = 0.5*( (dH/dM) - (rHy/rM) );
 532     b   = 0.5*( (dH/dM) + (rHy/rM) );
 533     c   = rHx / (2*dM*rM);
 534
 535   } else {
 536     /* this is the general construction */
 537
 538     bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
 539     aijk= get_angle      (nrang, angles, param->AI, param->AJ, param->AK);
 540     aijl= get_angle      (nrang, angles, param->AI, param->AJ, param->AL);
 541     akjl= get_angle      (nrang, angles, param->AK, param->AJ, param->AL);
 542     bError = bError ||
 543       (bij==NOTSET) || (aijk==NOTSET) || (aijl==NOTSET) || (akjl==NOTSET);
 544
 545     pijk = cos(aijk)*bij;
 546     pijl = cos(aijl)*bij;
 547     a = ( pijk + (pijk*cos(akjl)-pijl) * cos(akjl) / sqr(sin(akjl)) ) / bjk;
 548     b = ( pijl + (pijl*cos(akjl)-pijk) * cos(akjl) / sqr(sin(akjl)) ) / bjl;
 549     c = - sqrt( sqr(bij) -
 550                 ( sqr(pijk) - 2*pijk*pijl*cos(akjl) + sqr(pijl) )
 551                 / sqr(sin(akjl)) )
 552       / ( bjk*bjl*sin(akjl) );
 553   }
 554
 555   param->C0 = a;
 556   param->C1 = b;
 557   if (bSwapParity)
 558     param->C2 = -c;
 559   else
 560     param->C2 =  c;
 561   if (debug)
 562     fprintf(debug,"params for vsite3out %u: %g %g %g\n",
 563             param->AI+1,param->C0,param->C1,param->C2);
 564   return bError;
 565 }
 566
 567 static bool calc_vsite4fd_param(t_param *param,
 568                                 int nrbond, t_mybonded *bonds,
 569                                 int nrang,  t_mybonded *angles)
 570 {
 571   /* i = virtual site          |    ,k
 572    * j = 1st bonded heavy atom | i-j-m
 573    * k,l,m = 2nd bonded atoms  |    `l
 574    */
 575
 576   bool bError;
 577   real bij,bjk,bjl,bjm,aijk,aijl,aijm,akjm,akjl;
 578   real pk,pl,pm,cosakl,cosakm,sinakl,sinakm,cl,cm;
 579
 580   bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
 581   bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
 582   bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
 583   bjm = get_bond_length(nrbond, bonds, param->AJ, param->AM);
 584   aijk= get_angle      (nrang, angles, param->AI, param->AJ, param->AK);
 585   aijl= get_angle      (nrang, angles, param->AI, param->AJ, param->AL);
 586   aijm= get_angle      (nrang, angles, param->AI, param->AJ, param->AM);
 587   akjm= get_angle      (nrang, angles, param->AK, param->AJ, param->AM);
 588   akjl= get_angle      (nrang, angles, param->AK, param->AJ, param->AL);
 589   bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) || (bjm==NOTSET) ||
 590     (aijk==NOTSET) || (aijl==NOTSET) || (aijm==NOTSET) || (akjm==NOTSET) ||
 591     (akjl==NOTSET);
 592
 593   if (!bError) {
 594     pk = bjk*sin(aijk);
 595     pl = bjl*sin(aijl);
 596     pm = bjm*sin(aijm);
 597     cosakl = (cos(akjl) - cos(aijk)*cos(aijl)) / (sin(aijk)*sin(aijl));
 598     cosakm = (cos(akjm) - cos(aijk)*cos(aijm)) / (sin(aijk)*sin(aijm));
 599     if ( cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1 ) {
 600       fprintf(stderr,"virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
 601               param->AI+1,RAD2DEG*aijk,RAD2DEG*aijl,RAD2DEG*aijm);
 602       gmx_fatal(FARGS,"invalid construction in calc_vsite4fd for atom %d: "
 603                   "cosakl=%f, cosakm=%f\n",param->AI+1,cosakl,cosakm);
 604     }
 605     sinakl = sqrt(1-sqr(cosakl));
 606     sinakm = sqrt(1-sqr(cosakm));
 607
 608     /* note: there is a '+' because of the way the sines are calculated */
 609     cl = -pk / ( pl*cosakl - pk + pl*sinakl*(pm*cosakm-pk)/(pm*sinakm) );
 610     cm = -pk / ( pm*cosakm - pk + pm*sinakm*(pl*cosakl-pk)/(pl*sinakl) );
 611
 612     param->C0 = cl;
 613     param->C1 = cm;
 614     param->C2 = -bij;
 615     if (debug)
 616       fprintf(debug,"params for vsite4fd %u: %g %g %g\n",
 617               param->AI+1,param->C0,param->C1,param->C2);
 618   }
 619
 620   return bError;
 621 }
 622
 623
 624 static bool
 625 calc_vsite4fdn_param(t_param *param,
 626                      int nrbond, t_mybonded *bonds,
 627                      int nrang,  t_mybonded *angles)
 628 {
 629     /* i = virtual site          |    ,k
 630     * j = 1st bonded heavy atom | i-j-m
 631     * k,l,m = 2nd bonded atoms  |    `l
 632     */
 633
 634     bool bError;
 635     real bij,bjk,bjl,bjm,aijk,aijl,aijm;
 636     real pk,pl,pm,a,b;
 637
 638     bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
 639     bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
 640     bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
 641     bjm = get_bond_length(nrbond, bonds, param->AJ, param->AM);
 642     aijk= get_angle      (nrang, angles, param->AI, param->AJ, param->AK);
 643     aijl= get_angle      (nrang, angles, param->AI, param->AJ, param->AL);
 644     aijm= get_angle      (nrang, angles, param->AI, param->AJ, param->AM);
 645
 646     bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) || (bjm==NOTSET) ||
 647         (aijk==NOTSET) || (aijl==NOTSET) || (aijm==NOTSET);
 648
 649     if (!bError) {
 650
 651         /* Calculate component of bond j-k along the direction i-j */
 652         pk = -bjk*cos(aijk);
 653
 654         /* Calculate component of bond j-l along the direction i-j */
 655         pl = -bjl*cos(aijl);
 656
 657         /* Calculate component of bond j-m along the direction i-j */
 658         pm = -bjm*cos(aijm);
 659
 660         if(fabs(pl)<1000*GMX_REAL_MIN || fabs(pm)<1000*GMX_REAL_MIN)
 661         {
 662             fprintf(stderr,"virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
 663                     param->AI+1,RAD2DEG*aijk,RAD2DEG*aijl,RAD2DEG*aijm);
 664             gmx_fatal(FARGS,"invalid construction in calc_vsite4fdn for atom %d: "
 665                       "pl=%f, pm=%f\n",param->AI+1,pl,pm);
 666         }
 667
 668         a = pk/pl;
 669         b = pk/pm;
 670
 671         param->C0 = a;
 672         param->C1 = b;
 673         param->C2 = bij;
 674
 675         if (debug)
 676             fprintf(debug,"params for vsite4fdn %u: %g %g %g\n",
 677                     param->AI+1,param->C0,param->C1,param->C2);
 678     }
 679
 680     return bError;
 681 }
 682
 683
 684
 685 int set_vsites(bool bVerbose, t_atoms *atoms, gpp_atomtype_t atype,
 686                 t_params plist[])
 687 {
 688   int i,j,ftype;
 689   int nvsite,nrbond,nrang,nridih,nrset;
 690   bool bFirst,bSet,bERROR;
 691   at2vsitebond_t *at2vb;
 692   t_mybonded *bonds;
 693   t_mybonded *angles;
 694   t_mybonded *idihs;
 695
 696   bFirst = TRUE;
 697   bERROR = TRUE;
 698   nvsite=0;
 699   if (debug)
 700     fprintf(debug, "\nCalculating parameters for virtual sites\n");
 701
 702   /* Make a reverse list to avoid ninteractions^2 operations */
 703   at2vb = make_at2vsitebond(atoms->nr,plist);
 704
 705   for(ftype=0; (ftype<F_NRE); ftype++)
 706     if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
 707       nrset=0;
 708       nvsite+=plist[ftype].nr;
 709       for(i=0; (i<plist[ftype].nr); i++) {
 710         /* check if all parameters are set */
 711         bSet=TRUE;
 712         for(j=0; j<NRFP(ftype) && bSet; j++)
 713           bSet = plist[ftype].param[i].c[j]!=NOTSET;
 714
 715         if (debug) {
 716           fprintf(debug,"bSet=%s ",bool_names[bSet]);
 717           print_param(debug,ftype,i,&plist[ftype].param[i]);
 718         }
 719         if (!bSet) {
 720           if (bVerbose && bFirst) {
 721             fprintf(stderr,"Calculating parameters for virtual sites\n");
 722             bFirst=FALSE;
 723           }
 724
 725           nrbond=nrang=nridih=0;
 726           bonds = NULL;
 727           angles= NULL;
 728           idihs = NULL;
 729           nrset++;
 730           /* now set the vsite parameters: */
 731           get_bondeds(NRAL(ftype), plist[ftype].param[i].a, at2vb, plist,
 732                       &nrbond, &bonds, &nrang,  &angles, &nridih, &idihs);
 733           if (debug) {
 734             fprintf(debug, "Found %d bonds, %d angles and %d idihs "
 735                     "for virtual site %u (%s)\n",nrbond,nrang,nridih,
 736                     plist[ftype].param[i].AI+1,
 737                     interaction_function[ftype].longname);
 738             print_bad(debug, nrbond, bonds, nrang, angles, nridih, idihs);
 739           } /* debug */
 740           switch(ftype) {
 741           case F_VSITE3:
 742             bERROR =
 743               calc_vsite3_param(atype, &(plist[ftype].param[i]), atoms,
 744                                 nrbond, bonds, nrang, angles);
 745             break;
 746           case F_VSITE3FD:
 747             bERROR =
 748               calc_vsite3fd_param(&(plist[ftype].param[i]),
 749                                   nrbond, bonds, nrang, angles);
 750             break;
 751           case F_VSITE3FAD:
 752             bERROR =
 753               calc_vsite3fad_param(&(plist[ftype].param[i]),
 754                                    nrbond, bonds, nrang, angles);
 755             break;
 756           case F_VSITE3OUT:
 757             bERROR =
 758               calc_vsite3out_param(atype, &(plist[ftype].param[i]), atoms,
 759                                    nrbond, bonds, nrang, angles);
 760             break;
 761           case F_VSITE4FD:
 762             bERROR =
 763               calc_vsite4fd_param(&(plist[ftype].param[i]),
 764                                   nrbond, bonds, nrang, angles);
 765             break;
 766           case F_VSITE4FDN:
 767           bERROR =
 768               calc_vsite4fdn_param(&(plist[ftype].param[i]),
 769                                nrbond, bonds, nrang, angles);
 770           break;
 771           default:
 772             gmx_fatal(FARGS,"Automatic parameter generation not supported "
 773                         "for %s atom %d",
 774                         interaction_function[ftype].longname,
 775                         plist[ftype].param[i].AI+1);
 776           } /* switch */
 777           if (bERROR)
 778             gmx_fatal(FARGS,"Automatic parameter generation not supported "
 779                         "for %s atom %d for this bonding configuration",
 780                         interaction_function[ftype].longname,
 781                         plist[ftype].param[i].AI+1);
 782           sfree(bonds);
 783           sfree(angles);
 784           sfree(idihs);
 785         } /* if bSet */
 786       } /* for i */
 787       if (debug && plist[ftype].nr)
 788         fprintf(stderr,"Calculated parameters for %d out of %d %s atoms\n",
 789                 nrset,plist[ftype].nr,interaction_function[ftype].longname);
 790     } /* if IF_VSITE */
 791
 792   done_at2vsitebond(atoms->nr,at2vb);
 793
 794   return nvsite;
 795 }
 796
 797 void set_vsites_ptype(bool bVerbose, gmx_moltype_t *molt)
 798 {
 799   int ftype,i;
 800   int nra,nrd;
 801   t_ilist *il;
 802   t_iatom *ia,avsite;
 803
 804   if (bVerbose)
 805     fprintf(stderr,"Setting particle type to V for virtual sites\n");
 806   if (debug)
 807     fprintf(stderr,"checking %d functypes\n",F_NRE);
 808   for(ftype=0; ftype<F_NRE; ftype++) {
 809     il = &molt->ilist[ftype];
 810     if (interaction_function[ftype].flags & IF_VSITE) {
 811       nra    = interaction_function[ftype].nratoms;
 812       nrd    = il->nr;
 813       ia     = il->iatoms;
 814
 815       if (debug && nrd)
 816         fprintf(stderr,"doing %d %s virtual sites\n",
 817                 (nrd / (nra+1)),interaction_function[ftype].longname);
 818
 819       for(i=0; (i<nrd); ) {
 820         /* The virtual site */
 821         avsite = ia[1];
 822         molt->atoms.atom[avsite].ptype = eptVSite;
 823
 824         i  += nra+1;
 825         ia += nra+1;
 826       }
 827     }
 828   }
 829
 830 }
 831
 832 typedef struct {
 833   int ftype,parnr;
 834 } t_pindex;
 835
 836 static void check_vsite_constraints(t_params *plist,
 837                                     int cftype, int vsite_type[])
 838 {
 839   int      i,k,n;
 840   atom_id  atom;
 841   t_params *ps;
 842
 843   n=0;
 844   ps = &(plist[cftype]);
 845   for(i=0; (i<ps->nr); i++)
 846     for(k=0; k<2; k++) {
 847       atom = ps->param[i].a[k];
 848       if (vsite_type[atom]!=NOTSET) {
 849         fprintf(stderr,"ERROR: Cannot have constraint (%u-%u) with virtual site (%u)\n",
 850                 ps->param[i].AI+1, ps->param[i].AJ+1, atom+1);
 851         n++;
 852       }
 853     }
 854   if (n)
 855     gmx_fatal(FARGS,"There were %d virtual sites involved in constraints",n);
 856 }
 857
 858 static void clean_vsite_bonds(t_params *plist, t_pindex pindex[],
 859                             int cftype, int vsite_type[])
 860 {
 861   int      ftype,i,j,parnr,k,l,m,n,nvsite,nOut,kept_i,vsnral,vsitetype;
 862   int      nconverted,nremoved;
 863   atom_id  atom,oatom,constr,at1,at2;
 864   atom_id  vsiteatoms[MAXATOMLIST];
 865   bool     bKeep,bRemove,bUsed,bPresent,bThisFD,bThisOUT,bAllFD,bFirstTwo;
 866   t_params *ps;
 867
 868   if (cftype == F_CONNBONDS)
 869     return;
 870
 871   ps = &(plist[cftype]);
 872   vsnral=0;
 873   kept_i=0;
 874   nconverted=0;
 875   nremoved=0;
 876   nOut=0;
 877   for(i=0; (i<ps->nr); i++) { /* for all bonds in the plist */
 878     bKeep=FALSE;
 879     bRemove=FALSE;
 880     bAllFD=TRUE;
 881     /* check if all virtual sites are constructed from the same atoms */
 882     nvsite=0;
 883     if(debug)
 884       fprintf(debug,"constr %u %u:",ps->param[i].AI+1,ps->param[i].AJ+1);
 885     for(k=0; (k<2) && !bKeep && !bRemove; k++) {
 886       /* for all atoms in the bond */
 887       atom = ps->param[i].a[k];
 888       if (vsite_type[atom]!=NOTSET) {
 889         if(debug) {
 890           fprintf(debug," d%d[%d: %d %d %d]",k,atom+1,
 891                   plist[pindex[atom].ftype].param[pindex[atom].parnr].AJ+1,
 892                   plist[pindex[atom].ftype].param[pindex[atom].parnr].AK+1,
 893                   plist[pindex[atom].ftype].param[pindex[atom].parnr].AL+1);
 894         }
 895         nvsite++;
 896         bThisFD = ( (pindex[atom].ftype == F_VSITE3FD ) ||
 897                 (pindex[atom].ftype == F_VSITE3FAD) ||
 898                 (pindex[atom].ftype == F_VSITE4FD ) ||
 899                 (pindex[atom].ftype == F_VSITE4FDN ) );
 900         bThisOUT= ( (pindex[atom].ftype == F_VSITE3OUT) &&
 901                     (interaction_function[cftype].flags & IF_CONSTRAINT) );
 902         bAllFD = bAllFD && bThisFD;
 903         if (bThisFD || bThisOUT) {
 904           if(debug)fprintf(debug," %s",bThisOUT?"out":"fd");
 905           oatom = ps->param[i].a[1-k]; /* the other atom */
 906           if ( vsite_type[oatom]==NOTSET &&
 907                oatom==plist[pindex[atom].ftype].param[pindex[atom].parnr].AJ ){
 908             /* if the other atom isn't a vsite, and it is AI */
 909             bRemove=TRUE;
 910             if (bThisOUT)
 911               nOut++;
 912             if(debug)fprintf(debug," D-AI");
 913           }
 914         }
 915         if (!bRemove) {
 916           if (nvsite==1) {
 917             /* if this is the first vsite we encounter then
 918                store construction atoms */
 919             vsnral=NRAL(pindex[atom].ftype)-1;
 920             for(m=0; (m<vsnral); m++)
 921               vsiteatoms[m]=
 922                 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
 923           } else {
 924             /* if it is not the first then
 925                check if this vsite is constructed from the same atoms */
 926             if (vsnral == NRAL(pindex[atom].ftype)-1 )
 927               for(m=0; (m<vsnral) && !bKeep; m++) {
 928                 bPresent=FALSE;
 929                 constr=
 930                   plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
 931                 for(n=0; (n<vsnral) && !bPresent; n++)
 932                   if (constr == vsiteatoms[n])
 933                     bPresent=TRUE;
 934                 if (!bPresent) {
 935                   bKeep=TRUE;
 936                   if(debug)fprintf(debug," !present");
 937                 }
 938               }
 939             else {
 940               bKeep=TRUE;
 941               if(debug)fprintf(debug," !same#at");
 942             }
 943           }
 944         }
 945       }
 946     }
 947
 948     if (bRemove)
 949       bKeep=FALSE;
 950     else {
 951       /* if we have no virtual sites in this bond, keep it */
 952       if (nvsite==0) {
 953         if (debug)fprintf(debug," no vsite");
 954         bKeep=TRUE;
 955       }
 956
 957       /* check if all non-vsite atoms are used in construction: */
 958       bFirstTwo=TRUE;
 959       for(k=0; (k<2) && !bKeep; k++) { /* for all atoms in the bond */
 960         atom = ps->param[i].a[k];
 961         if (vsite_type[atom]==NOTSET) {
 962           bUsed=FALSE;
 963           for(m=0; (m<vsnral) && !bUsed; m++)
 964             if (atom == vsiteatoms[m]) {
 965               bUsed=TRUE;
 966               bFirstTwo = bFirstTwo && m<2;
 967             }
 968           if (!bUsed) {
 969             bKeep=TRUE;
 970             if(debug)fprintf(debug," !used");
 971           }
 972         }
 973       }
 974
 975       if ( ! ( bAllFD && bFirstTwo ) )
 976         /* check if all constructing atoms are constrained together */
 977         for (m=0; m<vsnral && !bKeep; m++) { /* all constr. atoms */
 978           at1 = vsiteatoms[m];
 979           at2 = vsiteatoms[(m+1) % vsnral];
 980           bPresent=FALSE;
 981           for (ftype=0; ftype<F_NRE; ftype++)
 982             if ( interaction_function[ftype].flags & IF_CONSTRAINT )
 983               for (j=0; (j<plist[ftype].nr) && !bPresent; j++)
 984                 /* all constraints until one matches */
 985                 bPresent = ( ( (plist[ftype].param[j].AI == at1) &&
 986                                (plist[ftype].param[j].AJ == at2) ) ||
 987                              ( (plist[ftype].param[j].AI == at2) &&
 988                                (plist[ftype].param[j].AJ == at1) ) );
 989           if (!bPresent) {
 990             bKeep=TRUE;
 991             if(debug)fprintf(debug," !bonded");
 992           }
 993         }
 994     }
 995
 996     if ( bKeep ) {
 997       if(debug)fprintf(debug," keeping");
 998       /* now copy the bond to the new array */
 999       memcpy(&(ps->param[kept_i]),
1000              &(ps->param[i]),(size_t)sizeof(ps->param[0]));
1001       kept_i++;
1002     } else if (IS_CHEMBOND(cftype)) {
1003       srenew(plist[F_CONNBONDS].param,plist[F_CONNBONDS].nr+1);
1004       memcpy(&(plist[F_CONNBONDS].param[plist[F_CONNBONDS].nr]),
1005              &(ps->param[i]),(size_t)sizeof(plist[F_CONNBONDS].param[0]));
1006       plist[F_CONNBONDS].nr++;
1007       nconverted++;
1008     } else
1009       nremoved++;
1010     if(debug)fprintf(debug,"\n");
1011   }
1012
1013   if (nremoved)
1014     fprintf(stderr,"Removed   %4d %15ss with virtual sites, %5d left\n",
1015             nremoved, interaction_function[cftype].longname, kept_i);
1016   if (nconverted)
1017     fprintf(stderr,"Converted %4d %15ss with virtual sites to connections, %5d left\n",
1018             nconverted, interaction_function[cftype].longname, kept_i);
1019   if (nOut)
1020     fprintf(stderr,"Warning: removed %d %ss with vsite with %s construction\n"
1021             "         This vsite construction does not guarantee constant "
1022             "bond-length\n"
1023             "         If the constructions were generated by pdb2gmx ignore "
1024             "this warning\n",
1025             nOut, interaction_function[cftype].longname,
1026             interaction_function[F_VSITE3OUT].longname );
1027   ps->nr=kept_i;
1028 }
1029
1030 static void clean_vsite_angles(t_params *plist, t_pindex pindex[],
1031                                int cftype, int vsite_type[],
1032                                at2vsitecon_t *at2vc)
1033 {
1034   int      i,j,parnr,k,l,m,n,nvsite,kept_i,vsnral,vsitetype;
1035   atom_id  atom,constr,at1,at2;
1036   atom_id  vsiteatoms[MAXATOMLIST];
1037   bool     bKeep,bUsed,bPresent,bAll3FAD,bFirstTwo;
1038   t_params *ps;
1039
1040   ps = &(plist[cftype]);
1041   vsnral=0;
1042   kept_i=0;
1043   for(i=0; (i<ps->nr); i++) { /* for all angles in the plist */
1044     bKeep=FALSE;
1045     bAll3FAD=TRUE;
1046     /* check if all virtual sites are constructed from the same atoms */
1047     nvsite=0;
1048     for(k=0; (k<3) && !bKeep; k++) { /* for all atoms in the angle */
1049       atom = ps->param[i].a[k];
1050       if (vsite_type[atom]!=NOTSET) {
1051         nvsite++;
1052         bAll3FAD = bAll3FAD && (pindex[atom].ftype == F_VSITE3FAD);
1053         if (nvsite==1) {
1054           /* store construction atoms of first vsite */
1055           vsnral=NRAL(pindex[atom].ftype)-1;
1056           for(m=0; (m<vsnral); m++)
1057             vsiteatoms[m]=
1058               plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1059         } else
1060           /* check if this vsite is constructed from the same atoms */
1061           if (vsnral == NRAL(pindex[atom].ftype)-1 )
1062             for(m=0; (m<vsnral) && !bKeep; m++) {
1063               bPresent=FALSE;
1064               constr=
1065                 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1066               for(n=0; (n<vsnral) && !bPresent; n++)
1067                 if (constr == vsiteatoms[n])
1068                   bPresent=TRUE;
1069               if (!bPresent)
1070                 bKeep=TRUE;
1071             }
1072           else
1073             bKeep=TRUE;
1074       }
1075     }
1076
1077     /* keep all angles with no virtual sites in them or
1078        with virtual sites with more than 3 constr. atoms */
1079     if ( nvsite == 0 && vsnral > 3 )
1080       bKeep=TRUE;
1081
1082     /* check if all non-vsite atoms are used in construction: */
1083     bFirstTwo=TRUE;
1084     for(k=0; (k<3) && !bKeep; k++) { /* for all atoms in the angle */
1085       atom = ps->param[i].a[k];
1086       if (vsite_type[atom]==NOTSET) {
1087         bUsed=FALSE;
1088         for(m=0; (m<vsnral) && !bUsed; m++)
1089           if (atom == vsiteatoms[m]) {
1090             bUsed=TRUE;
1091             bFirstTwo = bFirstTwo && m<2;
1092           }
1093         if (!bUsed)
1094           bKeep=TRUE;
1095       }
1096     }
1097
1098     if ( ! ( bAll3FAD && bFirstTwo ) )
1099       /* check if all constructing atoms are constrained together */
1100       for (m=0; m<vsnral && !bKeep; m++) { /* all constr. atoms */
1101         at1 = vsiteatoms[m];
1102         at2 = vsiteatoms[(m+1) % vsnral];
1103         bPresent=FALSE;
1104         for(j=0; j<at2vc[at1].nr; j++) {
1105           if (at2vc[at1].aj[j] == at2)
1106             bPresent = TRUE;
1107         }
1108         if (!bPresent)
1109           bKeep=TRUE;
1110       }
1111
1112     if ( bKeep ) {
1113       /* now copy the angle to the new array */
1114       memcpy(&(ps->param[kept_i]),
1115              &(ps->param[i]),(size_t)sizeof(ps->param[0]));
1116       kept_i++;
1117     }
1118   }
1119
1120   if (ps->nr != kept_i)
1121     fprintf(stderr,"Removed   %4d %15ss with virtual sites, %5d left\n",
1122             ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1123   ps->nr=kept_i;
1124 }
1125
1126 static void clean_vsite_dihs(t_params *plist, t_pindex pindex[],
1127                            int cftype, int vsite_type[])
1128 {
1129   int      ftype,i,parnr,k,l,m,n,nvsite,kept_i,vsnral;
1130   atom_id  atom,constr;
1131   atom_id  vsiteatoms[3];
1132   bool     bKeep,bUsed,bPresent;
1133   t_params *ps;
1134
1135   ps = &(plist[cftype]);
1136
1137   vsnral=0;
1138   kept_i=0;
1139   for(i=0; (i<ps->nr); i++) { /* for all dihedrals in the plist */
1140     bKeep=FALSE;
1141     /* check if all virtual sites are constructed from the same atoms */
1142     nvsite=0;
1143     for(k=0; (k<4) && !bKeep; k++) { /* for all atoms in the dihedral */
1144       atom = ps->param[i].a[k];
1145       if (vsite_type[atom]!=NOTSET) {
1146         nvsite++;
1147         if (nvsite==1) {
1148           /* store construction atoms of first vsite */
1149           vsnral=NRAL(pindex[atom].ftype)-1;
1150           for(m=0; (m<vsnral); m++)
1151             vsiteatoms[m]=
1152               plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1153           if (debug) {
1154             fprintf(debug,"dih w. vsite: %u %u %u %u\n",
1155                     ps->param[i].AI+1,ps->param[i].AJ+1,
1156                     ps->param[i].AK+1,ps->param[i].AL+1);
1157             fprintf(debug,"vsite %u from: %u %u %u\n",
1158                     atom+1,vsiteatoms[0]+1,vsiteatoms[1]+1,vsiteatoms[2]+1);
1159           }
1160         } else
1161           /* check if this vsite is constructed from the same atoms */
1162           if (vsnral == NRAL(pindex[atom].ftype)-1 )
1163             for(m=0; (m<vsnral) && !bKeep; m++) {
1164               bPresent=FALSE;
1165               constr=
1166                 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1167               for(n=0; (n<vsnral) && !bPresent; n++)
1168                 if (constr == vsiteatoms[n])
1169                   bPresent=TRUE;
1170               if (!bPresent)
1171                 bKeep=TRUE;
1172             }
1173       }
1174     }
1175
1176     /* keep all dihedrals with no virtual sites in them */
1177     if (nvsite==0)
1178       bKeep=TRUE;
1179
1180     /* check if all atoms in dihedral are either virtual sites, or used in
1181        construction of virtual sites. If so, keep it, if not throw away: */
1182     for(k=0; (k<4) && !bKeep; k++) { /* for all atoms in the dihedral */
1183       atom = ps->param[i].a[k];
1184       if (vsite_type[atom]==NOTSET) {
1185         bUsed=FALSE;
1186         for(m=0; (m<vsnral) && !bUsed; m++)
1187           if (atom == vsiteatoms[m])
1188             bUsed=TRUE;
1189         if (!bUsed) {
1190           bKeep=TRUE;
1191           if (debug) fprintf(debug,"unused atom in dih: %u\n",atom+1);
1192         }
1193       }
1194     }
1195
1196     if ( bKeep ) {
1197       memcpy(&(ps->param[kept_i]),
1198              &(ps->param[i]),(size_t)sizeof(ps->param[0]));
1199       kept_i++;
1200     }
1201   }
1202
1203   if (ps->nr != kept_i)
1204     fprintf(stderr,"Removed   %4d %15ss with virtual sites, %5d left\n",
1205             ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1206   ps->nr=kept_i;
1207 }
1208
1209 void clean_vsite_bondeds(t_params *plist, int natoms, bool bRmVSiteBds)
1210 {
1211   int i,k,nvsite,ftype,vsite,parnr;
1212   int *vsite_type;
1213   t_pindex *pindex;
1214   at2vsitecon_t *at2vc;
1215
1216   pindex=0; /* avoid warnings */
1217   /* make vsite_type array */
1218   snew(vsite_type,natoms);
1219   for(i=0; i<natoms; i++)
1220     vsite_type[i]=NOTSET;
1221   nvsite=0;
1222   for(ftype=0; ftype<F_NRE; ftype++)
1223     if (interaction_function[ftype].flags & IF_VSITE) {
1224       nvsite+=plist[ftype].nr;
1225       i = 0;
1226       while (i < plist[ftype].nr) {
1227         vsite = plist[ftype].param[i].AI;
1228         if ( vsite_type[vsite] == NOTSET)
1229           vsite_type[vsite] = ftype;
1230         else
1231           gmx_fatal(FARGS,"multiple vsite constructions for atom %d",vsite+1);
1232         if (ftype == F_VSITEN) {
1233           while (i < plist[ftype].nr && plist[ftype].param[i].AI == vsite)
1234             i++;
1235         } else {
1236           i++;
1237         }
1238       }
1239     }
1240
1241   /* the rest only if we have virtual sites: */
1242   if (nvsite) {
1243     fprintf(stderr,"Cleaning up constraints %swith virtual sites\n",
1244             bRmVSiteBds?"and constant bonded interactions ":"");
1245
1246     /* Make a reverse list to avoid ninteractions^2 operations */
1247     at2vc = make_at2vsitecon(natoms,plist);
1248
1249     snew(pindex,natoms);
1250     for(ftype=0; ftype<F_NRE; ftype++) {
1251       if ((interaction_function[ftype].flags & IF_VSITE) &&
1252           ftype != F_VSITEN) {
1253         for (parnr=0; (parnr<plist[ftype].nr); parnr++) {
1254           k=plist[ftype].param[parnr].AI;
1255           pindex[k].ftype=ftype;
1256           pindex[k].parnr=parnr;
1257         }
1258       }
1259     }
1260
1261     if (debug)
1262       for(i=0; i<natoms; i++)
1263         fprintf(debug,"atom %d vsite_type %s\n",i,
1264                 vsite_type[i]==NOTSET ? "NOTSET" :
1265                 interaction_function[vsite_type[i]].name);
1266
1267     /* remove things with vsite atoms */
1268     for(ftype=0; ftype<F_NRE; ftype++)
1269       if ( ( ( interaction_function[ftype].flags & IF_BOND ) && bRmVSiteBds ) ||
1270            ( interaction_function[ftype].flags & IF_CONSTRAINT ) ) {
1271         if (interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT) )
1272           clean_vsite_bonds (plist, pindex, ftype, vsite_type);
1273         else if (interaction_function[ftype].flags & IF_ATYPE)
1274           clean_vsite_angles(plist, pindex, ftype, vsite_type, at2vc);
1275         else if ( (ftype==F_PDIHS) || (ftype==F_IDIHS) )
1276           clean_vsite_dihs  (plist, pindex, ftype, vsite_type);
1277       }
1278     /* check if we have constraints left with virtual sites in them */
1279     for(ftype=0; ftype<F_NRE; ftype++)
1280       if (interaction_function[ftype].flags & IF_CONSTRAINT)
1281         check_vsite_constraints(plist, ftype, vsite_type);
1282
1283     done_at2vsitecon(natoms,at2vc);
1284   }
1285   sfree(pindex);
1286   sfree(vsite_type);
1287 }