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