src/kernel/vsite_parm.c

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