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