src/kernel/toppush.c

   1 /*
   2  * toppush.c
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <ctype.h>
  41 #include <math.h>
  42
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "string2.h"
  47 #include "names.h"
  48 #include "toputil.h"
  49 #include "toppush.h"
  50 #include "topdirs.h"
  51 #include "readir.h"
  52 #include "symtab.h"
  53 #include "gmx_fatal.h"
  54 #include "gpp_atomtype.h"
  55 #include "gpp_bond_atomtype.h"
  56
  57 void generate_nbparams(int comb,int ftype,t_params *plist,t_atomtype atype)
  58 {
  59   int   i,j,k=-1,nf;
  60   int   nr,nrfp;
  61   real  c,bi,bj,ci,cj,ci0,ci1,ci2,cj0,cj1,cj2;
  62   char  errbuf[256];
  63
  64   /* Lean mean shortcuts */
  65   nr   = get_atomtype_ntypes(atype);
  66   nrfp = NRFP(ftype);
  67   snew(plist->param,nr*nr);
  68   plist->nr=nr*nr;
  69
  70   /* Fill the matrix with force parameters */
  71   switch(ftype) {
  72   case F_LJ:
  73     switch (comb) {
  74     case eCOMB_GEOMETRIC:
  75       /* Gromos rules */
  76       for(i=k=0; (i<nr); i++) {
  77         for(j=0; (j<nr); j++,k++) {
  78           for(nf=0; (nf<nrfp); nf++) {
  79             ci = get_atomtype_nbparam(i,nf,atype);
  80             cj = get_atomtype_nbparam(j,nf,atype);
  81             c = sqrt(ci * cj);
  82             plist->param[k].c[nf]      = c;
  83           }
  84         }
  85       }
  86       break;
  87
  88     case eCOMB_ARITHMETIC:
  89       /* c0 and c1 are epsilon and sigma */
  90       for (i=k=0; (i < nr); i++)
  91         for (j=0; (j < nr); j++,k++) {
  92           ci0 = get_atomtype_nbparam(i,0,atype);
  93           cj0 = get_atomtype_nbparam(j,0,atype);
  94           ci1 = get_atomtype_nbparam(i,1,atype);
  95           cj1 = get_atomtype_nbparam(j,1,atype);
  96           plist->param[k].c[0] = (ci0+cj0)*0.5;
  97           plist->param[k].c[1] = sqrt(ci1*cj1);
  98         }
  99
 100       break;
 101     case eCOMB_GEOM_SIG_EPS:
 102       /* c0 and c1 are epsilon and sigma */
 103       for (i=k=0; (i < nr); i++)
 104         for (j=0; (j < nr); j++,k++) {
 105           ci0 = get_atomtype_nbparam(i,0,atype);
 106           cj0 = get_atomtype_nbparam(j,0,atype);
 107           ci1 = get_atomtype_nbparam(i,1,atype);
 108           cj1 = get_atomtype_nbparam(j,1,atype);
 109           plist->param[k].c[0] = sqrt(ci0*cj0);
 110           plist->param[k].c[1] = sqrt(ci1*cj1);
 111         }
 112
 113       break;
 114     default:
 115       gmx_fatal(FARGS,"No such combination rule %d",comb);
 116     }
 117     if (plist->nr != k)
 118       gmx_incons("Topology processing, generate nb parameters");
 119     break;
 120
 121   case F_BHAM:
 122     /* Buckingham rules */
 123     for(i=k=0; (i<nr); i++)
 124       for(j=0; (j<nr); j++,k++) {
 125         ci0 = get_atomtype_nbparam(i,0,atype);
 126         cj0 = get_atomtype_nbparam(j,0,atype);
 127         ci2 = get_atomtype_nbparam(i,2,atype);
 128         cj2 = get_atomtype_nbparam(j,2,atype);
 129         bi = get_atomtype_nbparam(i,1,atype);
 130         bj = get_atomtype_nbparam(j,1,atype);
 131         plist->param[k].c[0] = sqrt(ci0 * cj0);
 132         if ((bi == 0) || (bj == 0))
 133           plist->param[k].c[1] = 0;
 134         else
 135           plist->param[k].c[1] = 2.0/(1/bi+1/bj);
 136         plist->param[k].c[2] = sqrt(ci2 * cj2);
 137       }
 138
 139     break;
 140   default:
 141     sprintf(errbuf,"Invalid nonbonded type %s",
 142             interaction_function[ftype].longname);
 143     warning_error(errbuf);
 144   }
 145 }
 146
 147 static void realloc_nb_params(t_atomtype at,
 148                               t_nbparam ***nbparam,t_nbparam ***pair)
 149 {
 150   /* Add space in the non-bonded parameters matrix */
 151   int atnr = get_atomtype_ntypes(at);
 152   srenew(*nbparam,atnr);
 153   snew((*nbparam)[atnr-1],atnr);
 154   if (pair) {
 155     srenew(*pair,atnr);
 156     snew((*pair)[atnr-1],atnr);
 157   }
 158 }
 159
 160 static void copy_B_from_A(int ftype,double *c)
 161 {
 162   int nrfpA,nrfpB,i;
 163
 164   nrfpA = NRFPA(ftype);
 165   nrfpB = NRFPB(ftype);
 166
 167   /* Copy the B parameters from the first nrfpB A parameters */
 168   for(i=0; (i<nrfpB); i++) {
 169     c[nrfpA+i] = c[i];
 170   }
 171 }
 172
 173 void push_at (t_symtab *symtab, t_atomtype at, t_bond_atomtype bat,
 174               char *line,int nb_funct,
 175               t_nbparam ***nbparam,t_nbparam ***pair)
 176 {
 177   typedef struct {
 178     char *entry;
 179     int  ptype;
 180   } t_xlate;
 181   t_xlate xl[eptNR] = {
 182     { "A",   eptAtom },
 183     { "N",   eptNucleus },
 184     { "S",   eptShell },
 185     { "B",   eptBond },
 186     { "V",   eptVSite },
 187   };
 188
 189   int    nr,i,nfields,j,pt,nfp0=-1;
 190   int    batype_nr,nread;
 191   char   type[STRLEN],btype[STRLEN],ptype[STRLEN];
 192   double m,q;
 193   double c[MAXFORCEPARAM];
 194   double radius,vol,surftens;
 195   char   tmpfield[12][100]; /* Max 12 fields of width 100 */
 196   char   errbuf[256];
 197   t_atom  *atom;
 198   t_param *param;
 199   int    atomnr;
 200   bool   have_atomic_number;
 201   bool   have_bonded_type;
 202
 203   snew(atom,1);
 204   snew(param,1);
 205
 206   /* First assign input line to temporary array */
 207   nfields=sscanf(line,"%s%s%s%s%s%s%s%s%s%s%s%s",
 208                  tmpfield[0],tmpfield[1],tmpfield[2],tmpfield[3],tmpfield[4],tmpfield[5],
 209                  tmpfield[6],tmpfield[7],tmpfield[8],tmpfield[9],tmpfield[10],tmpfield[11]);
 210
 211   /* Comments on optional fields in the atomtypes section:
 212    *
 213    * The force field format is getting a bit old. For OPLS-AA we needed
 214    * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
 215    * we also needed the atomic numbers.
 216    * To avoid making all old or user-generated force fields unusable we
 217    * have introduced both these quantities as optional columns, and do some
 218    * acrobatics to check whether they are present or not.
 219    * This will all look much nicer when we switch to XML... sigh.
 220    *
 221    * Field 0 (mandatory) is the nonbonded type name. (string)
 222    * Field 1 (optional)  is the bonded type (string)
 223    * Field 2 (optional)  is the atomic number (int)
 224    * Field 3 (mandatory) is the mass (numerical)
 225    * Field 4 (mandatory) is the charge (numerical)
 226    * Field 5 (mandatory) is the particle type (single character)
 227    * This is followed by a number of nonbonded parameters.
 228    *
 229    * The safest way to identify the format is the particle type field.
 230    *
 231    * So, here is what we do:
 232    *
 233    * A. Read in the first six fields as strings
 234    * B. If field 3 (starting from 0) is a single char, we have neither
 235    *    bonded_type or atomic numbers.
 236    * C. If field 5 is a single char we have both.
 237    * D. If field 4 is a single char we check field 1. If this begins with
 238    *    an alphabetical character we have bonded types, otherwise atomic numbers.
 239    */
 240
 241   if(nfields<6)
 242   {
 243       too_few();
 244       return;
 245   }
 246
 247   if( (strlen(tmpfield[5])==1) && isalpha(tmpfield[5][0]) )
 248   {
 249       have_bonded_type   = TRUE;
 250       have_atomic_number = TRUE;
 251   }
 252   else if( (strlen(tmpfield[3])==1) && isalpha(tmpfield[3][0]) )
 253   {
 254       have_bonded_type   = FALSE;
 255       have_atomic_number = FALSE;
 256   }
 257   else
 258   {
 259       have_bonded_type   = ( isalpha(tmpfield[1][0]) != 0 );
 260       have_atomic_number = !have_bonded_type;
 261   }
 262
 263   /* optional fields */
 264   surftens = -1;
 265   vol      =  0;
 266   radius   =  0;
 267   atomnr   = -1;
 268
 269   switch (nb_funct) {
 270
 271   case F_LJ:
 272     nfp0 = 2;
 273
 274     if( have_atomic_number )
 275     {
 276         if ( have_bonded_type )
 277         {
 278             nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf",
 279                            type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],
 280                            &radius,&vol,&surftens);
 281             if(nread < 8)
 282             {
 283                 too_few();
 284                 return;
 285             }
 286         }
 287         else
 288         {
 289             /* have_atomic_number && !have_bonded_type */
 290             nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf",
 291                            type,&atomnr,&m,&q,ptype,&c[0],&c[1],
 292                            &radius,&vol,&surftens);
 293             if(nread < 7)
 294             {
 295                 too_few();
 296                 return;
 297             }
 298         }
 299     }
 300     else
 301     {
 302         if ( have_bonded_type )
 303         {
 304             /* !have_atomic_number && have_bonded_type */
 305             nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf",
 306                            type,btype,&m,&q,ptype,&c[0],&c[1],
 307                            &radius,&vol,&surftens);
 308             if(nread < 7)
 309             {
 310                 too_few();
 311                 return;
 312             }
 313         }
 314         else
 315         {
 316             /* !have_atomic_number && !have_bonded_type */
 317             nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf",
 318                            type,&m,&q,ptype,&c[0],&c[1],
 319                            &radius,&vol,&surftens);
 320             if(nread < 6)
 321             {
 322                 too_few();
 323                 return;
 324             }
 325         }
 326     }
 327
 328     if( !have_bonded_type )
 329     {
 330         strcpy(btype,type);
 331     }
 332
 333     if( !have_atomic_number )
 334     {
 335         atomnr = -1;
 336     }
 337
 338     break;
 339
 340   case F_BHAM:
 341     nfp0 = 3;
 342
 343       if( have_atomic_number )
 344       {
 345           if ( have_bonded_type )
 346           {
 347               nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
 348                              type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
 349                              &radius,&vol,&surftens);
 350               if(nread < 9)
 351               {
 352                   too_few();
 353                   return;
 354               }
 355           }
 356           else
 357           {
 358               /* have_atomic_number && !have_bonded_type */
 359               nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
 360                              type,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
 361                              &radius,&vol,&surftens);
 362               if(nread < 8)
 363               {
 364                   too_few();
 365                   return;
 366               }
 367           }
 368       }
 369       else
 370       {
 371           if ( have_bonded_type )
 372           {
 373               /* !have_atomic_number && have_bonded_type */
 374               nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
 375                              type,btype,&m,&q,ptype,&c[0],&c[1],&c[2],
 376                              &radius,&vol,&surftens);
 377               if(nread < 8)
 378               {
 379                   too_few();
 380                   return;
 381               }
 382           }
 383           else
 384           {
 385               /* !have_atomic_number && !have_bonded_type */
 386               nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
 387                              type,&m,&q,ptype,&c[0],&c[1],&c[2],
 388                              &radius,&vol,&surftens);
 389               if(nread < 7)
 390               {
 391                   too_few();
 392                   return;
 393               }
 394           }
 395       }
 396
 397       if( !have_bonded_type )
 398       {
 399           strcpy(btype,type);
 400       }
 401
 402       if( !have_atomic_number )
 403       {
 404           atomnr = -1;
 405       }
 406
 407       break;
 408
 409       default:
 410           gmx_fatal(FARGS,"Invalid function type %d in push_at %s %d",nb_funct,
 411                     __FILE__,__LINE__);
 412   }
 413   for(j=nfp0; (j<MAXFORCEPARAM); j++)
 414     c[j]=0.0;
 415
 416   if(strlen(type)==1 && isdigit(type[0]))
 417      gmx_fatal(FARGS,"Atom type names can't be single digits.");
 418
 419   if(strlen(btype)==1 && isdigit(btype[0]))
 420      gmx_fatal(FARGS,"Bond atom type names can't be single digits.");
 421
 422   /* Hack to read old topologies */
 423   if (strcasecmp(ptype,"D") == 0)
 424     sprintf(ptype,"V");
 425   for(j=0; (j<eptNR); j++)
 426     if (strcasecmp(ptype,xl[j].entry) == 0)
 427       break;
 428   if (j == eptNR)
 429     gmx_fatal(FARGS,"Invalid particle type %s on line %s",
 430                 ptype,line);
 431   pt=xl[j].ptype;
 432   if (debug)
 433     fprintf(debug,"ptype: %s\n",ptype_str[pt]);
 434
 435   atom->q = q;
 436   atom->m = m;
 437   atom->ptype = pt;
 438   for (i=0; (i<MAXFORCEPARAM); i++)
 439     param->c[i] = c[i];
 440
 441   if ((batype_nr = get_bond_atomtype_type(btype,bat)) == NOTSET)
 442     add_bond_atomtype(bat,symtab,btype);
 443   batype_nr = get_bond_atomtype_type(btype,bat);
 444
 445   if ((nr = get_atomtype_type(type,at)) != NOTSET) {
 446     sprintf(errbuf,"Overriding atomtype %s",type);
 447     warning(errbuf);
 448     if ((nr = set_atomtype(nr,at,symtab,atom,type,param,batype_nr,
 449                            radius,vol,surftens,atomnr)) == NOTSET)
 450       gmx_fatal(FARGS,"Replacing atomtype %s failed",type);
 451   }
 452   else if ((nr = add_atomtype(at,symtab,atom,type,param,
 453                               batype_nr,radius,vol,
 454                               surftens,atomnr)) == NOTSET)
 455     gmx_fatal(FARGS,"Adding atomtype %s failed",type);
 456   else {
 457     /* Add space in the non-bonded parameters matrix */
 458     realloc_nb_params(at,nbparam,pair);
 459   }
 460   sfree(atom);
 461   sfree(param);
 462 }
 463
 464 static void push_bondtype(t_params *       bt,
 465                           t_param *        b,
 466                           int              nral,
 467                           int              ftype,
 468                           bool             bAllowRepeat,
 469                           char *           line)
 470 {
 471         int  i,j;
 472         bool bTest,bFound,bCont,bId;
 473         int  nr   = bt->nr;
 474         int  nrfp = NRFP(ftype);
 475         char  errbuf[256];
 476
 477     /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
 478          are on directly _adjacent_ lines.
 479      */
 480
 481     /* First check if our atomtypes are _identical_ (not reversed) to the previous
 482          entry. If they are not identical we search for earlier duplicates. If they are
 483          we can skip it, since we already searched for the first line
 484          in this group.
 485      */
 486
 487     bFound=FALSE;
 488     bCont =FALSE;
 489
 490     if(bAllowRepeat && nr>1)
 491     {
 492         for (j=0,bCont=TRUE; (j < nral); j++)
 493         {
 494             bCont=bCont && (b->a[j] == bt->param[nr-2].a[j]);
 495         }
 496     }
 497
 498     /* Search for earlier duplicates if this entry was not a continuation
 499          from the previous line.
 500          */
 501     if(!bCont)
 502     {
 503         bFound=FALSE;
 504         for (i=0; (i < nr); i++) {
 505             bTest = TRUE;
 506             for (j=0; (j < nral); j++)
 507                 bTest=(bTest && (b->a[j] == bt->param[i].a[j]));
 508             if (!bTest) {
 509                 bTest=TRUE;
 510                 for(j=0; (j<nral); j++)
 511                     bTest=(bTest && (b->a[nral-1-j] == bt->param[i].a[j]));
 512             }
 513             if (bTest) {
 514                 if (!bFound) {
 515                     bId = TRUE;
 516                     for (j=0; (j < nrfp); j++)
 517                         bId = bId && (bt->param[i].c[j] == b->c[j]);
 518                     if (!bId) {
 519                         sprintf(errbuf,"Overriding %s parameters.%s",
 520                                 interaction_function[ftype].longname,
 521                                 (ftype==F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : "");
 522                         warning(errbuf);
 523                         fprintf(stderr,"  old:");
 524                         for (j=0; (j < nrfp); j++)
 525                             fprintf(stderr," %g",bt->param[i].c[j]);
 526                         fprintf(stderr," \n  new: %s\n\n",line);
 527                     }
 528                 }
 529                 /* Overwrite it! */
 530                 for (j=0; (j < nrfp); j++)
 531                     bt->param[i].c[j] = b->c[j];
 532                 bFound=TRUE;
 533             }
 534         }
 535         if (!bFound) {
 536             /* alloc */
 537             pr_alloc (2,bt);
 538
 539             /* fill the arrays up and down */
 540             memcpy(bt->param[bt->nr].c,  b->c,sizeof(b->c));
 541             memcpy(bt->param[bt->nr].a,  b->a,sizeof(b->a));
 542             memcpy(bt->param[bt->nr+1].c,b->c,sizeof(b->c));
 543
 544             for (j=0; (j < nral); j++)
 545                 bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
 546
 547             bt->nr += 2;
 548         }
 549     }
 550 }
 551
 552 void push_bt(directive d,t_params bt[],int nral,
 553              t_atomtype at,
 554              t_bond_atomtype bat,char *line)
 555 {
 556   const char *formal[MAXATOMLIST+1] = {
 557     "%s",
 558     "%s%s",
 559     "%s%s%s",
 560     "%s%s%s%s",
 561     "%s%s%s%s%s",
 562     "%s%s%s%s%s%s"
 563   };
 564   const char *formnl[MAXATOMLIST+1] = {
 565     "%*s",
 566     "%*s%*s",
 567     "%*s%*s%*s",
 568     "%*s%*s%*s%*s",
 569     "%*s%*s%*s%*s%*s",
 570     "%*s%*s%*s%*s%*s%*s"
 571   };
 572   const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
 573   int      i,ft,ftype,nn,nrfp,nrfpA,nrfpB;
 574   char     f1[STRLEN];
 575   char     alc[MAXATOMLIST+1][20];
 576   /* One force parameter more, so we can check if we read too many */
 577   double   c[MAXFORCEPARAM+1];
 578   t_param  p;
 579   char  errbuf[256];
 580
 581   if ((bat && at) || (!bat && !at))
 582     gmx_incons("You should pass either bat or at to push_bt");
 583
 584   /* Make format string (nral ints+functype) */
 585   if ((nn=sscanf(line,formal[nral],
 586                  alc[0],alc[1],alc[2],alc[3],alc[4],alc[5])) != nral+1) {
 587     sprintf(errbuf,"Not enough atomtypes (%d instead of %d)",nn-1,nral);
 588     warning_error(errbuf);
 589     return;
 590   }
 591
 592   ft    = atoi(alc[nral]);
 593   ftype = ifunc_index(d,ft);
 594   nrfp  = NRFP(ftype);
 595   nrfpA = interaction_function[ftype].nrfpA;
 596   nrfpB = interaction_function[ftype].nrfpB;
 597   strcpy(f1,formnl[nral]);
 598   strcat(f1,formlf);
 599   if ((nn=sscanf(line,f1,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5],&c[6],&c[7],&c[8],&c[9],&c[10],&c[11],&c[12]))
 600       != nrfp) {
 601     if (nn == nrfpA) {
 602       /* Copy the B-state from the A-state */
 603       copy_B_from_A(ftype,c);
 604     } else {
 605       if (nn < nrfpA) {
 606       warning_error("Not enough parameters");
 607       } else if (nn > nrfpA && nn < nrfp) {
 608         warning_error("Too many parameters or not enough parameters for topology B");
 609       } else if (nn > nrfp) {
 610         warning_error("Too many parameters");
 611       }
 612       for(i=nn; (i<nrfp); i++)
 613         c[i] = 0.0;
 614     }
 615   }
 616   for(i=0; (i<nral); i++) {
 617     if (at && ((p.a[i]=get_atomtype_type(alc[i],at)) == NOTSET))
 618       gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
 619     else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
 620       gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
 621   }
 622   for(i=0; (i<nrfp); i++)
 623     p.c[i]=c[i];
 624   push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line);
 625 }
 626
 627
 628 void push_dihedraltype(directive d,t_params bt[],
 629                        t_bond_atomtype bat,char *line)
 630 {
 631   const char *formal[MAXATOMLIST+1] = {
 632     "%s",
 633     "%s%s",
 634     "%s%s%s",
 635     "%s%s%s%s",
 636     "%s%s%s%s%s",
 637     "%s%s%s%s%s%s"
 638   };
 639   const char *formnl[MAXATOMLIST+1] = {
 640     "%*s",
 641     "%*s%*s",
 642     "%*s%*s%*s",
 643     "%*s%*s%*s%*s",
 644     "%*s%*s%*s%*s%*s",
 645     "%*s%*s%*s%*s%*s%*s"
 646   };
 647   const char *formlf[MAXFORCEPARAM] = {
 648     "%lf",
 649     "%lf%lf",
 650     "%lf%lf%lf",
 651     "%lf%lf%lf%lf",
 652     "%lf%lf%lf%lf%lf",
 653     "%lf%lf%lf%lf%lf%lf",
 654     "%lf%lf%lf%lf%lf%lf%lf",
 655     "%lf%lf%lf%lf%lf%lf%lf%lf",
 656     "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 657     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 658     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 659     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 660   };
 661   int      i,ft,ftype,nn,nrfp,nrfpA,nrfpB,nral;
 662   char     f1[STRLEN];
 663   char     alc[MAXATOMLIST+1][20];
 664   double   c[MAXFORCEPARAM];
 665   t_param  p;
 666   bool     bAllowRepeat;
 667   char  errbuf[256];
 668
 669   /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
 670    *
 671    * We first check for 2 atoms with the 3th column being an integer
 672    * defining the type. If this isn't the case, we try it with 4 atoms
 673    * and the 5th column defining the dihedral type.
 674    */
 675   nn=sscanf(line,formal[4],alc[0],alc[1],alc[2],alc[3],alc[4]);
 676   if(nn>=3 && strlen(alc[2])==1 && isdigit(alc[2][0])) {
 677     nral=2;
 678     ft    = atoi(alc[nral]);
 679     /* Move atom types around a bit and use 'X' for wildcard atoms
 680      * to create a 4-atom dihedral definition with arbitrary atoms in
 681      * position 1 and 4.
 682      */
 683     if(alc[2][0]=='2') {
 684       /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
 685       strcpy(alc[3],alc[1]);
 686       sprintf(alc[2],"X");
 687       sprintf(alc[1],"X");
 688       /* alc[0] stays put */
 689     } else {
 690       /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
 691       sprintf(alc[3],"X");
 692       strcpy(alc[2],alc[1]);
 693       strcpy(alc[1],alc[0]);
 694       sprintf(alc[0],"X");
 695     }
 696   } else if(nn==5 && strlen(alc[4])==1 && isdigit(alc[4][0])) {
 697     nral=4;
 698     ft    = atoi(alc[nral]);
 699   } else {
 700     sprintf(errbuf,"Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)",nn-1);
 701     warning_error(errbuf);
 702     return;
 703   }
 704
 705         if(ft == 9)
 706         {
 707                 /* Previously, we have always overwritten parameters if e.g. a torsion
 708                  with the same atomtypes occurs on multiple lines. However, CHARMM and
 709                  some other force fields specify multiple dihedrals over some bonds,
 710                  including cosines with multiplicity 6 and somethimes even higher.
 711                  Thus, they cannot be represented with Ryckaert-Bellemans terms.
 712                  To add support for these force fields, Dihedral type 4 is identical to
 713                  normal proper dihedrals, but repeated entries are allowed.
 714                  */
 715                 bAllowRepeat = TRUE;
 716                 ft = 1;
 717         }
 718         else
 719         {
 720                 bAllowRepeat = FALSE;
 721         }
 722
 723
 724   ftype = ifunc_index(d,ft);
 725   nrfp  = NRFP(ftype);
 726   nrfpA = interaction_function[ftype].nrfpA;
 727   nrfpB = interaction_function[ftype].nrfpB;
 728
 729   strcpy(f1,formnl[nral]);
 730   strcat(f1,formlf[nrfp-1]);
 731
 732   /* Check number of parameters given */
 733   if ((nn=sscanf(line,f1,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5],&c[6],&c[7],&c[8],&c[9],&c[10],&c[11]))
 734       != nrfp) {
 735     if (nn == nrfpA) {
 736       /* Copy the B-state from the A-state */
 737       copy_B_from_A(ftype,c);
 738     } else {
 739       if (nn < nrfpA) {
 740       warning_error("Not enough parameters");
 741       } else if (nn > nrfpA && nn < nrfp) {
 742         warning_error("Too many parameters or not enough parameters for topology B");
 743       } else if (nn > nrfp) {
 744         warning_error("Too many parameters");
 745       }
 746       for(i=nn; (i<nrfp); i++)
 747         c[i] = 0.0;
 748     }
 749   }
 750
 751   for(i=0; (i<4); i++) {
 752     if(!strcmp(alc[i],"X"))
 753       p.a[i]=-1;
 754     else {
 755       if ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET)
 756         gmx_fatal(FARGS,"Unknown bond_atomtype %s",alc[i]);
 757     }
 758   }
 759   for(i=0; (i<nrfp); i++)
 760     p.c[i]=c[i];
 761   /* Always use 4 atoms here, since we created two wildcard atoms
 762    * if there wasn't of them 4 already.
 763    */
 764   push_bondtype (&(bt[ftype]),&p,4,ftype,bAllowRepeat,line);
 765 }
 766
 767
 768 void push_nbt(directive d,t_nbparam **nbt,t_atomtype atype,
 769               char *pline,int nb_funct)
 770 {
 771   /* swap the atoms */
 772   const char *form2="%*s%*s%*s%lf%lf";
 773   const char *form3="%*s%*s%*s%lf%lf%lf";
 774   const char *form4="%*s%*s%*s%lf%lf%lf%lf";
 775   char    a0[80],a1[80];
 776   int     i,f,n,ftype,atnr,nrfp;
 777   double  c[4];
 778   real    cr[4],sig6;
 779   atom_id ai,aj;
 780   t_nbparam *nbp;
 781   bool    bId;
 782   char  errbuf[256];
 783
 784   if (sscanf (pline,"%s%s%d",a0,a1,&f) != 3) {
 785     too_few();
 786     return;
 787   }
 788
 789   ftype=ifunc_index(d,f);
 790
 791   if (ftype != nb_funct) {
 792     sprintf(errbuf,"Trying to add %s while the default nonbond type is %s",
 793             interaction_function[ftype].longname,
 794             interaction_function[nb_funct].longname);
 795     warning_error(errbuf);
 796     return;
 797   }
 798
 799   /* Get the force parameters */
 800   nrfp = NRFP(ftype);
 801   if (ftype == F_LJ14) {
 802     n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
 803     if (n < 2) {
 804       too_few();
 805       return;
 806     }
 807     /* When the B topology parameters are not set,
 808      * copy them from topology A
 809      */
 810     for(i=n; i<nrfp; i++)
 811       c[i] = c[i-2];
 812   }
 813   else if (ftype == F_LJC14_Q) {
 814     n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
 815     if (n < 4) {
 816       too_few();
 817       return;
 818     }
 819   }
 820   else if (nrfp == 2) {
 821     if (sscanf(pline,form2,&c[0],&c[1]) != 2) {
 822       too_few();
 823       return;
 824     }
 825   }
 826   else if (nrfp == 3) {
 827     if (sscanf(pline,form3,&c[0],&c[1],&c[2]) != 3) {
 828       too_few();
 829       return;
 830     }
 831   }
 832   else {
 833     gmx_fatal(FARGS,"Number of force parameters for nonbonded interactions is %d"
 834                 " in file %s, line %d",nrfp,__FILE__,__LINE__);
 835   }
 836   for(i=0; (i<nrfp); i++)
 837     cr[i] = c[i];
 838
 839   /* Put the parameters in the matrix */
 840   if ((ai = get_atomtype_type (a0,atype)) == NOTSET)
 841     gmx_fatal(FARGS,"Atomtype %s not found",a0);
 842   if ((aj = get_atomtype_type (a1,atype)) == NOTSET)
 843     gmx_fatal(FARGS,"Atomtype %s not found",a1);
 844   nbp = &(nbt[max(ai,aj)][min(ai,aj)]);
 845
 846   if (nbp->bSet) {
 847     bId = TRUE;
 848     for (i=0; i<nrfp; i++)
 849       bId = bId && (nbp->c[i] == cr[i]);
 850     if (!bId) {
 851       sprintf(errbuf,"Overriding non-bonded parameters,");
 852       warning(errbuf);
 853       fprintf(stderr,"  old:");
 854       for (i=0; i<nrfp; i++)
 855         fprintf(stderr," %g",nbp->c[i]);
 856       fprintf(stderr," new\n%s\n",pline);
 857     }
 858   }
 859   nbp->bSet = TRUE;
 860   for (i=0; i<nrfp; i++)
 861     nbp->c[i] = cr[i];
 862 }
 863
 864 static void push_atom_now(t_symtab *symtab,t_atoms *at,int atomnr,
 865                           int atomicnumber,
 866                           int type,char *ctype,int ptype,
 867                           int resnumber,int cgnumber,
 868                           char *resname,char *name,real m0,real q0,
 869                           int typeB,char *ctypeB,real mB,real qB)
 870 {
 871   int j,resnr_diff=-1;
 872   int nr = at->nr;
 873
 874   if (((nr==0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
 875     gmx_fatal(FARGS,"Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)",atomnr,at->nr);
 876   if (nr)
 877     resnr_diff = resnumber - (at->atom[at->nr-1].resnr+1);
 878   if (((nr==0) && (resnumber != 1)) ||
 879       (nr && (resnr_diff != 0) && (resnr_diff != 1)))
 880     gmx_fatal(FARGS,"Residue numbers in the .top are not numbered consecutively from 1 (rather, resnumber = %d and resnr_diff = %d)",resnumber,resnr_diff);
 881
 882   /* New atom instance
 883    * get new space for arrays
 884    */
 885   srenew(at->atom,nr+1);
 886   srenew(at->atomname,nr+1);
 887   srenew(at->atomtype,nr+1);
 888   srenew(at->atomtypeB,nr+1);
 889
 890   if (resnumber > at->nres) {
 891     at->nres=resnumber;
 892     srenew(at->resname,resnumber);
 893     at->resname[resnumber-1] = put_symtab(symtab,resname);
 894   }
 895   /* fill the list */
 896   at->atom[nr].type  = type;
 897   at->atom[nr].ptype = ptype;
 898   at->atom[nr].q     = q0;
 899   at->atom[nr].m     = m0;
 900   at->atom[nr].typeB = typeB;
 901   at->atom[nr].qB    = qB;
 902   at->atom[nr].mB    = mB;
 903
 904   at->atom[nr].resnr = resnumber-1;
 905   at->atom[nr].atomnumber = atomicnumber;
 906   at->atomname[nr] = put_symtab(symtab,name);
 907   at->atomtype[nr] = put_symtab(symtab,ctype);
 908   at->atomtypeB[nr] = put_symtab(symtab,ctypeB);
 909   at->nr++;
 910 }
 911
 912 void push_cg(t_block *block, int *lastindex, int index, int a)
 913 {
 914   if (debug)
 915     fprintf (debug,"Index %d, Atom %d\n",index,a);
 916
 917   if (((block->nr) && (*lastindex != index)) || (!block->nr)) {
 918     /* add a new block */
 919     block->nr++;
 920     srenew(block->index,block->nr+1);
 921   }
 922   block->index[block->nr] = a + 1;
 923   *lastindex = index;
 924 }
 925
 926 void push_atom(t_symtab *symtab,t_block *cgs,
 927                t_atoms *at,t_atomtype atype,char *line,int *lastcg)
 928 {
 929   int           nr,ptype;
 930   int           resnumber,cgnumber,atomnr,type,typeB,nscan;
 931   char          id[STRLEN],ctype[STRLEN],ctypeB[STRLEN],
 932                 resname[STRLEN],name[STRLEN],check[STRLEN];
 933   double        m,q,mb,qb;
 934   real          m0,q0,mB,qB;
 935
 936   /* Make a shortcut for writing in this molecule  */
 937   nr = at->nr;
 938
 939   /* Fixed parameters */
 940   if (sscanf(line,"%s%s%d%s%s%d",
 941              id,ctype,&resnumber,resname,name,&cgnumber) != 6) {
 942     too_few();
 943     return;
 944   }
 945   sscanf(id,"%d",&atomnr);
 946   if ((type  = get_atomtype_type(ctype,atype)) == NOTSET)
 947     gmx_fatal(FARGS,"Atomtype %s not found",ctype);
 948   ptype = get_atomtype_ptype(type,atype);
 949
 950   /* Set default from type */
 951   q0    = get_atomtype_qA(type,atype);
 952   m0    = get_atomtype_massA(type,atype);
 953   typeB = type;
 954   qB    = q0;
 955   mB    = m0;
 956
 957   /* Optional parameters */
 958   nscan = sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
 959                  &q,&m,ctypeB,&qb,&mb,check);
 960
 961   /* Nasty switch that falls thru all the way down! */
 962   if (nscan > 0) {
 963     q0 = qB = q;
 964     if (nscan > 1) {
 965       m0 = mB = m;
 966       if (nscan > 2) {
 967         typeB=get_atomtype_type(ctypeB,atype);
 968         qB = get_atomtype_qA(typeB,atype);
 969         mB = get_atomtype_massA(typeB,atype);
 970         if (nscan > 3) {
 971           qB = qb;
 972           if (nscan > 4) {
 973             mB = mb;
 974             if (nscan > 5)
 975               warning_error("Too many parameters");
 976           }
 977         }
 978       }
 979     }
 980   }
 981   if (debug)
 982     fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
 983
 984   push_cg(cgs,lastcg,cgnumber,nr);
 985
 986   push_atom_now(symtab,at,atomnr,get_atomtype_atomnumber(type,atype),
 987                 type,ctype,ptype,resnumber,cgnumber,
 988                 resname,name,m0,q0,typeB,
 989                 typeB==type ? ctype : ctypeB,mB,qB);
 990 }
 991
 992 void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line)
 993 {
 994   char type[STRLEN];
 995   int nrexcl,i;
 996   t_molinfo *newmol;
 997
 998   if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
 999     warning_error("Expected a molecule type name and nrexcl");
1000   }
1001
1002   /* Test if this atomtype overwrites another */
1003   i = 0;
1004   while (i < *nmol) {
1005     if (strcasecmp(*((*mol)[i].name),type) == 0)
1006       gmx_fatal(FARGS,"moleculetype %s is redefined",type);
1007     i++;
1008   }
1009
1010   (*nmol)++;
1011   srenew(*mol,*nmol);
1012   newmol = &((*mol)[*nmol-1]);
1013   init_molinfo(newmol);
1014
1015   /* Fill in the values */
1016   newmol->name     = put_symtab(symtab,type);
1017   newmol->nrexcl   = nrexcl;
1018   newmol->excl_set = FALSE;
1019 }
1020
1021 static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
1022                               t_param *p,int c_start,bool bB,bool bGenPairs)
1023 {
1024   int      i,j,ti,tj,ntype;
1025   bool     bFound;
1026   t_param  *pi=NULL;
1027   int      nr   = bt[ftype].nr;
1028   int      nral = NRAL(ftype);
1029   int      nrfp = interaction_function[ftype].nrfpA;
1030   int      nrfpB= interaction_function[ftype].nrfpB;
1031
1032   if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1033     return TRUE;
1034
1035   bFound=FALSE;
1036   if(bGenPairs) {
1037     /* First test the generated-pair position to save
1038      * time when we have 1000*1000 entries for e.g. OPLS...
1039      */
1040     ntype=sqrt(nr);
1041     if (bB) {
1042       ti=at->atom[p->a[0]].typeB;
1043       tj=at->atom[p->a[1]].typeB;
1044     } else {
1045       ti=at->atom[p->a[0]].type;
1046       tj=at->atom[p->a[1]].type;
1047     }
1048     pi=&(bt[ftype].param[ntype*ti+tj]);
1049     bFound=((ti == pi->a[0]) && (tj == pi->a[1]));
1050   }
1051
1052   /* Search explicitly if we didnt find it */
1053   if(!bFound) {
1054     for (i=0; ((i < nr) && !bFound); i++) {
1055       pi=&(bt[ftype].param[i]);
1056       if (bB)
1057         for (j=0; ((j < nral) &&
1058                    (at->atom[p->a[j]].typeB == pi->a[j])); j++);
1059       else
1060         for (j=0; ((j < nral) &&
1061                    (at->atom[p->a[j]].type == pi->a[j])); j++);
1062       bFound=(j==nral);
1063     }
1064   }
1065
1066   if (bFound) {
1067     if (bB) {
1068       if (nrfp+nrfpB > MAXFORCEPARAM) {
1069         gmx_incons("Too many force parameters");
1070       }
1071       for (j=c_start; (j < nrfpB); j++)
1072         p->c[nrfp+j] = pi->c[j];
1073     }
1074     else
1075       for (j=c_start; (j < nrfp); j++)
1076         p->c[j] = pi->c[j];
1077   }
1078   else {
1079     for (j=c_start; (j < nrfp); j++)
1080       p->c[j] = 0.0;
1081   }
1082   return bFound;
1083 }
1084
1085
1086 static bool default_params(int ftype,t_params bt[],
1087                            t_atoms *at,t_atomtype atype,
1088                            t_param *p,bool bB,
1089                            t_param **param_def,
1090                int *nparam_def)
1091 {
1092   int      i,j,nparam_found;
1093   bool     bFound,bSame;
1094   t_param  *pi=NULL;
1095   t_param  *pj=NULL;
1096   int      nr   = bt[ftype].nr;
1097   int      nral = NRAL(ftype);
1098   int      nrfpA= interaction_function[ftype].nrfpA;
1099   int      nrfpB= interaction_function[ftype].nrfpB;
1100
1101   if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1102     return TRUE;
1103
1104
1105   /* We allow wildcards now. The first type (with or without wildcards) that
1106    * fits is used, so you should probably put the wildcarded bondtypes
1107    * at the end of each section.
1108    */
1109   bFound=FALSE;
1110   nparam_found=0;
1111   /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1112    * special case for this. Check for B state outside loop to speed it up.
1113    */
1114   if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS)
1115   {
1116           if (bB)
1117           {
1118                   for (i=0; ((i < nr) && !bFound); i++)
1119                   {
1120                           pi=&(bt[ftype].param[i]);
1121                           bFound=
1122                           (
1123                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1124                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1125                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1126                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1127                            );
1128                   }
1129           }
1130           else
1131           {
1132                   /* State A */
1133                   for (i=0; ((i < nr) && !bFound); i++)
1134                   {
1135                           pi=&(bt[ftype].param[i]);
1136                           bFound=
1137                           (
1138                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].type,atype)==pi->AI)) &&
1139                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].type,atype)==pi->AJ)) &&
1140                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].type,atype)==pi->AK)) &&
1141                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].type,atype)==pi->AL))
1142                            );
1143                   }
1144           }
1145           /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1146            * The rules in that case is that additional matches HAVE to be on adjacent lines!
1147            */
1148           if(bFound==TRUE)
1149           {
1150                   nparam_found++;
1151                   bSame = TRUE;
1152                   /* Continue from current i value */
1153                   for(j=i ; j<nr && bSame ; j++)
1154                   {
1155                           pj=&(bt[ftype].param[j]);
1156                           bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1157                           if(bSame)
1158                                   nparam_found++;
1159                           /* nparam_found will be increased as long as the numbers match */
1160                   }
1161           }
1162   } else { /* Not a dihedral */
1163           for (i=0; ((i < nr) && !bFound); i++) {
1164                   pi=&(bt[ftype].param[i]);
1165                   if (bB)
1166                           for (j=0; ((j < nral) &&
1167                                                  (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
1168                                   ;
1169                   else
1170                           for (j=0; ((j < nral) &&
1171                                                  (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
1172                                   ;
1173                   bFound=(j==nral);
1174           }
1175           if(bFound)
1176                   nparam_found=1;
1177   }
1178
1179   *param_def = pi;
1180   *nparam_def = nparam_found;
1181
1182   return bFound;
1183 }
1184
1185
1186
1187 void push_bond(directive d,t_params bondtype[],t_params bond[],
1188                t_atoms *at,t_atomtype atype,char *line,
1189                bool bBonded,bool bGenPairs,real fudgeQQ,
1190                bool bZero,bool *bWarn_copy_A_B)
1191 {
1192   const char *aaformat[MAXATOMLIST]= {
1193     "%d%d",
1194     "%d%d%d",
1195     "%d%d%d%d",
1196     "%d%d%d%d%d",
1197     "%d%d%d%d%d%d"
1198   };
1199   const char *asformat[MAXATOMLIST]= {
1200     "%*s%*s",
1201     "%*s%*s%*s",
1202     "%*s%*s%*s%*s",
1203     "%*s%*s%*s%*s%*s",
1204     "%*s%*s%*s%*s%*s%*s"
1205   };
1206   const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1207   int      nr,i,j,nral,nread,ftype;
1208   char     format[STRLEN];
1209   /* One force parameter more, so we can check if we read too many */
1210   double   cc[MAXFORCEPARAM+1];
1211   int      aa[MAXATOMLIST+1];
1212   t_param  param,paramB,*param_defA,*param_defB;
1213   bool     bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
1214   int      nparam_defA,nparam_defB;
1215   char  errbuf[256];
1216
1217   nparam_defA=nparam_defB=0;
1218
1219   ftype = ifunc_index(d,1);
1220   nral  = NRAL(ftype);
1221   for(j=0; j<MAXATOMLIST; j++)
1222     aa[j]=NOTSET;
1223   bDef = (NRFP(ftype) > 0);
1224
1225   nread = sscanf(line,aaformat[nral-1],
1226                  &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1227   if (nread < nral) {
1228     too_few();
1229     return;
1230   } else if (nread == nral)
1231     ftype = ifunc_index(d,1);
1232   else {
1233     /* this is a hack to allow for virtual sites with swapped parity */
1234     bSwapParity = (aa[nral]<0);
1235     if (bSwapParity)
1236       aa[nral] = -aa[nral];
1237     ftype = ifunc_index(d,aa[nral]);
1238     if (bSwapParity)
1239       switch(ftype) {
1240       case F_VSITE3FAD:
1241       case F_VSITE3OUT:
1242         break;
1243       default:
1244         gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1245                     interaction_function[F_VSITE3FAD].longname,
1246                     interaction_function[F_VSITE3OUT].longname);
1247       }
1248   }
1249
1250   /* Check for double atoms and atoms out of bounds */
1251   for(i=0; (i<nral); i++) {
1252     if ( aa[i] < 1 || aa[i] > at->nr )
1253       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1254                 "Atom index (%d) in %s out of bounds (1-%d).\n"
1255                 "This probably means that you have inserted topology section \"%s\"\n"
1256                 "in a part belonging to a different molecule than you intended to.\n"
1257                 "In that case move the \"%s\" section to the right molecule.",
1258                 get_warning_file(),get_warning_line(),
1259                 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1260     for(j=i+1; (j<nral); j++)
1261       if (aa[i] == aa[j]) {
1262         sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1263         warning(errbuf);
1264       }
1265   }
1266   if (ftype == F_SETTLE)
1267     if (aa[0]+2 > at->nr)
1268       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1269                   "             Atom index (%d) in %s out of bounds (1-%d)\n"
1270                   "             Settle works on atoms %d, %d and %d",
1271                   get_warning_file(),get_warning_line(),
1272                   aa[0],dir2str(d),at->nr,
1273                   aa[0],aa[0]+1,aa[0]+2);
1274
1275   /* default force parameters  */
1276   for(j=0; (j<MAXATOMLIST); j++)
1277     param.a[j] = aa[j]-1;
1278   for(j=0; (j<MAXFORCEPARAM); j++)
1279     param.c[j] = 0.0;
1280
1281   /* Get force params for normal and free energy perturbation
1282    * studies, as determined by types!
1283    */
1284   if (bBonded) {
1285     bFoundA = default_params(ftype,bondtype,at,atype,&param,FALSE,&param_defA,&nparam_defA);
1286     if (bFoundA) {
1287       /* Copy the A-state and B-state default parameters */
1288       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1289         param.c[j] = param_defA->c[j];
1290     }
1291     bFoundB = default_params(ftype,bondtype,at,atype,&param,TRUE,&param_defB,&nparam_defB);
1292     if (bFoundB) {
1293       /* Copy only the B-state default parameters */
1294       for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1295         param.c[j] = param_defB->c[j];
1296     }
1297   } else if (ftype == F_LJ14) {
1298     bFoundA = default_nb_params(ftype, bondtype,at,&param,0,FALSE,bGenPairs);
1299     bFoundB = default_nb_params(ftype, bondtype,at,&param,0,TRUE, bGenPairs);
1300   } else if (ftype == F_LJC14_Q) {
1301     param.c[0] = fudgeQQ;
1302     /* Fill in the A-state charges as default parameters */
1303     param.c[1] = at->atom[param.a[0]].q;
1304     param.c[2] = at->atom[param.a[1]].q;
1305     /* The default LJ parameters are the standard 1-4 parameters */
1306     bFoundA = default_nb_params(F_LJ14,bondtype,at,&param,3,FALSE,bGenPairs);
1307     bFoundB = TRUE;
1308   } else if (ftype == F_LJC_PAIRS_NB) {
1309     /* Defaults are not supported here */
1310     bFoundA = FALSE;
1311     bFoundB = TRUE;
1312   } else {
1313     gmx_incons("Unknown function type in push_bond");
1314   }
1315
1316   if (nread > nral) {
1317           /* Manually specified parameters - in this case we discard multiple torsion info! */
1318
1319     strcpy(format,asformat[nral-1]);
1320     strcat(format,ccformat);
1321
1322     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1323                    &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1324
1325     if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1326       {
1327         /* We only have to issue a warning if these atoms are perturbed! */
1328         bPert = FALSE;
1329         for(j=0; (j<nral); j++)
1330           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1331
1332         if (bPert && *bWarn_copy_A_B)
1333           {
1334             sprintf(errbuf,
1335                     "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1336             warning(errbuf);
1337             *bWarn_copy_A_B = FALSE;
1338           }
1339
1340         /* If only the A parameters were specified, copy them to the B state */
1341         /* The B-state parameters correspond to the first nrfpB
1342          * A-state parameters.
1343          */
1344         for(j=0; (j<NRFPB(ftype)); j++)
1345           cc[nread++] = cc[j];
1346       }
1347
1348     /* If nread was 0 or EOF, no parameters were read => use defaults.
1349      * If nread was nrfpA we copied above so nread=nrfp.
1350      * If nread was nrfp we are cool.
1351      * For F_LJC14_Q we allow supplying fudgeQQ only.
1352      * Anything else is an error!
1353      */
1354     if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1355         !(ftype == F_LJC14_Q && nread == 1))
1356       {
1357         gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1358                   nread,NRFPA(ftype),NRFP(ftype),
1359                   interaction_function[ftype].longname);
1360       }
1361
1362     for(j=0; (j<nread); j++)
1363       param.c[j]=cc[j];
1364
1365     /* Check whether we have to use the defaults */
1366     if (nread == NRFP(ftype))
1367       bDef = FALSE;
1368   }
1369   else
1370   {
1371     nread = 0;
1372   }
1373   /* nread now holds the number of force parameters read! */
1374
1375
1376
1377
1378         if (bDef)
1379         {
1380           /* Use defaults */
1381                 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1382                 if(ftype==F_PDIHS)
1383                 {
1384                         if ((nparam_defA != nparam_defB) || ((nparam_defA>1 || nparam_defB>1) && (param_defA!=param_defB)))
1385                         {
1386                                 sprintf(errbuf,
1387                                                 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1388                                                 "Please specify perturbed parameters manually for this torsion in your topology!");
1389                                 warning_error(errbuf);
1390                         }
1391                 }
1392
1393           if (nread > 0 && nread < NRFPA(ftype))
1394           {
1395                   /* Issue an error, do not use defaults */
1396                   sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1397                   warning_error(errbuf);
1398           }
1399
1400       if (nread == 0 || nread == EOF)
1401           {
1402                   if (!bFoundA)
1403                   {
1404                           if (interaction_function[ftype].flags & IF_VSITE)
1405                           {
1406                                   /* set them to NOTSET, will be calculated later */
1407                                   for(j=0; (j<MAXFORCEPARAM); j++)
1408                                           param.c[j] = NOTSET;
1409
1410                                   if (bSwapParity)
1411                                           param.C1 = -1; /* flag to swap parity of vsite construction */
1412                           }
1413                           else
1414                           {
1415                                   if (bZero)
1416                                   {
1417                                           fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1418                                                           interaction_function[ftype].longname);
1419                                   }
1420                                   else
1421                                   {
1422                                           sprintf(errbuf,"No default %s types",interaction_function[ftype].longname);
1423                                           warning_error(errbuf);
1424                                   }
1425                           }
1426                   }
1427                   else
1428                   {
1429                           if (bSwapParity)
1430                           {
1431                                   switch(ftype)
1432                                   {
1433                                           case F_VSITE3FAD:
1434                                                   param.C0 = 360 - param.C0;
1435                                                   break;
1436                                           case F_VSITE3OUT:
1437                                                   param.C2 = -param.C2;
1438                                                   break;
1439                                   }
1440                           }
1441                   }
1442                   if (!bFoundB)
1443                   {
1444                           /* We only have to issue a warning if these atoms are perturbed! */
1445                           bPert = FALSE;
1446                           for(j=0; (j<nral); j++)
1447                                   bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1448
1449                           if (bPert)
1450                           {
1451                                   sprintf(errbuf,"No default %s types for perturbed atoms, "
1452                                                   "using normal values",interaction_function[ftype].longname);
1453                                   warning(errbuf);
1454                           }
1455                   }
1456           }
1457   }
1458
1459         if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1460                 && param.c[5]!=param.c[2])
1461                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1462                                   "             %s multiplicity can not be perturbed %f!=%f",
1463                                   get_warning_file(),get_warning_line(),
1464                                   interaction_function[ftype].longname,
1465                                   param.c[2],param.c[5]);
1466
1467         if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1468                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1469                                   "             %s table number can not be perturbed %d!=%d",
1470                                   get_warning_file(),get_warning_line(),
1471                                   interaction_function[ftype].longname,
1472                                   (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1473
1474         /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1475         if (ftype==F_RBDIHS) {
1476                 nr=0;
1477                 for(i=0;i<NRFP(ftype);i++) {
1478                         if(param.c[i]!=0)
1479                                 nr++;
1480                 }
1481                 if(nr==0)
1482                         return;
1483         }
1484
1485         /* Put the values in the appropriate arrays */
1486         push_bondnow (&bond[ftype],&param);
1487
1488   if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1489     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1490                 "             %s table number can not be perturbed %d!=%d",
1491                 get_warning_file(),get_warning_line(),
1492                 interaction_function[ftype].longname,
1493                 (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1494
1495   /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1496   if (ftype==F_RBDIHS) {
1497     nr=0;
1498     for(i=0;i<NRFP(ftype);i++) {
1499       if(param.c[i]!=0)
1500         nr++;
1501     }
1502     if(nr==0)
1503       return;
1504   }
1505
1506   /* Put the values in the appropriate arrays */
1507   add_param_to_list (&bond[ftype],&param);
1508
1509   /* Push additional torsions from FF for ftype==9 if we have them.
1510    * We have already checked that the A/B states do not differ in this case,
1511    * so we do not have to double-check that again, or the vsite stuff.
1512    * In addition, those torsions cannot be automatically perturbed.
1513    */
1514   if(bDef && ftype==F_PDIHS)
1515   {
1516     for(i=1;i<nparam_defA;i++)
1517     {
1518       /* Advance pointer! */
1519       param_defA++;
1520       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1521         param.c[j] = param_defA->c[j];
1522       /* And push the next term for this torsion */
1523       push_bondnow (&bond[ftype],&param);
1524     }
1525   }
1526 }
1527
1528
1529
1530 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1531                   t_atoms *at,t_atomtype atype,char *line)
1532 {
1533   char *ptr;
1534   int  type,ftype,j,n,ret,nj,a;
1535   int  *atc=NULL;
1536   double *weight=NULL,weight_tot;
1537   t_param param;
1538
1539   /* default force parameters  */
1540   for(j=0; (j<MAXATOMLIST); j++)
1541     param.a[j] = NOTSET;
1542   for(j=0; (j<MAXFORCEPARAM); j++)
1543     param.c[j] = 0.0;
1544
1545   ptr = line;
1546   ret = sscanf(ptr,"%d%n",&a,&n);
1547   ptr += n;
1548   if (ret == 0)
1549     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1550               "             Expected an atom index in section \"%s\"",
1551                 get_warning_file(),get_warning_line(),
1552               dir2str(d));
1553
1554   param.a[0] = a - 1;
1555
1556   ret = sscanf(ptr,"%d%n",&type,&n);
1557   ptr += n;
1558   ftype = ifunc_index(d,type);
1559
1560   weight_tot = 0;
1561   nj = 0;
1562   do {
1563     ret = sscanf(ptr,"%d%n",&a,&n);
1564     ptr += n;
1565     if (ret > 0) {
1566       if (nj % 20 == 0) {
1567         srenew(atc,nj+20);
1568         srenew(weight,nj+20);
1569       }
1570       atc[nj] = a - 1;
1571       switch (type) {
1572       case 1:
1573         weight[nj] = 1;
1574         break;
1575       case 2:
1576         /* Here we use the A-state mass as a parameter.
1577          * Note that the B-state mass has no influence.
1578          */
1579         weight[nj] = at->atom[atc[nj]].m;
1580         break;
1581       case 3:
1582         weight[nj] = -1;
1583         ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1584         ptr += n;
1585         if (weight[nj] < 0)
1586           gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1587                     "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1588                     get_warning_file(),get_warning_line(),
1589                     nj+1,atc[nj]+1);
1590         break;
1591       default:
1592         gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1593       }
1594       weight_tot += weight[nj];
1595       nj++;
1596     }
1597   } while (ret > 0);
1598
1599   if (nj == 0)
1600     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1601               "             Expected more than one atom index in section \"%s\"",
1602               get_warning_file(),get_warning_line(),
1603               dir2str(d));
1604
1605   if (weight_tot == 0)
1606      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1607                "             The total mass of the construting atoms is zero",
1608                get_warning_file(),get_warning_line());
1609
1610   for(j=0; j<nj; j++) {
1611     param.a[1] = atc[j];
1612     param.c[0] = nj;
1613     param.c[1] = weight[j]/weight_tot;
1614     /* Put the values in the appropriate arrays */
1615     add_param_to_list (&bond[ftype],&param);
1616   }
1617
1618   sfree(atc);
1619   sfree(weight);
1620 }
1621
1622 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1623               int *nrcopies)
1624 {
1625   int  i,copies;
1626   char type[STRLEN];
1627
1628   *nrcopies=0;
1629   if (sscanf(pline,"%s%d",type,&copies) != 2) {
1630     too_few();
1631     return;
1632   }
1633
1634   /* search moleculename */
1635   for (i=0; ((i<nrmols) && strcasecmp(type,*(mols[i].name))); i++)
1636     ;
1637
1638   if (i<nrmols) {
1639     *nrcopies        = copies;
1640     *whichmol        = i;
1641   } else
1642     gmx_fatal(FARGS,"No such moleculetype %s",type);
1643 }
1644
1645 void init_block2(t_block2 *b2, int natom)
1646 {
1647   int i;
1648
1649   b2->nr=natom;
1650   snew(b2->nra,b2->nr);
1651   snew(b2->a,b2->nr);
1652   for(i=0; (i<b2->nr); i++)
1653     b2->a[i]=NULL;
1654 }
1655
1656 void done_block2(t_block2 *b2)
1657 {
1658   int i;
1659
1660   if (b2->nr) {
1661     for(i=0; (i<b2->nr); i++)
1662       sfree(b2->a[i]);
1663     sfree(b2->a);
1664     sfree(b2->nra);
1665     b2->nr=0;
1666   }
1667 }
1668
1669 void push_excl(char *line, t_block2 *b2)
1670 {
1671   int  i,j;
1672   int  n;
1673   char base[STRLEN],format[STRLEN];
1674
1675   if (sscanf(line,"%d",&i)==0)
1676     return;
1677
1678   if ((1 <= i) && (i <= b2->nr))
1679     i--;
1680   else {
1681     if (debug)
1682       fprintf(debug,"Unbound atom %d\n",i-1);
1683     return;
1684   }
1685   strcpy(base,"%*d");
1686   do {
1687     strcpy(format,base);
1688     strcat(format,"%d");
1689     n=sscanf(line,format,&j);
1690     if (n == 1) {
1691       if ((1 <= j) && (j <= b2->nr)) {
1692         j--;
1693         srenew(b2->a[i],++(b2->nra[i]));
1694         b2->a[i][b2->nra[i]-1]=j;
1695         /* also add the reverse exclusion! */
1696         srenew(b2->a[j],++(b2->nra[j]));
1697         b2->a[j][b2->nra[j]-1]=i;
1698         strcat(base,"%*d");
1699       }
1700       else
1701         gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
1702     }
1703   } while (n == 1);
1704 }
1705
1706 void b_to_b2(t_blocka *b, t_block2 *b2)
1707 {
1708   int     i;
1709   atom_id j,a;
1710
1711   for(i=0; (i<b->nr); i++)
1712     for(j=b->index[i]; (j<b->index[i+1]); j++) {
1713       a=b->a[j];
1714       srenew(b2->a[i],++b2->nra[i]);
1715       b2->a[i][b2->nra[i]-1]=a;
1716     }
1717 }
1718
1719 void b2_to_b(t_block2 *b2, t_blocka *b)
1720 {
1721   int     i,nra;
1722   atom_id j;
1723
1724   nra=0;
1725   for(i=0; (i<b2->nr); i++) {
1726     b->index[i]=nra;
1727     for(j=0; (j<b2->nra[i]); j++)
1728       b->a[nra+j]=b2->a[i][j];
1729     nra+=b2->nra[i];
1730   }
1731   /* terminate list */
1732   b->index[i]=nra;
1733 }
1734
1735 static int icomp(const void *v1, const void *v2)
1736 {
1737   return (*((atom_id *) v1))-(*((atom_id *) v2));
1738 }
1739
1740 void merge_excl(t_blocka *excl, t_block2 *b2)
1741 {
1742   int     i,k;
1743   atom_id j;
1744   int     nra;
1745
1746   if (!b2->nr)
1747     return;
1748   else if (b2->nr != excl->nr) {
1749     gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
1750                 b2->nr,excl->nr);
1751   }
1752   else if (debug)
1753     fprintf(debug,"Entering merge_excl\n");
1754
1755   /* First copy all entries from excl to b2 */
1756   b_to_b2(excl,b2);
1757
1758   /* Count and sort the exclusions */
1759   nra=0;
1760   for(i=0; (i<b2->nr); i++) {
1761     if (b2->nra[i] > 0) {
1762       /* remove double entries */
1763       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
1764       k=1;
1765       for(j=1; (j<b2->nra[i]); j++)
1766         if (b2->a[i][j]!=b2->a[i][k-1]) {
1767           b2->a[i][k]=b2->a[i][j];
1768           k++;
1769         }
1770       b2->nra[i]=k;
1771       nra+=b2->nra[i];
1772     }
1773   }
1774   excl->nra=nra;
1775   srenew(excl->a,excl->nra);
1776
1777   b2_to_b(b2,excl);
1778 }
1779
1780 int add_atomtype_decoupled(t_symtab *symtab,t_atomtype at,
1781                            t_nbparam ***nbparam,t_nbparam ***pair)
1782 {
1783   t_atom  atom;
1784   t_param param;
1785   int     i,nr;
1786
1787   /* Define an atom type with all parameters set to zero (no interactions) */
1788   atom.q = 0.0;
1789   atom.m = 0.0;
1790   /* Type for decoupled atoms could be anything,
1791    * this should be changed automatically later when required.
1792    */
1793   atom.ptype = eptAtom;
1794   for (i=0; (i<MAXFORCEPARAM); i++)
1795     param.c[i] = 0.0;
1796
1797   nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0);
1798
1799   /* Add space in the non-bonded parameters matrix */
1800   realloc_nb_params(at,nbparam,pair);
1801
1802   return nr;
1803 }
1804
1805 static void convert_pairs_to_pairsQ(t_params *plist,
1806                                     real fudgeQQ,t_atoms *atoms)
1807 {
1808   t_param *param;
1809   int  i;
1810   real v,w;
1811
1812   /* Copy the pair list to the pairQ list */
1813   plist[F_LJC14_Q] = plist[F_LJ14];
1814   /* Empty the pair list */
1815   plist[F_LJ14].nr    = 0;
1816   plist[F_LJ14].param = NULL;
1817   param = plist[F_LJC14_Q].param;
1818   for(i=0; i<plist[F_LJC14_Q].nr; i++) {
1819     v = param[i].c[0];
1820     w = param[i].c[1];
1821     param[i].c[0] = fudgeQQ;
1822     param[i].c[1] = atoms->atom[param[i].a[0]].q;
1823     param[i].c[2] = atoms->atom[param[i].a[1]].q;
1824     param[i].c[3] = v;
1825     param[i].c[4] = w;
1826   }
1827 }
1828
1829 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
1830 {
1831   int n,ntype,i,j,k;
1832   t_atom *atom;
1833   t_blocka *excl;
1834   bool bExcl;
1835   t_param param;
1836
1837   n = mol->atoms.nr;
1838   atom = mol->atoms.atom;
1839
1840   ntype = sqrt(nbp->nr);
1841
1842   for (i=0; i<MAXATOMLIST; i++)
1843     param.a[i] = NOTSET;
1844   for (i=0; i<MAXFORCEPARAM; i++)
1845     param.c[i] = NOTSET;
1846
1847   /* Add a pair interaction for all non-excluded atom pairs */
1848   excl = &mol->excls;
1849   for(i=0; i<n; i++) {
1850     for(j=i+1; j<n; j++) {
1851       bExcl = FALSE;
1852       for(k=excl->index[i]; k<excl->index[i+1]; k++) {
1853         if (excl->a[k] == j)
1854           bExcl = TRUE;
1855       }
1856       if (!bExcl) {
1857         if (nb_funct != F_LJ)
1858           gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
1859         param.a[0] = i;
1860         param.a[1] = j;
1861         param.c[0] = atom[i].q;
1862         param.c[1] = atom[j].q;
1863         param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
1864         param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
1865         add_param_to_list(&mol->plist[F_LJC_PAIRS_NB],&param);
1866       }
1867     }
1868   }
1869 }
1870
1871 static void set_excl_all(t_blocka *excl)
1872 {
1873   int nat,i,j,k;
1874
1875   /* Get rid of the current exclusions and exclude all atom pairs */
1876   nat = excl->nr;
1877   excl->nra = nat*nat;
1878   srenew(excl->a,excl->nra);
1879   k = 0;
1880   for(i=0; i<nat; i++) {
1881     excl->index[i] = k;
1882     for(j=0; j<nat; j++)
1883       excl->a[k++] = j;
1884   }
1885   excl->index[nat] = k;
1886 }
1887
1888 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
1889                            int couple_lam0,int couple_lam1)
1890 {
1891   int i;
1892
1893   for(i=0; i<atoms->nr; i++) {
1894     if (couple_lam0 != ecouplamVDWQ)
1895       atoms->atom[i].q     = 0.0;
1896     if (couple_lam0 == ecouplamNONE)
1897       atoms->atom[i].type  = atomtype_decouple;
1898     if (couple_lam1 != ecouplamVDWQ)
1899       atoms->atom[i].qB    = 0.0;
1900     if (couple_lam1 == ecouplamNONE)
1901       atoms->atom[i].typeB = atomtype_decouple;
1902   }
1903 }
1904
1905 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
1906                             int couple_lam0,int couple_lam1,
1907                             bool bCoupleIntra,int nb_funct,t_params *nbp)
1908 {
1909   convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
1910   if (!bCoupleIntra) {
1911     generate_LJCpairsNB(mol,nb_funct,nbp);
1912     set_excl_all(&mol->excls);
1913   }
1914   decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
1915 }