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         add_param_to_list (&bond[ftype],&param);
1487
1488   /* Push additional torsions from FF for ftype==9 if we have them.
1489    * We have already checked that the A/B states do not differ in this case,
1490    * so we do not have to double-check that again, or the vsite stuff.
1491    * In addition, those torsions cannot be automatically perturbed.
1492    */
1493   if(bDef && ftype==F_PDIHS)
1494   {
1495     for(i=1;i<nparam_defA;i++)
1496     {
1497       /* Advance pointer! */
1498       param_defA++;
1499       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1500         param.c[j] = param_defA->c[j];
1501       /* And push the next term for this torsion */
1502       add_param_to_list (&bond[ftype],&param);
1503     }
1504   }
1505 }
1506
1507
1508
1509 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1510                   t_atoms *at,t_atomtype atype,char *line)
1511 {
1512   char *ptr;
1513   int  type,ftype,j,n,ret,nj,a;
1514   int  *atc=NULL;
1515   double *weight=NULL,weight_tot;
1516   t_param param;
1517
1518   /* default force parameters  */
1519   for(j=0; (j<MAXATOMLIST); j++)
1520     param.a[j] = NOTSET;
1521   for(j=0; (j<MAXFORCEPARAM); j++)
1522     param.c[j] = 0.0;
1523
1524   ptr = line;
1525   ret = sscanf(ptr,"%d%n",&a,&n);
1526   ptr += n;
1527   if (ret == 0)
1528     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1529               "             Expected an atom index in section \"%s\"",
1530                 get_warning_file(),get_warning_line(),
1531               dir2str(d));
1532
1533   param.a[0] = a - 1;
1534
1535   ret = sscanf(ptr,"%d%n",&type,&n);
1536   ptr += n;
1537   ftype = ifunc_index(d,type);
1538
1539   weight_tot = 0;
1540   nj = 0;
1541   do {
1542     ret = sscanf(ptr,"%d%n",&a,&n);
1543     ptr += n;
1544     if (ret > 0) {
1545       if (nj % 20 == 0) {
1546         srenew(atc,nj+20);
1547         srenew(weight,nj+20);
1548       }
1549       atc[nj] = a - 1;
1550       switch (type) {
1551       case 1:
1552         weight[nj] = 1;
1553         break;
1554       case 2:
1555         /* Here we use the A-state mass as a parameter.
1556          * Note that the B-state mass has no influence.
1557          */
1558         weight[nj] = at->atom[atc[nj]].m;
1559         break;
1560       case 3:
1561         weight[nj] = -1;
1562         ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1563         ptr += n;
1564         if (weight[nj] < 0)
1565           gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1566                     "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1567                     get_warning_file(),get_warning_line(),
1568                     nj+1,atc[nj]+1);
1569         break;
1570       default:
1571         gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1572       }
1573       weight_tot += weight[nj];
1574       nj++;
1575     }
1576   } while (ret > 0);
1577
1578   if (nj == 0)
1579     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1580               "             Expected more than one atom index in section \"%s\"",
1581               get_warning_file(),get_warning_line(),
1582               dir2str(d));
1583
1584   if (weight_tot == 0)
1585      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1586                "             The total mass of the construting atoms is zero",
1587                get_warning_file(),get_warning_line());
1588
1589   for(j=0; j<nj; j++) {
1590     param.a[1] = atc[j];
1591     param.c[0] = nj;
1592     param.c[1] = weight[j]/weight_tot;
1593     /* Put the values in the appropriate arrays */
1594     add_param_to_list (&bond[ftype],&param);
1595   }
1596
1597   sfree(atc);
1598   sfree(weight);
1599 }
1600
1601 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1602               int *nrcopies)
1603 {
1604   int  i,copies;
1605   char type[STRLEN];
1606
1607   *nrcopies=0;
1608   if (sscanf(pline,"%s%d",type,&copies) != 2) {
1609     too_few();
1610     return;
1611   }
1612
1613   /* search moleculename */
1614   for (i=0; ((i<nrmols) && strcasecmp(type,*(mols[i].name))); i++)
1615     ;
1616
1617   if (i<nrmols) {
1618     *nrcopies        = copies;
1619     *whichmol        = i;
1620   } else
1621     gmx_fatal(FARGS,"No such moleculetype %s",type);
1622 }
1623
1624 void init_block2(t_block2 *b2, int natom)
1625 {
1626   int i;
1627
1628   b2->nr=natom;
1629   snew(b2->nra,b2->nr);
1630   snew(b2->a,b2->nr);
1631   for(i=0; (i<b2->nr); i++)
1632     b2->a[i]=NULL;
1633 }
1634
1635 void done_block2(t_block2 *b2)
1636 {
1637   int i;
1638
1639   if (b2->nr) {
1640     for(i=0; (i<b2->nr); i++)
1641       sfree(b2->a[i]);
1642     sfree(b2->a);
1643     sfree(b2->nra);
1644     b2->nr=0;
1645   }
1646 }
1647
1648 void push_excl(char *line, t_block2 *b2)
1649 {
1650   int  i,j;
1651   int  n;
1652   char base[STRLEN],format[STRLEN];
1653
1654   if (sscanf(line,"%d",&i)==0)
1655     return;
1656
1657   if ((1 <= i) && (i <= b2->nr))
1658     i--;
1659   else {
1660     if (debug)
1661       fprintf(debug,"Unbound atom %d\n",i-1);
1662     return;
1663   }
1664   strcpy(base,"%*d");
1665   do {
1666     strcpy(format,base);
1667     strcat(format,"%d");
1668     n=sscanf(line,format,&j);
1669     if (n == 1) {
1670       if ((1 <= j) && (j <= b2->nr)) {
1671         j--;
1672         srenew(b2->a[i],++(b2->nra[i]));
1673         b2->a[i][b2->nra[i]-1]=j;
1674         /* also add the reverse exclusion! */
1675         srenew(b2->a[j],++(b2->nra[j]));
1676         b2->a[j][b2->nra[j]-1]=i;
1677         strcat(base,"%*d");
1678       }
1679       else
1680         gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
1681     }
1682   } while (n == 1);
1683 }
1684
1685 void b_to_b2(t_blocka *b, t_block2 *b2)
1686 {
1687   int     i;
1688   atom_id j,a;
1689
1690   for(i=0; (i<b->nr); i++)
1691     for(j=b->index[i]; (j<b->index[i+1]); j++) {
1692       a=b->a[j];
1693       srenew(b2->a[i],++b2->nra[i]);
1694       b2->a[i][b2->nra[i]-1]=a;
1695     }
1696 }
1697
1698 void b2_to_b(t_block2 *b2, t_blocka *b)
1699 {
1700   int     i,nra;
1701   atom_id j;
1702
1703   nra=0;
1704   for(i=0; (i<b2->nr); i++) {
1705     b->index[i]=nra;
1706     for(j=0; (j<b2->nra[i]); j++)
1707       b->a[nra+j]=b2->a[i][j];
1708     nra+=b2->nra[i];
1709   }
1710   /* terminate list */
1711   b->index[i]=nra;
1712 }
1713
1714 static int icomp(const void *v1, const void *v2)
1715 {
1716   return (*((atom_id *) v1))-(*((atom_id *) v2));
1717 }
1718
1719 void merge_excl(t_blocka *excl, t_block2 *b2)
1720 {
1721   int     i,k;
1722   atom_id j;
1723   int     nra;
1724
1725   if (!b2->nr)
1726     return;
1727   else if (b2->nr != excl->nr) {
1728     gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
1729                 b2->nr,excl->nr);
1730   }
1731   else if (debug)
1732     fprintf(debug,"Entering merge_excl\n");
1733
1734   /* First copy all entries from excl to b2 */
1735   b_to_b2(excl,b2);
1736
1737   /* Count and sort the exclusions */
1738   nra=0;
1739   for(i=0; (i<b2->nr); i++) {
1740     if (b2->nra[i] > 0) {
1741       /* remove double entries */
1742       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
1743       k=1;
1744       for(j=1; (j<b2->nra[i]); j++)
1745         if (b2->a[i][j]!=b2->a[i][k-1]) {
1746           b2->a[i][k]=b2->a[i][j];
1747           k++;
1748         }
1749       b2->nra[i]=k;
1750       nra+=b2->nra[i];
1751     }
1752   }
1753   excl->nra=nra;
1754   srenew(excl->a,excl->nra);
1755
1756   b2_to_b(b2,excl);
1757 }
1758
1759 int add_atomtype_decoupled(t_symtab *symtab,t_atomtype at,
1760                            t_nbparam ***nbparam,t_nbparam ***pair)
1761 {
1762   t_atom  atom;
1763   t_param param;
1764   int     i,nr;
1765
1766   /* Define an atom type with all parameters set to zero (no interactions) */
1767   atom.q = 0.0;
1768   atom.m = 0.0;
1769   /* Type for decoupled atoms could be anything,
1770    * this should be changed automatically later when required.
1771    */
1772   atom.ptype = eptAtom;
1773   for (i=0; (i<MAXFORCEPARAM); i++)
1774     param.c[i] = 0.0;
1775
1776   nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0);
1777
1778   /* Add space in the non-bonded parameters matrix */
1779   realloc_nb_params(at,nbparam,pair);
1780
1781   return nr;
1782 }
1783
1784 static void convert_pairs_to_pairsQ(t_params *plist,
1785                                     real fudgeQQ,t_atoms *atoms)
1786 {
1787   t_param *param;
1788   int  i;
1789   real v,w;
1790
1791   /* Copy the pair list to the pairQ list */
1792   plist[F_LJC14_Q] = plist[F_LJ14];
1793   /* Empty the pair list */
1794   plist[F_LJ14].nr    = 0;
1795   plist[F_LJ14].param = NULL;
1796   param = plist[F_LJC14_Q].param;
1797   for(i=0; i<plist[F_LJC14_Q].nr; i++) {
1798     v = param[i].c[0];
1799     w = param[i].c[1];
1800     param[i].c[0] = fudgeQQ;
1801     param[i].c[1] = atoms->atom[param[i].a[0]].q;
1802     param[i].c[2] = atoms->atom[param[i].a[1]].q;
1803     param[i].c[3] = v;
1804     param[i].c[4] = w;
1805   }
1806 }
1807
1808 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
1809 {
1810   int n,ntype,i,j,k;
1811   t_atom *atom;
1812   t_blocka *excl;
1813   bool bExcl;
1814   t_param param;
1815
1816   n = mol->atoms.nr;
1817   atom = mol->atoms.atom;
1818
1819   ntype = sqrt(nbp->nr);
1820
1821   for (i=0; i<MAXATOMLIST; i++)
1822     param.a[i] = NOTSET;
1823   for (i=0; i<MAXFORCEPARAM; i++)
1824     param.c[i] = NOTSET;
1825
1826   /* Add a pair interaction for all non-excluded atom pairs */
1827   excl = &mol->excls;
1828   for(i=0; i<n; i++) {
1829     for(j=i+1; j<n; j++) {
1830       bExcl = FALSE;
1831       for(k=excl->index[i]; k<excl->index[i+1]; k++) {
1832         if (excl->a[k] == j)
1833           bExcl = TRUE;
1834       }
1835       if (!bExcl) {
1836         if (nb_funct != F_LJ)
1837           gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
1838         param.a[0] = i;
1839         param.a[1] = j;
1840         param.c[0] = atom[i].q;
1841         param.c[1] = atom[j].q;
1842         param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
1843         param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
1844         add_param_to_list(&mol->plist[F_LJC_PAIRS_NB],&param);
1845       }
1846     }
1847   }
1848 }
1849
1850 static void set_excl_all(t_blocka *excl)
1851 {
1852   int nat,i,j,k;
1853
1854   /* Get rid of the current exclusions and exclude all atom pairs */
1855   nat = excl->nr;
1856   excl->nra = nat*nat;
1857   srenew(excl->a,excl->nra);
1858   k = 0;
1859   for(i=0; i<nat; i++) {
1860     excl->index[i] = k;
1861     for(j=0; j<nat; j++)
1862       excl->a[k++] = j;
1863   }
1864   excl->index[nat] = k;
1865 }
1866
1867 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
1868                            int couple_lam0,int couple_lam1)
1869 {
1870   int i;
1871
1872   for(i=0; i<atoms->nr; i++) {
1873     if (couple_lam0 != ecouplamVDWQ)
1874       atoms->atom[i].q     = 0.0;
1875     if (couple_lam0 == ecouplamNONE)
1876       atoms->atom[i].type  = atomtype_decouple;
1877     if (couple_lam1 != ecouplamVDWQ)
1878       atoms->atom[i].qB    = 0.0;
1879     if (couple_lam1 == ecouplamNONE)
1880       atoms->atom[i].typeB = atomtype_decouple;
1881   }
1882 }
1883
1884 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
1885                             int couple_lam0,int couple_lam1,
1886                             bool bCoupleIntra,int nb_funct,t_params *nbp)
1887 {
1888   convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
1889   if (!bCoupleIntra) {
1890     generate_LJCpairsNB(mol,nb_funct,nbp);
1891     set_excl_all(&mol->excls);
1892   }
1893   decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
1894 }