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