src/kernel/toppush.c

   1 /*
   2  * $Id$
   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     }
 536     if (!bFound) {
 537         /* alloc */
 538         pr_alloc (2,bt);
 539
 540         /* fill the arrays up and down */
 541         memcpy(bt->param[bt->nr].c,  b->c,sizeof(b->c));
 542         memcpy(bt->param[bt->nr].a,  b->a,sizeof(b->a));
 543         memcpy(bt->param[bt->nr+1].c,b->c,sizeof(b->c));
 544
 545         for (j=0; (j < nral); j++)
 546             bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
 547
 548         bt->nr += 2;
 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         if ((typeB = get_atomtype_type(ctypeB,atype)) == NOTSET) {
 968           gmx_fatal(FARGS,"Atomtype %s not found",ctypeB);
 969         }
 970         qB = get_atomtype_qA(typeB,atype);
 971         mB = get_atomtype_massA(typeB,atype);
 972         if (nscan > 3) {
 973           qB = qb;
 974           if (nscan > 4) {
 975             mB = mb;
 976             if (nscan > 5)
 977               warning_error("Too many parameters");
 978           }
 979         }
 980       }
 981     }
 982   }
 983   if (debug)
 984     fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
 985
 986   push_cg(cgs,lastcg,cgnumber,nr);
 987
 988   push_atom_now(symtab,at,atomnr,get_atomtype_atomnumber(type,atype),
 989                 type,ctype,ptype,resnumber,cgnumber,
 990                 resname,name,m0,q0,typeB,
 991                 typeB==type ? ctype : ctypeB,mB,qB);
 992 }
 993
 994 void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line)
 995 {
 996   char type[STRLEN];
 997   int nrexcl,i;
 998   t_molinfo *newmol;
 999
1000   if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
1001     warning_error("Expected a molecule type name and nrexcl");
1002   }
1003
1004   /* Test if this atomtype overwrites another */
1005   i = 0;
1006   while (i < *nmol) {
1007     if (strcasecmp(*((*mol)[i].name),type) == 0)
1008       gmx_fatal(FARGS,"moleculetype %s is redefined",type);
1009     i++;
1010   }
1011
1012   (*nmol)++;
1013   srenew(*mol,*nmol);
1014   newmol = &((*mol)[*nmol-1]);
1015   init_molinfo(newmol);
1016
1017   /* Fill in the values */
1018   newmol->name     = put_symtab(symtab,type);
1019   newmol->nrexcl   = nrexcl;
1020   newmol->excl_set = FALSE;
1021 }
1022
1023 static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
1024                               t_param *p,int c_start,bool bB,bool bGenPairs)
1025 {
1026   int      i,j,ti,tj,ntype;
1027   bool     bFound;
1028   t_param  *pi=NULL;
1029   int      nr   = bt[ftype].nr;
1030   int      nral = NRAL(ftype);
1031   int      nrfp = interaction_function[ftype].nrfpA;
1032   int      nrfpB= interaction_function[ftype].nrfpB;
1033
1034   if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1035     return TRUE;
1036
1037   bFound=FALSE;
1038   if(bGenPairs) {
1039     /* First test the generated-pair position to save
1040      * time when we have 1000*1000 entries for e.g. OPLS...
1041      */
1042     ntype=sqrt(nr);
1043     if (bB) {
1044       ti=at->atom[p->a[0]].typeB;
1045       tj=at->atom[p->a[1]].typeB;
1046     } else {
1047       ti=at->atom[p->a[0]].type;
1048       tj=at->atom[p->a[1]].type;
1049     }
1050     pi=&(bt[ftype].param[ntype*ti+tj]);
1051     bFound=((ti == pi->a[0]) && (tj == pi->a[1]));
1052   }
1053
1054   /* Search explicitly if we didnt find it */
1055   if(!bFound) {
1056     for (i=0; ((i < nr) && !bFound); i++) {
1057       pi=&(bt[ftype].param[i]);
1058       if (bB)
1059         for (j=0; ((j < nral) &&
1060                    (at->atom[p->a[j]].typeB == pi->a[j])); j++);
1061       else
1062         for (j=0; ((j < nral) &&
1063                    (at->atom[p->a[j]].type == pi->a[j])); j++);
1064       bFound=(j==nral);
1065     }
1066   }
1067
1068   if (bFound) {
1069     if (bB) {
1070       if (nrfp+nrfpB > MAXFORCEPARAM) {
1071         gmx_incons("Too many force parameters");
1072       }
1073       for (j=c_start; (j < nrfpB); j++)
1074         p->c[nrfp+j] = pi->c[j];
1075     }
1076     else
1077       for (j=c_start; (j < nrfp); j++)
1078         p->c[j] = pi->c[j];
1079   }
1080   else {
1081     for (j=c_start; (j < nrfp); j++)
1082       p->c[j] = 0.0;
1083   }
1084   return bFound;
1085 }
1086
1087
1088 static bool default_params(int ftype,t_params bt[],
1089                            t_atoms *at,t_atomtype atype,
1090                            t_param *p,bool bB,
1091                            t_param **param_def,
1092                int *nparam_def)
1093 {
1094   int      i,j,nparam_found;
1095   bool     bFound,bSame;
1096   t_param  *pi=NULL;
1097   t_param  *pj=NULL;
1098   int      nr   = bt[ftype].nr;
1099   int      nral = NRAL(ftype);
1100   int      nrfpA= interaction_function[ftype].nrfpA;
1101   int      nrfpB= interaction_function[ftype].nrfpB;
1102
1103   if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1104     return TRUE;
1105
1106
1107   /* We allow wildcards now. The first type (with or without wildcards) that
1108    * fits is used, so you should probably put the wildcarded bondtypes
1109    * at the end of each section.
1110    */
1111   bFound=FALSE;
1112   nparam_found=0;
1113   /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1114    * special case for this. Check for B state outside loop to speed it up.
1115    */
1116   if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS)
1117   {
1118           if (bB)
1119           {
1120                   for (i=0; ((i < nr) && !bFound); i++)
1121                   {
1122                           pi=&(bt[ftype].param[i]);
1123                           bFound=
1124                           (
1125                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1126                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1127                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1128                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1129                            );
1130                   }
1131           }
1132           else
1133           {
1134                   /* State A */
1135                   for (i=0; ((i < nr) && !bFound); i++)
1136                   {
1137                           pi=&(bt[ftype].param[i]);
1138                           bFound=
1139                           (
1140                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].type,atype)==pi->AI)) &&
1141                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].type,atype)==pi->AJ)) &&
1142                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].type,atype)==pi->AK)) &&
1143                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].type,atype)==pi->AL))
1144                            );
1145                   }
1146           }
1147           /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1148            * The rules in that case is that additional matches HAVE to be on adjacent lines!
1149            */
1150           if(bFound==TRUE)
1151           {
1152                   nparam_found++;
1153                   bSame = TRUE;
1154                   /* Continue from current i value */
1155                   for(j=i+1 ; j<nr && bSame ; j+=2)
1156                   {
1157                           pj=&(bt[ftype].param[j]);
1158                           bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1159                           if(bSame)
1160                                   nparam_found++;
1161                           /* nparam_found will be increased as long as the numbers match */
1162                   }
1163           }
1164   } else { /* Not a dihedral */
1165           for (i=0; ((i < nr) && !bFound); i++) {
1166                   pi=&(bt[ftype].param[i]);
1167                   if (bB)
1168                           for (j=0; ((j < nral) &&
1169                                                  (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
1170                                   ;
1171                   else
1172                           for (j=0; ((j < nral) &&
1173                                                  (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
1174                                   ;
1175                   bFound=(j==nral);
1176           }
1177           if(bFound)
1178                   nparam_found=1;
1179   }
1180
1181   *param_def = pi;
1182   *nparam_def = nparam_found;
1183
1184   return bFound;
1185 }
1186
1187
1188
1189 void push_bondnow(t_params *bond, t_param *b)
1190 {
1191   int j;
1192
1193   if (debug) {
1194     fprintf(debug,"push_bondnow: nr = %d\n",bond->nr);
1195     fflush(debug);
1196   }
1197   /* allocate one position extra */
1198   pr_alloc (1,bond);
1199
1200   /* fill the arrays */
1201   for (j=0; (j < MAXFORCEPARAM); j++)
1202     bond->param[bond->nr].c[j]   = b->c[j];
1203   for (j=0; (j < MAXATOMLIST); j++)
1204     bond->param[bond->nr].a[j]   = b->a[j];
1205
1206   bond->nr++;
1207 }
1208
1209 void push_bond(directive d,t_params bondtype[],t_params bond[],
1210                t_atoms *at,t_atomtype atype,char *line,
1211                bool bBonded,bool bGenPairs,real fudgeQQ,
1212                bool bZero,bool *bWarn_copy_A_B)
1213 {
1214   const char *aaformat[MAXATOMLIST]= {
1215     "%d%d",
1216     "%d%d%d",
1217     "%d%d%d%d",
1218     "%d%d%d%d%d",
1219     "%d%d%d%d%d%d"
1220   };
1221   const char *asformat[MAXATOMLIST]= {
1222     "%*s%*s",
1223     "%*s%*s%*s",
1224     "%*s%*s%*s%*s",
1225     "%*s%*s%*s%*s%*s",
1226     "%*s%*s%*s%*s%*s%*s"
1227   };
1228   const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1229   int      nr,i,j,nral,nread,ftype;
1230   char     format[STRLEN];
1231   /* One force parameter more, so we can check if we read too many */
1232   double   cc[MAXFORCEPARAM+1];
1233   int      aa[MAXATOMLIST+1];
1234   t_param  param,paramB,*param_defA,*param_defB;
1235   bool     bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
1236   int      nparam_defA,nparam_defB;
1237   char  errbuf[256];
1238
1239   nparam_defA=nparam_defB=0;
1240
1241   ftype = ifunc_index(d,1);
1242   nral  = NRAL(ftype);
1243   for(j=0; j<MAXATOMLIST; j++)
1244     aa[j]=NOTSET;
1245   bDef = (NRFP(ftype) > 0);
1246
1247   nread = sscanf(line,aaformat[nral-1],
1248                  &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1249   if (nread < nral) {
1250     too_few();
1251     return;
1252   } else if (nread == nral)
1253     ftype = ifunc_index(d,1);
1254   else {
1255     /* this is a hack to allow for virtual sites with swapped parity */
1256     bSwapParity = (aa[nral]<0);
1257     if (bSwapParity)
1258       aa[nral] = -aa[nral];
1259     ftype = ifunc_index(d,aa[nral]);
1260     if (bSwapParity)
1261       switch(ftype) {
1262       case F_VSITE3FAD:
1263       case F_VSITE3OUT:
1264         break;
1265       default:
1266         gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1267                     interaction_function[F_VSITE3FAD].longname,
1268                     interaction_function[F_VSITE3OUT].longname);
1269       }
1270   }
1271
1272   /* Check for double atoms and atoms out of bounds */
1273   for(i=0; (i<nral); i++) {
1274     if ( aa[i] < 1 || aa[i] > at->nr )
1275       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1276                 "Atom index (%d) in %s out of bounds (1-%d).\n"
1277                 "This probably means that you have inserted topology section \"%s\"\n"
1278                 "in a part belonging to a different molecule than you intended to.\n"
1279                 "In that case move the \"%s\" section to the right molecule.",
1280                 get_warning_file(),get_warning_line(),
1281                 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1282     for(j=i+1; (j<nral); j++)
1283       if (aa[i] == aa[j]) {
1284         sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1285         warning(errbuf);
1286       }
1287   }
1288   if (ftype == F_SETTLE)
1289     if (aa[0]+2 > at->nr)
1290       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1291                   "             Atom index (%d) in %s out of bounds (1-%d)\n"
1292                   "             Settle works on atoms %d, %d and %d",
1293                   get_warning_file(),get_warning_line(),
1294                   aa[0],dir2str(d),at->nr,
1295                   aa[0],aa[0]+1,aa[0]+2);
1296
1297   /* default force parameters  */
1298   for(j=0; (j<MAXATOMLIST); j++)
1299     param.a[j] = aa[j]-1;
1300   for(j=0; (j<MAXFORCEPARAM); j++)
1301     param.c[j] = 0.0;
1302
1303   /* Get force params for normal and free energy perturbation
1304    * studies, as determined by types!
1305    */
1306   if (bBonded) {
1307     bFoundA = default_params(ftype,bondtype,at,atype,&param,FALSE,&param_defA,&nparam_defA);
1308     if (bFoundA) {
1309       /* Copy the A-state and B-state default parameters */
1310       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1311         param.c[j] = param_defA->c[j];
1312     }
1313     bFoundB = default_params(ftype,bondtype,at,atype,&param,TRUE,&param_defB,&nparam_defB);
1314     if (bFoundB) {
1315       /* Copy only the B-state default parameters */
1316       for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1317         param.c[j] = param_defB->c[j];
1318     }
1319   } else if (ftype == F_LJ14) {
1320     bFoundA = default_nb_params(ftype, bondtype,at,&param,0,FALSE,bGenPairs);
1321     bFoundB = default_nb_params(ftype, bondtype,at,&param,0,TRUE, bGenPairs);
1322   } else if (ftype == F_LJC14_Q) {
1323     param.c[0] = fudgeQQ;
1324     /* Fill in the A-state charges as default parameters */
1325     param.c[1] = at->atom[param.a[0]].q;
1326     param.c[2] = at->atom[param.a[1]].q;
1327     /* The default LJ parameters are the standard 1-4 parameters */
1328     bFoundA = default_nb_params(F_LJ14,bondtype,at,&param,3,FALSE,bGenPairs);
1329     bFoundB = TRUE;
1330   } else if (ftype == F_LJC_PAIRS_NB) {
1331     /* Defaults are not supported here */
1332     bFoundA = FALSE;
1333     bFoundB = TRUE;
1334   } else {
1335     gmx_incons("Unknown function type in push_bond");
1336   }
1337
1338   if (nread > nral) {
1339           /* Manually specified parameters - in this case we discard multiple torsion info! */
1340
1341     strcpy(format,asformat[nral-1]);
1342     strcat(format,ccformat);
1343
1344     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1345                    &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1346
1347     if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1348       {
1349         /* We only have to issue a warning if these atoms are perturbed! */
1350         bPert = FALSE;
1351         for(j=0; (j<nral); j++)
1352           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1353
1354         if (bPert && *bWarn_copy_A_B)
1355           {
1356             sprintf(errbuf,
1357                     "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1358             warning(errbuf);
1359             *bWarn_copy_A_B = FALSE;
1360           }
1361
1362         /* If only the A parameters were specified, copy them to the B state */
1363         /* The B-state parameters correspond to the first nrfpB
1364          * A-state parameters.
1365          */
1366         for(j=0; (j<NRFPB(ftype)); j++)
1367           cc[nread++] = cc[j];
1368       }
1369
1370     /* If nread was 0 or EOF, no parameters were read => use defaults.
1371      * If nread was nrfpA we copied above so nread=nrfp.
1372      * If nread was nrfp we are cool.
1373      * For F_LJC14_Q we allow supplying fudgeQQ only.
1374      * Anything else is an error!
1375      */
1376     if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1377         !(ftype == F_LJC14_Q && nread == 1))
1378       {
1379         gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1380                   nread,NRFPA(ftype),NRFP(ftype),
1381                   interaction_function[ftype].longname);
1382       }
1383
1384     for(j=0; (j<nread); j++)
1385       param.c[j]=cc[j];
1386
1387     /* Check whether we have to use the defaults */
1388     if (nread == NRFP(ftype))
1389       bDef = FALSE;
1390   }
1391   else
1392   {
1393     nread = 0;
1394   }
1395   /* nread now holds the number of force parameters read! */
1396
1397
1398
1399
1400         if (bDef)
1401         {
1402           /* Use defaults */
1403                 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1404                 if(ftype==F_PDIHS)
1405                 {
1406                         if ((nparam_defA != nparam_defB) || ((nparam_defA>1 || nparam_defB>1) && (param_defA!=param_defB)))
1407                         {
1408                                 sprintf(errbuf,
1409                                                 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1410                                                 "Please specify perturbed parameters manually for this torsion in your topology!");
1411                                 warning_error(errbuf);
1412                         }
1413                 }
1414
1415           if (nread > 0 && nread < NRFPA(ftype))
1416           {
1417                   /* Issue an error, do not use defaults */
1418                   sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1419                   warning_error(errbuf);
1420           }
1421
1422       if (nread == 0 || nread == EOF)
1423           {
1424                   if (!bFoundA)
1425                   {
1426                           if (interaction_function[ftype].flags & IF_VSITE)
1427                           {
1428                                   /* set them to NOTSET, will be calculated later */
1429                                   for(j=0; (j<MAXFORCEPARAM); j++)
1430                                           param.c[j] = NOTSET;
1431
1432                                   if (bSwapParity)
1433                                           param.C1 = -1; /* flag to swap parity of vsite construction */
1434                           }
1435                           else
1436                           {
1437                                   if (bZero)
1438                                   {
1439                                           fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1440                                                           interaction_function[ftype].longname);
1441                                   }
1442                                   else
1443                                   {
1444                                           sprintf(errbuf,"No default %s types",interaction_function[ftype].longname);
1445                                           warning_error(errbuf);
1446                                   }
1447                           }
1448                   }
1449                   else
1450                   {
1451                           if (bSwapParity)
1452                           {
1453                                   switch(ftype)
1454                                   {
1455                                           case F_VSITE3FAD:
1456                                                   param.C0 = 360 - param.C0;
1457                                                   break;
1458                                           case F_VSITE3OUT:
1459                                                   param.C2 = -param.C2;
1460                                                   break;
1461                                   }
1462                           }
1463                   }
1464                   if (!bFoundB)
1465                   {
1466                           /* We only have to issue a warning if these atoms are perturbed! */
1467                           bPert = FALSE;
1468                           for(j=0; (j<nral); j++)
1469                                   bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1470
1471                           if (bPert)
1472                           {
1473                                   sprintf(errbuf,"No default %s types for perturbed atoms, "
1474                                                   "using normal values",interaction_function[ftype].longname);
1475                                   warning(errbuf);
1476                           }
1477                   }
1478           }
1479   }
1480
1481         if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1482                 && param.c[5]!=param.c[2])
1483                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1484                                   "             %s multiplicity can not be perturbed %f!=%f",
1485                                   get_warning_file(),get_warning_line(),
1486                                   interaction_function[ftype].longname,
1487                                   param.c[2],param.c[5]);
1488
1489         if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1490                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1491                                   "             %s table number can not be perturbed %d!=%d",
1492                                   get_warning_file(),get_warning_line(),
1493                                   interaction_function[ftype].longname,
1494                                   (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1495
1496         /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1497         if (ftype==F_RBDIHS) {
1498                 nr=0;
1499                 for(i=0;i<NRFP(ftype);i++) {
1500                         if(param.c[i]!=0)
1501                                 nr++;
1502                 }
1503                 if(nr==0)
1504                         return;
1505         }
1506
1507         /* Put the values in the appropriate arrays */
1508         push_bondnow (&bond[ftype],&param);
1509
1510         /* Push additional torsions from FF for ftype==9 if we have them.
1511          * We have already checked that the A/B states do not differ in this case,
1512          * so we do not have to double-check that again, or the vsite stuff.
1513          * In addition, those torsions cannot be automatically perturbed.
1514          */
1515         if(bDef && ftype==F_PDIHS)
1516         {
1517                 for(i=1;i<nparam_defA;i++)
1518                 {
1519                         /* Advance pointer! */
1520                         param_defA+=2;
1521                         for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1522                                 param.c[j] = param_defA->c[j];
1523                         /* And push the next term for this torsion */
1524                         push_bondnow (&bond[ftype],&param);
1525                 }
1526         }
1527 }
1528
1529
1530
1531 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1532                   t_atoms *at,t_atomtype atype,char *line)
1533 {
1534   char *ptr;
1535   int  type,ftype,j,n,ret,nj,a;
1536   int  *atc=NULL;
1537   double *weight=NULL,weight_tot;
1538   t_param param;
1539
1540   /* default force parameters  */
1541   for(j=0; (j<MAXATOMLIST); j++)
1542     param.a[j] = NOTSET;
1543   for(j=0; (j<MAXFORCEPARAM); j++)
1544     param.c[j] = 0.0;
1545
1546   ptr = line;
1547   ret = sscanf(ptr,"%d%n",&a,&n);
1548   ptr += n;
1549   if (ret == 0)
1550     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1551               "             Expected an atom index in section \"%s\"",
1552                 get_warning_file(),get_warning_line(),
1553               dir2str(d));
1554
1555   param.a[0] = a - 1;
1556
1557   ret = sscanf(ptr,"%d%n",&type,&n);
1558   ptr += n;
1559   ftype = ifunc_index(d,type);
1560
1561   weight_tot = 0;
1562   nj = 0;
1563   do {
1564     ret = sscanf(ptr,"%d%n",&a,&n);
1565     ptr += n;
1566     if (ret > 0) {
1567       if (nj % 20 == 0) {
1568         srenew(atc,nj+20);
1569         srenew(weight,nj+20);
1570       }
1571       atc[nj] = a - 1;
1572       switch (type) {
1573       case 1:
1574         weight[nj] = 1;
1575         break;
1576       case 2:
1577         /* Here we use the A-state mass as a parameter.
1578          * Note that the B-state mass has no influence.
1579          */
1580         weight[nj] = at->atom[atc[nj]].m;
1581         break;
1582       case 3:
1583         weight[nj] = -1;
1584         ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1585         ptr += n;
1586         if (weight[nj] < 0)
1587           gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1588                     "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1589                     get_warning_file(),get_warning_line(),
1590                     nj+1,atc[nj]+1);
1591         break;
1592       default:
1593         gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1594       }
1595       weight_tot += weight[nj];
1596       nj++;
1597     }
1598   } while (ret > 0);
1599
1600   if (nj == 0)
1601     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1602               "             Expected more than one atom index in section \"%s\"",
1603               get_warning_file(),get_warning_line(),
1604               dir2str(d));
1605
1606   if (weight_tot == 0)
1607      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1608                "             The total mass of the construting atoms is zero",
1609                get_warning_file(),get_warning_line());
1610
1611   for(j=0; j<nj; j++) {
1612     param.a[1] = atc[j];
1613     param.c[0] = nj;
1614     param.c[1] = weight[j]/weight_tot;
1615     /* Put the values in the appropriate arrays */
1616     push_bondnow (&bond[ftype],&param);
1617   }
1618
1619   sfree(atc);
1620   sfree(weight);
1621 }
1622
1623 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1624               int *nrcopies)
1625 {
1626   int  i,copies;
1627   char type[STRLEN];
1628
1629   *nrcopies=0;
1630   if (sscanf(pline,"%s%d",type,&copies) != 2) {
1631     too_few();
1632     return;
1633   }
1634
1635   /* search moleculename */
1636   for (i=0; ((i<nrmols) && strcasecmp(type,*(mols[i].name))); i++)
1637     ;
1638
1639   if (i<nrmols) {
1640     *nrcopies        = copies;
1641     *whichmol        = i;
1642   } else
1643     gmx_fatal(FARGS,"No such moleculetype %s",type);
1644 }
1645
1646 void init_block2(t_block2 *b2, int natom)
1647 {
1648   int i;
1649
1650   b2->nr=natom;
1651   snew(b2->nra,b2->nr);
1652   snew(b2->a,b2->nr);
1653   for(i=0; (i<b2->nr); i++)
1654     b2->a[i]=NULL;
1655 }
1656
1657 void done_block2(t_block2 *b2)
1658 {
1659   int i;
1660
1661   if (b2->nr) {
1662     for(i=0; (i<b2->nr); i++)
1663       sfree(b2->a[i]);
1664     sfree(b2->a);
1665     sfree(b2->nra);
1666     b2->nr=0;
1667   }
1668 }
1669
1670 void push_excl(char *line, t_block2 *b2)
1671 {
1672   int  i,j;
1673   int  n;
1674   char base[STRLEN],format[STRLEN];
1675
1676   if (sscanf(line,"%d",&i)==0)
1677     return;
1678
1679   if ((1 <= i) && (i <= b2->nr))
1680     i--;
1681   else {
1682     if (debug)
1683       fprintf(debug,"Unbound atom %d\n",i-1);
1684     return;
1685   }
1686   strcpy(base,"%*d");
1687   do {
1688     strcpy(format,base);
1689     strcat(format,"%d");
1690     n=sscanf(line,format,&j);
1691     if (n == 1) {
1692       if ((1 <= j) && (j <= b2->nr)) {
1693         j--;
1694         srenew(b2->a[i],++(b2->nra[i]));
1695         b2->a[i][b2->nra[i]-1]=j;
1696         /* also add the reverse exclusion! */
1697         srenew(b2->a[j],++(b2->nra[j]));
1698         b2->a[j][b2->nra[j]-1]=i;
1699         strcat(base,"%*d");
1700       }
1701       else
1702         gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
1703     }
1704   } while (n == 1);
1705 }
1706
1707 void b_to_b2(t_blocka *b, t_block2 *b2)
1708 {
1709   int     i;
1710   atom_id j,a;
1711
1712   for(i=0; (i<b->nr); i++)
1713     for(j=b->index[i]; (j<b->index[i+1]); j++) {
1714       a=b->a[j];
1715       srenew(b2->a[i],++b2->nra[i]);
1716       b2->a[i][b2->nra[i]-1]=a;
1717     }
1718 }
1719
1720 void b2_to_b(t_block2 *b2, t_blocka *b)
1721 {
1722   int     i,nra;
1723   atom_id j;
1724
1725   nra=0;
1726   for(i=0; (i<b2->nr); i++) {
1727     b->index[i]=nra;
1728     for(j=0; (j<b2->nra[i]); j++)
1729       b->a[nra+j]=b2->a[i][j];
1730     nra+=b2->nra[i];
1731   }
1732   /* terminate list */
1733   b->index[i]=nra;
1734 }
1735
1736 static int icomp(const void *v1, const void *v2)
1737 {
1738   return (*((atom_id *) v1))-(*((atom_id *) v2));
1739 }
1740
1741 void merge_excl(t_blocka *excl, t_block2 *b2)
1742 {
1743   int     i,k;
1744   atom_id j;
1745   int     nra;
1746
1747   if (!b2->nr)
1748     return;
1749   else if (b2->nr != excl->nr) {
1750     gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
1751                 b2->nr,excl->nr);
1752   }
1753   else if (debug)
1754     fprintf(debug,"Entering merge_excl\n");
1755
1756   /* First copy all entries from excl to b2 */
1757   b_to_b2(excl,b2);
1758
1759   /* Count and sort the exclusions */
1760   nra=0;
1761   for(i=0; (i<b2->nr); i++) {
1762     if (b2->nra[i] > 0) {
1763       /* remove double entries */
1764       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
1765       k=1;
1766       for(j=1; (j<b2->nra[i]); j++)
1767         if (b2->a[i][j]!=b2->a[i][k-1]) {
1768           b2->a[i][k]=b2->a[i][j];
1769           k++;
1770         }
1771       b2->nra[i]=k;
1772       nra+=b2->nra[i];
1773     }
1774   }
1775   excl->nra=nra;
1776   srenew(excl->a,excl->nra);
1777
1778   b2_to_b(b2,excl);
1779 }
1780
1781 int add_atomtype_decoupled(t_symtab *symtab,t_atomtype at,
1782                            t_nbparam ***nbparam,t_nbparam ***pair)
1783 {
1784   t_atom  atom;
1785   t_param param;
1786   int     i,nr;
1787
1788   /* Define an atom type with all parameters set to zero (no interactions) */
1789   atom.q = 0.0;
1790   atom.m = 0.0;
1791   /* Type for decoupled atoms could be anything,
1792    * this should be changed automatically later when required.
1793    */
1794   atom.ptype = eptAtom;
1795   for (i=0; (i<MAXFORCEPARAM); i++)
1796     param.c[i] = 0.0;
1797
1798   nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0);
1799
1800   /* Add space in the non-bonded parameters matrix */
1801   realloc_nb_params(at,nbparam,pair);
1802
1803   return nr;
1804 }
1805
1806 static void convert_pairs_to_pairsQ(t_params *plist,
1807                                     real fudgeQQ,t_atoms *atoms)
1808 {
1809   t_param *param;
1810   int  i;
1811   real v,w;
1812
1813   /* Copy the pair list to the pairQ list */
1814   plist[F_LJC14_Q] = plist[F_LJ14];
1815   /* Empty the pair list */
1816   plist[F_LJ14].nr    = 0;
1817   plist[F_LJ14].param = NULL;
1818   param = plist[F_LJC14_Q].param;
1819   for(i=0; i<plist[F_LJC14_Q].nr; i++) {
1820     v = param[i].c[0];
1821     w = param[i].c[1];
1822     param[i].c[0] = fudgeQQ;
1823     param[i].c[1] = atoms->atom[param[i].a[0]].q;
1824     param[i].c[2] = atoms->atom[param[i].a[1]].q;
1825     param[i].c[3] = v;
1826     param[i].c[4] = w;
1827   }
1828 }
1829
1830 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
1831 {
1832   int n,ntype,i,j,k;
1833   t_atom *atom;
1834   t_blocka *excl;
1835   bool bExcl;
1836   t_param param;
1837
1838   n = mol->atoms.nr;
1839   atom = mol->atoms.atom;
1840
1841   ntype = sqrt(nbp->nr);
1842
1843   for (i=0; i<MAXATOMLIST; i++)
1844     param.a[i] = NOTSET;
1845   for (i=0; i<MAXFORCEPARAM; i++)
1846     param.c[i] = NOTSET;
1847
1848   /* Add a pair interaction for all non-excluded atom pairs */
1849   excl = &mol->excls;
1850   for(i=0; i<n; i++) {
1851     for(j=i+1; j<n; j++) {
1852       bExcl = FALSE;
1853       for(k=excl->index[i]; k<excl->index[i+1]; k++) {
1854         if (excl->a[k] == j)
1855           bExcl = TRUE;
1856       }
1857       if (!bExcl) {
1858         if (nb_funct != F_LJ)
1859           gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
1860         param.a[0] = i;
1861         param.a[1] = j;
1862         param.c[0] = atom[i].q;
1863         param.c[1] = atom[j].q;
1864         param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
1865         param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
1866         push_bondnow(&mol->plist[F_LJC_PAIRS_NB],&param);
1867       }
1868     }
1869   }
1870 }
1871
1872 static void set_excl_all(t_blocka *excl)
1873 {
1874   int nat,i,j,k;
1875
1876   /* Get rid of the current exclusions and exclude all atom pairs */
1877   nat = excl->nr;
1878   excl->nra = nat*nat;
1879   srenew(excl->a,excl->nra);
1880   k = 0;
1881   for(i=0; i<nat; i++) {
1882     excl->index[i] = k;
1883     for(j=0; j<nat; j++)
1884       excl->a[k++] = j;
1885   }
1886   excl->index[nat] = k;
1887 }
1888
1889 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
1890                            int couple_lam0,int couple_lam1)
1891 {
1892   int i;
1893
1894   for(i=0; i<atoms->nr; i++) {
1895     if (couple_lam0 != ecouplamVDWQ)
1896       atoms->atom[i].q     = 0.0;
1897     if (couple_lam0 == ecouplamNONE)
1898       atoms->atom[i].type  = atomtype_decouple;
1899     if (couple_lam1 != ecouplamVDWQ)
1900       atoms->atom[i].qB    = 0.0;
1901     if (couple_lam1 == ecouplamNONE)
1902       atoms->atom[i].typeB = atomtype_decouple;
1903   }
1904 }
1905
1906 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
1907                             int couple_lam0,int couple_lam1,
1908                             bool bCoupleIntra,int nb_funct,t_params *nbp)
1909 {
1910   convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
1911   if (!bCoupleIntra) {
1912     generate_LJCpairsNB(mol,nb_funct,nbp);
1913     set_excl_all(&mol->excls);
1914   }
1915   decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
1916 }