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,gpp_atomtype_t 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(gpp_atomtype_t 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, gpp_atomtype_t at, t_bond_atomtype bat,
 174               char *line,int nb_funct,
 175               t_nbparam ***nbparam,t_nbparam ***pair)
 176 {
 177   typedef struct {
 178     const 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       = -1;
 266   radius    = -1;
 267   gb_radius = -1;
 268   atomnr    = -1;
 269   S_hct     = -1;
 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                 {
 549             bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
 550         }
 551
 552         bt->nr += 2;
 553     }
 554 }
 555
 556 void push_bt(directive d,t_params bt[],int nral,
 557              gpp_atomtype_t at,
 558              t_bond_atomtype bat,char *line)
 559 {
 560   const char *formal[MAXATOMLIST+1] = {
 561     "%s",
 562     "%s%s",
 563     "%s%s%s",
 564     "%s%s%s%s",
 565     "%s%s%s%s%s",
 566     "%s%s%s%s%s%s",
 567     "%s%s%s%s%s%s%s"
 568   };
 569   const char *formnl[MAXATOMLIST+1] = {
 570     "%*s",
 571     "%*s%*s",
 572     "%*s%*s%*s",
 573     "%*s%*s%*s%*s",
 574     "%*s%*s%*s%*s%*s",
 575     "%*s%*s%*s%*s%*s%*s",
 576     "%*s%*s%*s%*s%*s%*s%*s"
 577   };
 578   const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
 579   int      i,ft,ftype,nn,nrfp,nrfpA,nrfpB;
 580   char     f1[STRLEN];
 581   char     alc[MAXATOMLIST+1][20];
 582   /* One force parameter more, so we can check if we read too many */
 583   double   c[MAXFORCEPARAM+1];
 584   t_param  p;
 585   char  errbuf[256];
 586
 587   if ((bat && at) || (!bat && !at))
 588     gmx_incons("You should pass either bat or at to push_bt");
 589
 590   /* Make format string (nral ints+functype) */
 591   if ((nn=sscanf(line,formal[nral],
 592                  alc[0],alc[1],alc[2],alc[3],alc[4],alc[5])) != nral+1) {
 593     sprintf(errbuf,"Not enough atomtypes (%d instead of %d)",nn-1,nral);
 594     warning_error(errbuf);
 595     return;
 596   }
 597
 598   ft    = strtol(alc[nral],NULL,10);
 599   ftype = ifunc_index(d,ft);
 600   nrfp  = NRFP(ftype);
 601   nrfpA = interaction_function[ftype].nrfpA;
 602   nrfpB = interaction_function[ftype].nrfpB;
 603   strcpy(f1,formnl[nral]);
 604   strcat(f1,formlf);
 605   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]))
 606       != nrfp) {
 607     if (nn == nrfpA) {
 608       /* Copy the B-state from the A-state */
 609       copy_B_from_A(ftype,c);
 610     } else {
 611       if (nn < nrfpA) {
 612       warning_error("Not enough parameters");
 613       } else if (nn > nrfpA && nn < nrfp) {
 614         warning_error("Too many parameters or not enough parameters for topology B");
 615       } else if (nn > nrfp) {
 616         warning_error("Too many parameters");
 617       }
 618       for(i=nn; (i<nrfp); i++)
 619         c[i] = 0.0;
 620     }
 621   }
 622   for(i=0; (i<nral); i++) {
 623     if (at && ((p.a[i]=get_atomtype_type(alc[i],at)) == NOTSET))
 624       gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
 625     else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
 626       gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
 627   }
 628   for(i=0; (i<nrfp); i++)
 629     p.c[i]=c[i];
 630   push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line);
 631 }
 632
 633
 634 void push_dihedraltype(directive d,t_params bt[],
 635                        t_bond_atomtype bat,char *line)
 636 {
 637   const char *formal[MAXATOMLIST+1] = {
 638     "%s",
 639     "%s%s",
 640     "%s%s%s",
 641     "%s%s%s%s",
 642     "%s%s%s%s%s",
 643     "%s%s%s%s%s%s",
 644     "%s%s%s%s%s%s%s"
 645   };
 646   const char *formnl[MAXATOMLIST+1] = {
 647     "%*s",
 648     "%*s%*s",
 649     "%*s%*s%*s",
 650     "%*s%*s%*s%*s",
 651     "%*s%*s%*s%*s%*s",
 652         "%*s%*s%*s%*s%*s%*s",
 653         "%*s%*s%*s%*s%*s%*s%*s"
 654   };
 655   const char *formlf[MAXFORCEPARAM] = {
 656     "%lf",
 657     "%lf%lf",
 658     "%lf%lf%lf",
 659     "%lf%lf%lf%lf",
 660     "%lf%lf%lf%lf%lf",
 661     "%lf%lf%lf%lf%lf%lf",
 662     "%lf%lf%lf%lf%lf%lf%lf",
 663     "%lf%lf%lf%lf%lf%lf%lf%lf",
 664     "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 665     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 666     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 667     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 668   };
 669   int      i,ft,ftype,nn,nrfp,nrfpA,nrfpB,nral;
 670   char     f1[STRLEN];
 671   char     alc[MAXATOMLIST+1][20];
 672   double   c[MAXFORCEPARAM];
 673   t_param  p;
 674   bool     bAllowRepeat;
 675   char  errbuf[256];
 676
 677   /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
 678    *
 679    * We first check for 2 atoms with the 3th column being an integer
 680    * defining the type. If this isn't the case, we try it with 4 atoms
 681    * and the 5th column defining the dihedral type.
 682    */
 683   nn=sscanf(line,formal[4],alc[0],alc[1],alc[2],alc[3],alc[4]);
 684   if(nn>=3 && strlen(alc[2])==1 && isdigit(alc[2][0])) {
 685     nral=2;
 686     ft    = strtol(alc[nral],NULL,10);
 687     /* Move atom types around a bit and use 'X' for wildcard atoms
 688      * to create a 4-atom dihedral definition with arbitrary atoms in
 689      * position 1 and 4.
 690      */
 691     if(alc[2][0]=='2') {
 692       /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
 693       strcpy(alc[3],alc[1]);
 694       sprintf(alc[2],"X");
 695       sprintf(alc[1],"X");
 696       /* alc[0] stays put */
 697     } else {
 698       /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
 699       sprintf(alc[3],"X");
 700       strcpy(alc[2],alc[1]);
 701       strcpy(alc[1],alc[0]);
 702       sprintf(alc[0],"X");
 703     }
 704   } else if(nn==5 && strlen(alc[4])==1 && isdigit(alc[4][0])) {
 705     nral=4;
 706     ft    = strtol(alc[nral],NULL,10);
 707   } else {
 708     sprintf(errbuf,"Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)",nn-1);
 709     warning_error(errbuf);
 710     return;
 711   }
 712
 713         if(ft == 9)
 714         {
 715                 /* Previously, we have always overwritten parameters if e.g. a torsion
 716                  with the same atomtypes occurs on multiple lines. However, CHARMM and
 717                  some other force fields specify multiple dihedrals over some bonds,
 718                  including cosines with multiplicity 6 and somethimes even higher.
 719                  Thus, they cannot be represented with Ryckaert-Bellemans terms.
 720                  To add support for these force fields, Dihedral type 4 is identical to
 721                  normal proper dihedrals, but repeated entries are allowed.
 722                  */
 723                 bAllowRepeat = TRUE;
 724                 ft = 1;
 725         }
 726         else
 727         {
 728                 bAllowRepeat = FALSE;
 729         }
 730
 731
 732   ftype = ifunc_index(d,ft);
 733   nrfp  = NRFP(ftype);
 734   nrfpA = interaction_function[ftype].nrfpA;
 735   nrfpB = interaction_function[ftype].nrfpB;
 736
 737   strcpy(f1,formnl[nral]);
 738   strcat(f1,formlf[nrfp-1]);
 739
 740   /* Check number of parameters given */
 741   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]))
 742       != nrfp) {
 743     if (nn == nrfpA) {
 744       /* Copy the B-state from the A-state */
 745       copy_B_from_A(ftype,c);
 746     } else {
 747       if (nn < nrfpA) {
 748       warning_error("Not enough parameters");
 749       } else if (nn > nrfpA && nn < nrfp) {
 750         warning_error("Too many parameters or not enough parameters for topology B");
 751       } else if (nn > nrfp) {
 752         warning_error("Too many parameters");
 753       }
 754       for(i=nn; (i<nrfp); i++)
 755         c[i] = 0.0;
 756     }
 757   }
 758
 759   for(i=0; (i<4); i++) {
 760     if(!strcmp(alc[i],"X"))
 761       p.a[i]=-1;
 762     else {
 763       if ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET)
 764         gmx_fatal(FARGS,"Unknown bond_atomtype %s",alc[i]);
 765     }
 766   }
 767   for(i=0; (i<nrfp); i++)
 768     p.c[i]=c[i];
 769   /* Always use 4 atoms here, since we created two wildcard atoms
 770    * if there wasn't of them 4 already.
 771    */
 772   push_bondtype (&(bt[ftype]),&p,4,ftype,bAllowRepeat,line);
 773 }
 774
 775
 776 void push_nbt(directive d,t_nbparam **nbt,gpp_atomtype_t atype,
 777               char *pline,int nb_funct)
 778 {
 779   /* swap the atoms */
 780   const char *form2="%*s%*s%*s%lf%lf";
 781   const char *form3="%*s%*s%*s%lf%lf%lf";
 782   const char *form4="%*s%*s%*s%lf%lf%lf%lf";
 783   const char *form5="%*s%*s%*s%lf%lf%lf%lf%lf";
 784   char    a0[80],a1[80];
 785   int     i,f,n,ftype,atnr,nrfp;
 786   double  c[4],dum;
 787   real    cr[4],sig6;
 788   atom_id ai,aj;
 789   t_nbparam *nbp;
 790   bool    bId;
 791   char  errbuf[256];
 792
 793   if (sscanf (pline,"%s%s%d",a0,a1,&f) != 3) {
 794     too_few();
 795     return;
 796   }
 797
 798   ftype=ifunc_index(d,f);
 799
 800   if (ftype != nb_funct) {
 801     sprintf(errbuf,"Trying to add %s while the default nonbond type is %s",
 802             interaction_function[ftype].longname,
 803             interaction_function[nb_funct].longname);
 804     warning_error(errbuf);
 805     return;
 806   }
 807
 808   /* Get the force parameters */
 809   nrfp = NRFP(ftype);
 810   if (ftype == F_LJ14) {
 811     n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
 812     if (n < 2) {
 813       too_few();
 814       return;
 815     }
 816     /* When the B topology parameters are not set,
 817      * copy them from topology A
 818      */
 819     for(i=n; i<nrfp; i++)
 820       c[i] = c[i-2];
 821   }
 822   else if (ftype == F_LJC14_Q) {
 823     n = sscanf(pline,form5,&c[0],&c[1],&c[2],&c[3],&dum);
 824     if (n != 4) {
 825       incorrect_n_param();
 826       return;
 827     }
 828   }
 829   else if (nrfp == 2) {
 830     if (sscanf(pline,form3,&c[0],&c[1],&dum) != 2) {
 831       incorrect_n_param();
 832       return;
 833     }
 834   }
 835   else if (nrfp == 3) {
 836     if (sscanf(pline,form4,&c[0],&c[1],&c[2],&dum) != 3) {
 837       incorrect_n_param();
 838       return;
 839     }
 840   }
 841   else {
 842     gmx_fatal(FARGS,"Number of force parameters for nonbonded interactions is %d"
 843                 " in file %s, line %d",nrfp,__FILE__,__LINE__);
 844   }
 845   for(i=0; (i<nrfp); i++)
 846     cr[i] = c[i];
 847
 848   /* Put the parameters in the matrix */
 849   if ((ai = get_atomtype_type (a0,atype)) == NOTSET)
 850     gmx_fatal(FARGS,"Atomtype %s not found",a0);
 851   if ((aj = get_atomtype_type (a1,atype)) == NOTSET)
 852     gmx_fatal(FARGS,"Atomtype %s not found",a1);
 853   nbp = &(nbt[max(ai,aj)][min(ai,aj)]);
 854
 855   if (nbp->bSet) {
 856     bId = TRUE;
 857     for (i=0; i<nrfp; i++)
 858       bId = bId && (nbp->c[i] == cr[i]);
 859     if (!bId) {
 860       sprintf(errbuf,"Overriding non-bonded parameters,");
 861       warning(errbuf);
 862       fprintf(stderr,"  old:");
 863       for (i=0; i<nrfp; i++)
 864         fprintf(stderr," %g",nbp->c[i]);
 865       fprintf(stderr," new\n%s\n",pline);
 866     }
 867   }
 868   nbp->bSet = TRUE;
 869   for (i=0; i<nrfp; i++)
 870     nbp->c[i] = cr[i];
 871 }
 872
 873 void
 874 push_gb_params (gpp_atomtype_t at, char *line)
 875 {
 876         int nfield;
 877         int i,n,k,found,gfound;
 878         double radius,vol,surftens,gb_radius,S_hct;
 879         char atypename[STRLEN];
 880         char errbuf[STRLEN];
 881
 882         if( (nfield = sscanf(line,"%s%lf%lf%lf%lf%lf",atypename,&radius,&vol,&surftens,&gb_radius,&S_hct)) != 6)
 883     {
 884                 sprintf(errbuf,"Too few gb parameters for type %s\n",atypename);
 885                 warning(errbuf);
 886     }
 887
 888         /* Search for atomtype */
 889         found = 0;
 890         gfound = -1;
 891         for(i=0;i<get_atomtype_ntypes(at) && !found;i++)
 892         {
 893           if(gmx_strncasecmp(atypename,get_atomtype_name(i,at),STRLEN-1)==0)
 894                 {
 895                         found = i;
 896                         gfound = i;
 897                 }
 898     }
 899
 900         if(gfound==-1)
 901     {
 902                 printf("Couldn't find topology match for atomtype %s\n",atypename);
 903                 abort();
 904     }
 905
 906         set_atomtype_gbparam(at,found,radius,vol,surftens,gb_radius,S_hct);
 907 }
 908
 909 void
 910 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
 911               t_bond_atomtype bat, char *line)
 912 {
 913         const char *formal="%s%s%s%s%s%s%s%s";
 914
 915         int      i,j,ft,ftype,nn,nrfp,nrfpA,nrfpB;
 916         int      start;
 917         int      nxcmap,nycmap,ncmap,read_cmap,sl,nct;
 918         char     s[20],alc[MAXATOMLIST+1][20];
 919         t_param  p;
 920         bool     bAllowRepeat;
 921         char     errbuf[256];
 922
 923         /* Keep the compiler happy */
 924         read_cmap = 0;
 925         start     = 0;
 926
 927         if((nn=sscanf(line,formal,alc[0],alc[1],alc[2],alc[3],alc[4],alc[5],alc[6],alc[7])) != nral+3)
 928         {
 929                 sprintf(errbuf,"Incorrect number of atomtypes for cmap (%d instead of 5)",nn-1);
 930                 warning_error(errbuf);
 931                 return;
 932         }
 933
 934         /* Compute an offset for each line where the cmap parameters start
 935          * ie. where the atom types and grid spacing information ends
 936          */
 937         for(i=0;i<nn;i++)
 938         {
 939                 start += (int)strlen(alc[i]);
 940         }
 941
 942         /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */
 943         /* start is the position on the line where we start to read the actual cmap grid data from the itp file */
 944         start = start + nn -1;
 945
 946         ft     = strtol(alc[nral],NULL,10);
 947         nxcmap = strtol(alc[nral+1],NULL,10);
 948         nycmap = strtol(alc[nral+2],NULL,10);
 949
 950         /* Check for equal grid spacing in x and y dims */
 951         if(nxcmap!=nycmap)
 952         {
 953                 gmx_fatal(FARGS,"Not the same grid spacing in x and y for cmap grid: x=%d, y=%d",nxcmap,nycmap);
 954         }
 955
 956         ncmap  = nxcmap*nycmap;
 957         ftype = ifunc_index(d,ft);
 958         nrfpA = strtol(alc[6],NULL,10)*strtol(alc[6],NULL,10);
 959         nrfpB = strtol(alc[7],NULL,10)*strtol(alc[7],NULL,10);
 960         nrfp  = nrfpA+nrfpB;
 961
 962         /* Allocate memory for the CMAP grid */
 963         bt->ncmap+=nrfp;
 964         srenew(bt->cmap,bt->ncmap);
 965
 966         /* Read in CMAP parameters */
 967         sl = 0;
 968         for(i=0;i<ncmap;i++)
 969         {
 970                 nn=sscanf(line+start+sl,"%s",s);
 971                 sl+=strlen(s)+1;
 972                 bt->cmap[i+(bt->ncmap)-nrfp]=strtod(s,NULL);
 973
 974                 if(nn==1)
 975                 {
 976                         read_cmap++;
 977                 }
 978                 else
 979                 {
 980                         gmx_fatal(FARGS,"Error in reading cmap parameter for angle %s %s %s %s %s",alc[0],alc[1],alc[2],alc[3],alc[4]);
 981                 }
 982
 983         }
 984
 985         /* Check do that we got the number of parameters we expected */
 986         if(read_cmap==nrfpA)
 987         {
 988                 for(i=0;i<ncmap;i++)
 989                 {
 990                         bt->cmap[i+ncmap]=bt->cmap[i];
 991                 }
 992         }
 993         else
 994         {
 995                 if(read_cmap<nrfpA)
 996                 {
 997                         warning_error("Not enough cmap parameters");
 998                 }
 999                 else if(read_cmap > nrfpA && read_cmap < nrfp)
1000                 {
1001                         warning_error("Too many cmap parameters or not enough parameters for topology B");
1002                 }
1003                 else if(read_cmap>nrfp)
1004                 {
1005                         warning_error("Too many cmap parameters");
1006                 }
1007         }
1008
1009
1010         /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1011          * so we can safely assign them each time
1012          */
1013         bt->grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
1014         bt->nc           = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1015         nct              = (nral+1) * bt->nc;
1016
1017         /* Allocate memory for the cmap_types information */
1018         srenew(bt->cmap_types,nct);
1019
1020         for(i=0; (i<nral); i++)
1021         {
1022                 if(at && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1023                 {
1024                         gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
1025                 }
1026                 else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1027                 {
1028                         gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
1029                 }
1030
1031                 /* Assign a grid number to each cmap_type */
1032                 bt->cmap_types[bt->nct++]=get_bond_atomtype_type(alc[i],bat);
1033         }
1034
1035         /* Assign a type number to this cmap */
1036         bt->cmap_types[bt->nct++]=bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1037
1038         /* Check for the correct number of atoms (again) */
1039         if(bt->nct!=nct)
1040         {
1041                 gmx_fatal(FARGS,"Incorrect number of atom types (%d) in cmap type %d\n",nct,bt->nc);
1042         }
1043
1044         /* Is this correct?? */
1045         for(i=0;i<MAXFORCEPARAM;i++)
1046         {
1047                 p.c[i]=NOTSET;
1048         }
1049
1050         /* Push the bond to the bondlist */
1051         push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line);
1052 }
1053
1054
1055 static void push_atom_now(t_symtab *symtab,t_atoms *at,int atomnr,
1056                           int atomicnumber,
1057                           int type,char *ctype,int ptype,
1058                           char *resnumberic,int cgnumber,
1059                           char *resname,char *name,real m0,real q0,
1060                           int typeB,char *ctypeB,real mB,real qB)
1061 {
1062   int j,resind=0,resnr;
1063   unsigned char ric;
1064   int nr = at->nr;
1065
1066   if (((nr==0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1067     gmx_fatal(FARGS,"Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)",atomnr,at->nr);
1068
1069   j = strlen(resnumberic) - 1;
1070   if (isdigit(resnumberic[j])) {
1071     ric = ' ';
1072   } else {
1073     ric = resnumberic[j];
1074     if (j == 0 || !isdigit(resnumberic[j-1])) {
1075       gmx_fatal(FARGS,"Invalid residue number '%s' for atom %d",
1076                 resnumberic,atomnr);
1077     }
1078   }
1079   resnr = strtol(resnumberic,NULL,10);
1080
1081   if (nr > 0) {
1082     resind = at->atom[nr-1].resind;
1083   }
1084   if (nr == 0 || strcmp(resname,*at->resinfo[resind].name) != 0 ||
1085       resnr != at->resinfo[resind].nr ||
1086       ric   != at->resinfo[resind].ic) {
1087     if (nr == 0) {
1088       resind = 0;
1089     } else {
1090       resind++;
1091     }
1092     at->nres = resind + 1;
1093     srenew(at->resinfo,at->nres);
1094     at->resinfo[resind].name = put_symtab(symtab,resname);
1095     at->resinfo[resind].nr   = resnr;
1096     at->resinfo[resind].ic   = ric;
1097   } else {
1098     resind = at->atom[at->nr-1].resind;
1099   }
1100
1101   /* New atom instance
1102    * get new space for arrays
1103    */
1104   srenew(at->atom,nr+1);
1105   srenew(at->atomname,nr+1);
1106   srenew(at->atomtype,nr+1);
1107   srenew(at->atomtypeB,nr+1);
1108
1109   /* fill the list */
1110   at->atom[nr].type  = type;
1111   at->atom[nr].ptype = ptype;
1112   at->atom[nr].q     = q0;
1113   at->atom[nr].m     = m0;
1114   at->atom[nr].typeB = typeB;
1115   at->atom[nr].qB    = qB;
1116   at->atom[nr].mB    = mB;
1117
1118   at->atom[nr].resind = resind;
1119   at->atom[nr].atomnumber = atomicnumber;
1120   at->atomname[nr] = put_symtab(symtab,name);
1121   at->atomtype[nr] = put_symtab(symtab,ctype);
1122   at->atomtypeB[nr] = put_symtab(symtab,ctypeB);
1123   at->nr++;
1124 }
1125
1126 void push_cg(t_block *block, int *lastindex, int index, int a)
1127 {
1128   if (debug)
1129     fprintf (debug,"Index %d, Atom %d\n",index,a);
1130
1131   if (((block->nr) && (*lastindex != index)) || (!block->nr)) {
1132     /* add a new block */
1133     block->nr++;
1134     srenew(block->index,block->nr+1);
1135   }
1136   block->index[block->nr] = a + 1;
1137   *lastindex = index;
1138 }
1139
1140 void push_atom(t_symtab *symtab,t_block *cgs,
1141                t_atoms *at,gpp_atomtype_t atype,char *line,int *lastcg)
1142 {
1143   int           nr,ptype;
1144   int           cgnumber,atomnr,type,typeB,nscan;
1145   char          id[STRLEN],ctype[STRLEN],ctypeB[STRLEN],
1146                 resnumberic[STRLEN],resname[STRLEN],name[STRLEN],check[STRLEN];
1147   double        m,q,mb,qb;
1148   real          m0,q0,mB,qB;
1149
1150   /* Make a shortcut for writing in this molecule  */
1151   nr = at->nr;
1152
1153   /* Fixed parameters */
1154   if (sscanf(line,"%s%s%s%s%s%d",
1155              id,ctype,resnumberic,resname,name,&cgnumber) != 6) {
1156     too_few();
1157     return;
1158   }
1159   sscanf(id,"%d",&atomnr);
1160   if ((type  = get_atomtype_type(ctype,atype)) == NOTSET)
1161     gmx_fatal(FARGS,"Atomtype %s not found",ctype);
1162   ptype = get_atomtype_ptype(type,atype);
1163
1164   /* Set default from type */
1165   q0    = get_atomtype_qA(type,atype);
1166   m0    = get_atomtype_massA(type,atype);
1167   typeB = type;
1168   qB    = q0;
1169   mB    = m0;
1170
1171   /* Optional parameters */
1172   nscan = sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1173                  &q,&m,ctypeB,&qb,&mb,check);
1174
1175   /* Nasty switch that falls thru all the way down! */
1176   if (nscan > 0) {
1177     q0 = qB = q;
1178     if (nscan > 1) {
1179       m0 = mB = m;
1180       if (nscan > 2) {
1181         if ((typeB = get_atomtype_type(ctypeB,atype)) == NOTSET) {
1182           gmx_fatal(FARGS,"Atomtype %s not found",ctypeB);
1183         }
1184         qB = get_atomtype_qA(typeB,atype);
1185         mB = get_atomtype_massA(typeB,atype);
1186         if (nscan > 3) {
1187           qB = qb;
1188           if (nscan > 4) {
1189             mB = mb;
1190             if (nscan > 5)
1191               warning_error("Too many parameters");
1192           }
1193         }
1194       }
1195     }
1196   }
1197   if (debug)
1198     fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
1199
1200   push_cg(cgs,lastcg,cgnumber,nr);
1201
1202   push_atom_now(symtab,at,atomnr,get_atomtype_atomnumber(type,atype),
1203                 type,ctype,ptype,resnumberic,cgnumber,
1204                 resname,name,m0,q0,typeB,
1205                 typeB==type ? ctype : ctypeB,mB,qB);
1206 }
1207
1208 void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line)
1209 {
1210   char type[STRLEN];
1211   int nrexcl,i;
1212   t_molinfo *newmol;
1213
1214   if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
1215     warning_error("Expected a molecule type name and nrexcl");
1216   }
1217
1218   /* Test if this atomtype overwrites another */
1219   i = 0;
1220   while (i < *nmol) {
1221     if (strcasecmp(*((*mol)[i].name),type) == 0)
1222       gmx_fatal(FARGS,"moleculetype %s is redefined",type);
1223     i++;
1224   }
1225
1226   (*nmol)++;
1227   srenew(*mol,*nmol);
1228   newmol = &((*mol)[*nmol-1]);
1229   init_molinfo(newmol);
1230
1231   /* Fill in the values */
1232   newmol->name     = put_symtab(symtab,type);
1233   newmol->nrexcl   = nrexcl;
1234   newmol->excl_set = FALSE;
1235 }
1236
1237 static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
1238                               t_param *p,int c_start,bool bB,bool bGenPairs)
1239 {
1240   int      i,j,ti,tj,ntype;
1241   bool     bFound;
1242   t_param  *pi=NULL;
1243   int      nr   = bt[ftype].nr;
1244   int      nral = NRAL(ftype);
1245   int      nrfp = interaction_function[ftype].nrfpA;
1246   int      nrfpB= interaction_function[ftype].nrfpB;
1247
1248   if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1249     return TRUE;
1250
1251   bFound=FALSE;
1252   if(bGenPairs) {
1253     /* First test the generated-pair position to save
1254      * time when we have 1000*1000 entries for e.g. OPLS...
1255      */
1256     ntype=sqrt(nr);
1257     if (bB) {
1258       ti=at->atom[p->a[0]].typeB;
1259       tj=at->atom[p->a[1]].typeB;
1260     } else {
1261       ti=at->atom[p->a[0]].type;
1262       tj=at->atom[p->a[1]].type;
1263     }
1264     pi=&(bt[ftype].param[ntype*ti+tj]);
1265     bFound=((ti == pi->a[0]) && (tj == pi->a[1]));
1266   }
1267
1268   /* Search explicitly if we didnt find it */
1269   if(!bFound) {
1270     for (i=0; ((i < nr) && !bFound); i++) {
1271       pi=&(bt[ftype].param[i]);
1272       if (bB)
1273         for (j=0; ((j < nral) &&
1274                    (at->atom[p->a[j]].typeB == pi->a[j])); j++);
1275       else
1276         for (j=0; ((j < nral) &&
1277                    (at->atom[p->a[j]].type == pi->a[j])); j++);
1278       bFound=(j==nral);
1279     }
1280   }
1281
1282   if (bFound) {
1283     if (bB) {
1284       if (nrfp+nrfpB > MAXFORCEPARAM) {
1285         gmx_incons("Too many force parameters");
1286       }
1287       for (j=c_start; (j < nrfpB); j++)
1288         p->c[nrfp+j] = pi->c[j];
1289     }
1290     else
1291       for (j=c_start; (j < nrfp); j++)
1292         p->c[j] = pi->c[j];
1293   }
1294   else {
1295     for (j=c_start; (j < nrfp); j++)
1296       p->c[j] = 0.0;
1297   }
1298   return bFound;
1299 }
1300
1301 static bool default_cmap_params(int ftype, t_params bondtype[],
1302                                 t_atoms *at, gpp_atomtype_t atype,
1303                                 t_param *p, bool bB,
1304                                 int *cmap_type, int *nparam_def)
1305 {
1306         int i,j,nparam_found;
1307         int ct;
1308         bool bFound=FALSE;
1309
1310         nparam_found = 0;
1311         ct           = 0;
1312
1313         /* Match the current cmap angle against the list of cmap_types */
1314         for(i=0;i<bondtype->nct && !bFound;i+=6)
1315         {
1316                 if(bB)
1317                 {
1318
1319                 }
1320                 else
1321                 {
1322                         if(
1323                            (get_atomtype_batype(at->atom[p->a[0]].type,atype)==bondtype->cmap_types[i])   &&
1324                            (get_atomtype_batype(at->atom[p->a[1]].type,atype)==bondtype->cmap_types[i+1]) &&
1325                            (get_atomtype_batype(at->atom[p->a[2]].type,atype)==bondtype->cmap_types[i+2]) &&
1326                            (get_atomtype_batype(at->atom[p->a[3]].type,atype)==bondtype->cmap_types[i+3]) &&
1327                            (get_atomtype_batype(at->atom[p->a[4]].type,atype)==bondtype->cmap_types[i+4]))
1328                         {
1329                                 /* Found cmap torsion */
1330                                 bFound=TRUE;
1331                                 ct = bondtype->cmap_types[i+5];
1332                                 nparam_found=1;
1333                         }
1334                 }
1335         }
1336
1337         /* If we did not find a matching type for this cmap torsion */
1338         if(!bFound)
1339         {
1340                 gmx_fatal(FARGS,"Unknown cmap torsion between atoms %d %d %d %d %d\n",
1341                                   p->a[0]+1,p->a[1]+1,p->a[2]+1,p->a[3]+1,p->a[4]+1);
1342         }
1343
1344         *nparam_def = nparam_found;
1345         *cmap_type  = ct;
1346
1347         return bFound;
1348 }
1349
1350 static bool default_params(int ftype,t_params bt[],
1351                            t_atoms *at,gpp_atomtype_t atype,
1352                            t_param *p,bool bB,
1353                            t_param **param_def,
1354                int *nparam_def)
1355 {
1356   int      i,j,nparam_found;
1357   bool     bFound,bSame;
1358   t_param  *pi=NULL;
1359   t_param  *pj=NULL;
1360   int      nr   = bt[ftype].nr;
1361   int      nral = NRAL(ftype);
1362   int      nrfpA= interaction_function[ftype].nrfpA;
1363   int      nrfpB= interaction_function[ftype].nrfpB;
1364
1365   if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1366     return TRUE;
1367
1368
1369   /* We allow wildcards now. The first type (with or without wildcards) that
1370    * fits is used, so you should probably put the wildcarded bondtypes
1371    * at the end of each section.
1372    */
1373   bFound=FALSE;
1374   nparam_found=0;
1375   /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1376    * special case for this. Check for B state outside loop to speed it up.
1377    */
1378   if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS)
1379   {
1380           if (bB)
1381           {
1382                   for (i=0; ((i < nr) && !bFound); i++)
1383                   {
1384                           pi=&(bt[ftype].param[i]);
1385                           bFound=
1386                           (
1387                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1388                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1389                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1390                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1391                            );
1392                   }
1393           }
1394           else
1395           {
1396                   /* State A */
1397                   for (i=0; ((i < nr) && !bFound); i++)
1398                   {
1399                           pi=&(bt[ftype].param[i]);
1400                           bFound=
1401                           (
1402                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].type,atype)==pi->AI)) &&
1403                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].type,atype)==pi->AJ)) &&
1404                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].type,atype)==pi->AK)) &&
1405                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].type,atype)==pi->AL))
1406                            );
1407                   }
1408           }
1409           /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1410            * The rules in that case is that additional matches HAVE to be on adjacent lines!
1411            */
1412           if(bFound==TRUE)
1413           {
1414                   nparam_found++;
1415                   bSame = TRUE;
1416                   /* Continue from current i value */
1417                   for(j=i+1 ; j<nr && bSame ; j+=2)
1418                   {
1419                           pj=&(bt[ftype].param[j]);
1420                           bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1421                           if(bSame)
1422                                   nparam_found++;
1423                           /* nparam_found will be increased as long as the numbers match */
1424                   }
1425           }
1426   } else { /* Not a dihedral */
1427           for (i=0; ((i < nr) && !bFound); i++) {
1428                   pi=&(bt[ftype].param[i]);
1429                   if (bB)
1430                           for (j=0; ((j < nral) &&
1431                                                  (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
1432                                   ;
1433                   else
1434                           for (j=0; ((j < nral) &&
1435                                                  (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
1436                                   ;
1437                   bFound=(j==nral);
1438           }
1439           if(bFound)
1440                   nparam_found=1;
1441   }
1442
1443   *param_def = pi;
1444   *nparam_def = nparam_found;
1445
1446   return bFound;
1447 }
1448
1449
1450
1451 void push_bond(directive d,t_params bondtype[],t_params bond[],
1452                t_atoms *at,gpp_atomtype_t atype,char *line,
1453                bool bBonded,bool bGenPairs,real fudgeQQ,
1454                bool bZero,bool *bWarn_copy_A_B)
1455 {
1456   const char *aaformat[MAXATOMLIST]= {
1457     "%d%d",
1458     "%d%d%d",
1459     "%d%d%d%d",
1460     "%d%d%d%d%d",
1461     "%d%d%d%d%d%d",
1462     "%d%d%d%d%d%d%d"
1463   };
1464   const char *asformat[MAXATOMLIST]= {
1465     "%*s%*s",
1466     "%*s%*s%*s",
1467     "%*s%*s%*s%*s",
1468     "%*s%*s%*s%*s%*s",
1469     "%*s%*s%*s%*s%*s%*s",
1470     "%*s%*s%*s%*s%*s%*s%*s"
1471   };
1472   const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1473   int      nr,i,j,nral,nread,ftype;
1474   char     format[STRLEN];
1475   /* One force parameter more, so we can check if we read too many */
1476   double   cc[MAXFORCEPARAM+1];
1477   int      aa[MAXATOMLIST+1];
1478   t_param  param,paramB,*param_defA,*param_defB;
1479   bool     bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
1480   int      nparam_defA,nparam_defB;
1481   char  errbuf[256];
1482
1483   nparam_defA=nparam_defB=0;
1484
1485   ftype = ifunc_index(d,1);
1486   nral  = NRAL(ftype);
1487   for(j=0; j<MAXATOMLIST; j++)
1488     aa[j]=NOTSET;
1489   bDef = (NRFP(ftype) > 0);
1490
1491   nread = sscanf(line,aaformat[nral-1],
1492                  &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1493   if (nread < nral) {
1494     too_few();
1495     return;
1496   } else if (nread == nral)
1497     ftype = ifunc_index(d,1);
1498   else {
1499     /* this is a hack to allow for virtual sites with swapped parity */
1500     bSwapParity = (aa[nral]<0);
1501     if (bSwapParity)
1502       aa[nral] = -aa[nral];
1503     ftype = ifunc_index(d,aa[nral]);
1504     if (bSwapParity)
1505       switch(ftype) {
1506       case F_VSITE3FAD:
1507       case F_VSITE3OUT:
1508         break;
1509       default:
1510         gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1511                     interaction_function[F_VSITE3FAD].longname,
1512                     interaction_function[F_VSITE3OUT].longname);
1513       }
1514   }
1515
1516   /* Check for double atoms and atoms out of bounds */
1517   for(i=0; (i<nral); i++) {
1518     if ( aa[i] < 1 || aa[i] > at->nr )
1519       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1520                 "Atom index (%d) in %s out of bounds (1-%d).\n"
1521                 "This probably means that you have inserted topology section \"%s\"\n"
1522                 "in a part belonging to a different molecule than you intended to.\n"
1523                 "In that case move the \"%s\" section to the right molecule.",
1524                 get_warning_file(),get_warning_line(),
1525                 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1526     for(j=i+1; (j<nral); j++)
1527       if (aa[i] == aa[j]) {
1528         sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1529         warning(errbuf);
1530       }
1531   }
1532   if (ftype == F_SETTLE)
1533     if (aa[0]+2 > at->nr)
1534       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1535                   "             Atom index (%d) in %s out of bounds (1-%d)\n"
1536                   "             Settle works on atoms %d, %d and %d",
1537                   get_warning_file(),get_warning_line(),
1538                   aa[0],dir2str(d),at->nr,
1539                   aa[0],aa[0]+1,aa[0]+2);
1540
1541   /* default force parameters  */
1542   for(j=0; (j<MAXATOMLIST); j++)
1543     param.a[j] = aa[j]-1;
1544   for(j=0; (j<MAXFORCEPARAM); j++)
1545     param.c[j] = 0.0;
1546
1547   /* Get force params for normal and free energy perturbation
1548    * studies, as determined by types!
1549    */
1550   if (bBonded) {
1551     bFoundA = default_params(ftype,bondtype,at,atype,&param,FALSE,&param_defA,&nparam_defA);
1552     if (bFoundA) {
1553       /* Copy the A-state and B-state default parameters */
1554       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1555         param.c[j] = param_defA->c[j];
1556     }
1557     bFoundB = default_params(ftype,bondtype,at,atype,&param,TRUE,&param_defB,&nparam_defB);
1558     if (bFoundB) {
1559       /* Copy only the B-state default parameters */
1560       for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1561         param.c[j] = param_defB->c[j];
1562     }
1563   } else if (ftype == F_LJ14) {
1564     bFoundA = default_nb_params(ftype, bondtype,at,&param,0,FALSE,bGenPairs);
1565     bFoundB = default_nb_params(ftype, bondtype,at,&param,0,TRUE, bGenPairs);
1566   } else if (ftype == F_LJC14_Q) {
1567     param.c[0] = fudgeQQ;
1568     /* Fill in the A-state charges as default parameters */
1569     param.c[1] = at->atom[param.a[0]].q;
1570     param.c[2] = at->atom[param.a[1]].q;
1571     /* The default LJ parameters are the standard 1-4 parameters */
1572     bFoundA = default_nb_params(F_LJ14,bondtype,at,&param,3,FALSE,bGenPairs);
1573     bFoundB = TRUE;
1574   } else if (ftype == F_LJC_PAIRS_NB) {
1575     /* Defaults are not supported here */
1576     bFoundA = FALSE;
1577     bFoundB = TRUE;
1578   } else {
1579     gmx_incons("Unknown function type in push_bond");
1580   }
1581
1582   if (nread > nral) {
1583           /* Manually specified parameters - in this case we discard multiple torsion info! */
1584
1585     strcpy(format,asformat[nral-1]);
1586     strcat(format,ccformat);
1587
1588     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1589                    &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1590
1591     if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1592       {
1593         /* We only have to issue a warning if these atoms are perturbed! */
1594         bPert = FALSE;
1595         for(j=0; (j<nral); j++)
1596           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1597
1598         if (bPert && *bWarn_copy_A_B)
1599           {
1600             sprintf(errbuf,
1601                     "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1602             warning(errbuf);
1603             *bWarn_copy_A_B = FALSE;
1604           }
1605
1606         /* If only the A parameters were specified, copy them to the B state */
1607         /* The B-state parameters correspond to the first nrfpB
1608          * A-state parameters.
1609          */
1610         for(j=0; (j<NRFPB(ftype)); j++)
1611           cc[nread++] = cc[j];
1612       }
1613
1614     /* If nread was 0 or EOF, no parameters were read => use defaults.
1615      * If nread was nrfpA we copied above so nread=nrfp.
1616      * If nread was nrfp we are cool.
1617      * For F_LJC14_Q we allow supplying fudgeQQ only.
1618      * Anything else is an error!
1619      */
1620     if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1621         !(ftype == F_LJC14_Q && nread == 1))
1622       {
1623         gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1624                   nread,NRFPA(ftype),NRFP(ftype),
1625                   interaction_function[ftype].longname);
1626       }
1627
1628     for(j=0; (j<nread); j++)
1629       param.c[j]=cc[j];
1630
1631     /* Check whether we have to use the defaults */
1632     if (nread == NRFP(ftype))
1633       bDef = FALSE;
1634   }
1635   else
1636   {
1637     nread = 0;
1638   }
1639   /* nread now holds the number of force parameters read! */
1640
1641
1642
1643
1644         if (bDef)
1645         {
1646           /* Use defaults */
1647                 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1648                 if(ftype==F_PDIHS)
1649                 {
1650                         if ((nparam_defA != nparam_defB) || ((nparam_defA>1 || nparam_defB>1) && (param_defA!=param_defB)))
1651                         {
1652                                 sprintf(errbuf,
1653                                                 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1654                                                 "Please specify perturbed parameters manually for this torsion in your topology!");
1655                                 warning_error(errbuf);
1656                         }
1657                 }
1658
1659           if (nread > 0 && nread < NRFPA(ftype))
1660           {
1661                   /* Issue an error, do not use defaults */
1662                   sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1663                   warning_error(errbuf);
1664           }
1665
1666       if (nread == 0 || nread == EOF)
1667           {
1668                   if (!bFoundA)
1669                   {
1670                           if (interaction_function[ftype].flags & IF_VSITE)
1671                           {
1672                                   /* set them to NOTSET, will be calculated later */
1673                                   for(j=0; (j<MAXFORCEPARAM); j++)
1674                                           param.c[j] = NOTSET;
1675
1676                                   if (bSwapParity)
1677                                           param.C1 = -1; /* flag to swap parity of vsite construction */
1678                           }
1679                           else
1680                           {
1681                                   if (bZero)
1682                                   {
1683                                           fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1684                                                           interaction_function[ftype].longname);
1685                                   }
1686                                   else
1687                                   {
1688                                           sprintf(errbuf,"No default %s types",interaction_function[ftype].longname);
1689                                           warning_error(errbuf);
1690                                   }
1691                           }
1692                   }
1693                   else
1694                   {
1695                           if (bSwapParity)
1696                           {
1697                                   switch(ftype)
1698                                   {
1699                                           case F_VSITE3FAD:
1700                                                   param.C0 = 360 - param.C0;
1701                                                   break;
1702                                           case F_VSITE3OUT:
1703                                                   param.C2 = -param.C2;
1704                                                   break;
1705                                   }
1706                           }
1707                   }
1708                   if (!bFoundB)
1709                   {
1710                           /* We only have to issue a warning if these atoms are perturbed! */
1711                           bPert = FALSE;
1712                           for(j=0; (j<nral); j++)
1713                                   bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1714
1715                           if (bPert)
1716                           {
1717                                   sprintf(errbuf,"No default %s types for perturbed atoms, "
1718                                                   "using normal values",interaction_function[ftype].longname);
1719                                   warning(errbuf);
1720                           }
1721                   }
1722           }
1723   }
1724
1725         if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1726                 && param.c[5]!=param.c[2])
1727                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1728                                   "             %s multiplicity can not be perturbed %f!=%f",
1729                                   get_warning_file(),get_warning_line(),
1730                                   interaction_function[ftype].longname,
1731                                   param.c[2],param.c[5]);
1732
1733         if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1734                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1735                                   "             %s table number can not be perturbed %d!=%d",
1736                                   get_warning_file(),get_warning_line(),
1737                                   interaction_function[ftype].longname,
1738                                   (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1739
1740         /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1741         if (ftype==F_RBDIHS) {
1742                 nr=0;
1743                 for(i=0;i<NRFP(ftype);i++) {
1744                         if(param.c[i]!=0)
1745                                 nr++;
1746                 }
1747                 if(nr==0)
1748                         return;
1749         }
1750
1751         /* Put the values in the appropriate arrays */
1752         add_param_to_list (&bond[ftype],&param);
1753
1754         /* Push additional torsions from FF for ftype==9 if we have them.
1755          * We have already checked that the A/B states do not differ in this case,
1756          * so we do not have to double-check that again, or the vsite stuff.
1757          * In addition, those torsions cannot be automatically perturbed.
1758          */
1759         if(bDef && ftype==F_PDIHS)
1760         {
1761                 for(i=1;i<nparam_defA;i++)
1762                 {
1763                         /* Advance pointer! */
1764                         param_defA+=2;
1765                         for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1766                                 param.c[j] = param_defA->c[j];
1767                         /* And push the next term for this torsion */
1768                         add_param_to_list (&bond[ftype],&param);
1769                 }
1770         }
1771 }
1772
1773 void push_cmap(directive d, t_params bondtype[], t_params bond[],
1774                            t_atoms *at, gpp_atomtype_t atype, char *line,
1775                            bool *bWarn_copy_A_B)
1776 {
1777         const char *aaformat[MAXATOMLIST+1]=
1778         {
1779                 "%d",
1780                 "%d%d",
1781                 "%d%d%d",
1782                 "%d%d%d%d",
1783                 "%d%d%d%d%d",
1784                 "%d%d%d%d%d%d",
1785                 "%d%d%d%d%d%d%d"
1786         };
1787
1788         int i,j,nr,ftype,nral,nread,ncmap_params;
1789         int cmap_type;
1790         int aa[MAXATOMLIST+1];
1791         char  errbuf[256];
1792         bool bFound;
1793         t_param param,paramB,*param_defA,*param_defB;
1794
1795         ftype        = ifunc_index(d,1);
1796         nral         = NRAL(ftype);
1797         nr           = bondtype[ftype].nr;
1798         ncmap_params = 0;
1799
1800         nread = sscanf(line,aaformat[nral-1],
1801                                    &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1802
1803         if (nread < nral)
1804         {
1805                 too_few();
1806                 return;
1807         }
1808         else if (nread == nral)
1809         {
1810                 ftype = ifunc_index(d,1);
1811         }
1812
1813         /* Check for double atoms and atoms out of bounds */
1814         for(i=0;i<nral;i++)
1815         {
1816                 if(aa[i] < 1 || aa[i] > at->nr)
1817                 {
1818                         gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1819                                           "Atom index (%d) in %s out of bounds (1-%d).\n"
1820                                           "This probably means that you have inserted topology section \"%s\"\n"
1821                                           "in a part belonging to a different molecule than you intended to.\n"
1822                                           "In that case move the \"%s\" section to the right molecule.",
1823                                           get_warning_file(),get_warning_line(),
1824                                           aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1825                 }
1826
1827                 for(j=i+1; (j<nral); j++)
1828                 {
1829                         if (aa[i] == aa[j])
1830                         {
1831                                 sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1832                                 warning(errbuf);
1833                         }
1834                 }
1835         }
1836
1837         /* default force parameters  */
1838         for(j=0; (j<MAXATOMLIST); j++)
1839         {
1840                 param.a[j] = aa[j]-1;
1841         }
1842         for(j=0; (j<MAXFORCEPARAM); j++)
1843         {
1844                 param.c[j] = 0.0;
1845         }
1846
1847         /* Get the cmap type for this cmap angle */
1848         bFound = default_cmap_params(ftype,bondtype,at,atype,&param,FALSE,&cmap_type,&ncmap_params);
1849
1850         /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
1851         if(bFound && ncmap_params==1)
1852         {
1853                 /* Put the values in the appropriate arrays */
1854                 param.c[0]=cmap_type;
1855                 add_param_to_list(&bond[ftype],&param);
1856         }
1857         else
1858         {
1859                 /* This is essentially the same check as in default_cmap_params() done one more time */
1860                 gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
1861                                   param.a[0]+1,param.a[1]+1,param.a[2]+1,param.a[3]+1,param.a[4]+1);
1862         }
1863 }
1864
1865
1866
1867 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1868                   t_atoms *at,gpp_atomtype_t atype,char *line)
1869 {
1870   char *ptr;
1871   int  type,ftype,j,n,ret,nj,a;
1872   int  *atc=NULL;
1873   double *weight=NULL,weight_tot;
1874   t_param param;
1875
1876   /* default force parameters  */
1877   for(j=0; (j<MAXATOMLIST); j++)
1878     param.a[j] = NOTSET;
1879   for(j=0; (j<MAXFORCEPARAM); j++)
1880     param.c[j] = 0.0;
1881
1882   ptr = line;
1883   ret = sscanf(ptr,"%d%n",&a,&n);
1884   ptr += n;
1885   if (ret == 0)
1886     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1887               "             Expected an atom index in section \"%s\"",
1888                 get_warning_file(),get_warning_line(),
1889               dir2str(d));
1890
1891   param.a[0] = a - 1;
1892
1893   ret = sscanf(ptr,"%d%n",&type,&n);
1894   ptr += n;
1895   ftype = ifunc_index(d,type);
1896
1897   weight_tot = 0;
1898   nj = 0;
1899   do {
1900     ret = sscanf(ptr,"%d%n",&a,&n);
1901     ptr += n;
1902     if (ret > 0) {
1903       if (nj % 20 == 0) {
1904         srenew(atc,nj+20);
1905         srenew(weight,nj+20);
1906       }
1907       atc[nj] = a - 1;
1908       switch (type) {
1909       case 1:
1910         weight[nj] = 1;
1911         break;
1912       case 2:
1913         /* Here we use the A-state mass as a parameter.
1914          * Note that the B-state mass has no influence.
1915          */
1916         weight[nj] = at->atom[atc[nj]].m;
1917         break;
1918       case 3:
1919         weight[nj] = -1;
1920         ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1921         ptr += n;
1922         if (weight[nj] < 0)
1923           gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1924                     "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1925                     get_warning_file(),get_warning_line(),
1926                     nj+1,atc[nj]+1);
1927         break;
1928       default:
1929         gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1930       }
1931       weight_tot += weight[nj];
1932       nj++;
1933     }
1934   } while (ret > 0);
1935
1936   if (nj == 0)
1937     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1938               "             Expected more than one atom index in section \"%s\"",
1939               get_warning_file(),get_warning_line(),
1940               dir2str(d));
1941
1942   if (weight_tot == 0)
1943      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1944                "             The total mass of the construting atoms is zero",
1945                get_warning_file(),get_warning_line());
1946
1947   for(j=0; j<nj; j++) {
1948     param.a[1] = atc[j];
1949     param.c[0] = nj;
1950     param.c[1] = weight[j]/weight_tot;
1951     /* Put the values in the appropriate arrays */
1952     add_param_to_list (&bond[ftype],&param);
1953   }
1954
1955   sfree(atc);
1956   sfree(weight);
1957 }
1958
1959 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1960               int *nrcopies)
1961 {
1962   int  i,copies;
1963   char type[STRLEN];
1964
1965   *nrcopies=0;
1966   if (sscanf(pline,"%s%d",type,&copies) != 2) {
1967     too_few();
1968     return;
1969   }
1970
1971   /* search moleculename */
1972   for (i=0; ((i<nrmols) && strcasecmp(type,*(mols[i].name))); i++)
1973     ;
1974
1975   if (i<nrmols) {
1976     *nrcopies        = copies;
1977     *whichmol        = i;
1978   } else
1979     gmx_fatal(FARGS,"No such moleculetype %s",type);
1980 }
1981
1982 void init_block2(t_block2 *b2, int natom)
1983 {
1984   int i;
1985
1986   b2->nr=natom;
1987   snew(b2->nra,b2->nr);
1988   snew(b2->a,b2->nr);
1989   for(i=0; (i<b2->nr); i++)
1990     b2->a[i]=NULL;
1991 }
1992
1993 void done_block2(t_block2 *b2)
1994 {
1995   int i;
1996
1997   if (b2->nr) {
1998     for(i=0; (i<b2->nr); i++)
1999       sfree(b2->a[i]);
2000     sfree(b2->a);
2001     sfree(b2->nra);
2002     b2->nr=0;
2003   }
2004 }
2005
2006 void push_excl(char *line, t_block2 *b2)
2007 {
2008   int  i,j;
2009   int  n;
2010   char base[STRLEN],format[STRLEN];
2011
2012   if (sscanf(line,"%d",&i)==0)
2013     return;
2014
2015   if ((1 <= i) && (i <= b2->nr))
2016     i--;
2017   else {
2018     if (debug)
2019       fprintf(debug,"Unbound atom %d\n",i-1);
2020     return;
2021   }
2022   strcpy(base,"%*d");
2023   do {
2024     strcpy(format,base);
2025     strcat(format,"%d");
2026     n=sscanf(line,format,&j);
2027     if (n == 1) {
2028       if ((1 <= j) && (j <= b2->nr)) {
2029         j--;
2030         srenew(b2->a[i],++(b2->nra[i]));
2031         b2->a[i][b2->nra[i]-1]=j;
2032         /* also add the reverse exclusion! */
2033         srenew(b2->a[j],++(b2->nra[j]));
2034         b2->a[j][b2->nra[j]-1]=i;
2035         strcat(base,"%*d");
2036       }
2037       else
2038         gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
2039     }
2040   } while (n == 1);
2041 }
2042
2043 void b_to_b2(t_blocka *b, t_block2 *b2)
2044 {
2045   int     i;
2046   atom_id j,a;
2047
2048   for(i=0; (i<b->nr); i++)
2049     for(j=b->index[i]; (j<b->index[i+1]); j++) {
2050       a=b->a[j];
2051       srenew(b2->a[i],++b2->nra[i]);
2052       b2->a[i][b2->nra[i]-1]=a;
2053     }
2054 }
2055
2056 void b2_to_b(t_block2 *b2, t_blocka *b)
2057 {
2058   int     i,nra;
2059   atom_id j;
2060
2061   nra=0;
2062   for(i=0; (i<b2->nr); i++) {
2063     b->index[i]=nra;
2064     for(j=0; (j<b2->nra[i]); j++)
2065       b->a[nra+j]=b2->a[i][j];
2066     nra+=b2->nra[i];
2067   }
2068   /* terminate list */
2069   b->index[i]=nra;
2070 }
2071
2072 static int icomp(const void *v1, const void *v2)
2073 {
2074   return (*((atom_id *) v1))-(*((atom_id *) v2));
2075 }
2076
2077 void merge_excl(t_blocka *excl, t_block2 *b2)
2078 {
2079   int     i,k;
2080   atom_id j;
2081   int     nra;
2082
2083   if (!b2->nr)
2084     return;
2085   else if (b2->nr != excl->nr) {
2086     gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2087                 b2->nr,excl->nr);
2088   }
2089   else if (debug)
2090     fprintf(debug,"Entering merge_excl\n");
2091
2092   /* First copy all entries from excl to b2 */
2093   b_to_b2(excl,b2);
2094
2095   /* Count and sort the exclusions */
2096   nra=0;
2097   for(i=0; (i<b2->nr); i++) {
2098     if (b2->nra[i] > 0) {
2099       /* remove double entries */
2100       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
2101       k=1;
2102       for(j=1; (j<b2->nra[i]); j++)
2103         if (b2->a[i][j]!=b2->a[i][k-1]) {
2104           b2->a[i][k]=b2->a[i][j];
2105           k++;
2106         }
2107       b2->nra[i]=k;
2108       nra+=b2->nra[i];
2109     }
2110   }
2111   excl->nra=nra;
2112   srenew(excl->a,excl->nra);
2113
2114   b2_to_b(b2,excl);
2115 }
2116
2117 int add_atomtype_decoupled(t_symtab *symtab,gpp_atomtype_t at,
2118                            t_nbparam ***nbparam,t_nbparam ***pair)
2119 {
2120   t_atom  atom;
2121   t_param param;
2122   int     i,nr;
2123
2124   /* Define an atom type with all parameters set to zero (no interactions) */
2125   atom.q = 0.0;
2126   atom.m = 0.0;
2127   /* Type for decoupled atoms could be anything,
2128    * this should be changed automatically later when required.
2129    */
2130   atom.ptype = eptAtom;
2131   for (i=0; (i<MAXFORCEPARAM); i++)
2132     param.c[i] = 0.0;
2133
2134   nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0,0,0);
2135
2136   /* Add space in the non-bonded parameters matrix */
2137   realloc_nb_params(at,nbparam,pair);
2138
2139   return nr;
2140 }
2141
2142 static void convert_pairs_to_pairsQ(t_params *plist,
2143                                     real fudgeQQ,t_atoms *atoms)
2144 {
2145   t_param *param;
2146   int  i;
2147   real v,w;
2148
2149   /* Copy the pair list to the pairQ list */
2150   plist[F_LJC14_Q] = plist[F_LJ14];
2151   /* Empty the pair list */
2152   plist[F_LJ14].nr    = 0;
2153   plist[F_LJ14].param = NULL;
2154   param = plist[F_LJC14_Q].param;
2155   for(i=0; i<plist[F_LJC14_Q].nr; i++) {
2156     v = param[i].c[0];
2157     w = param[i].c[1];
2158     param[i].c[0] = fudgeQQ;
2159     param[i].c[1] = atoms->atom[param[i].a[0]].q;
2160     param[i].c[2] = atoms->atom[param[i].a[1]].q;
2161     param[i].c[3] = v;
2162     param[i].c[4] = w;
2163   }
2164 }
2165
2166 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
2167 {
2168   int n,ntype,i,j,k;
2169   t_atom *atom;
2170   t_blocka *excl;
2171   bool bExcl;
2172   t_param param;
2173
2174   n = mol->atoms.nr;
2175   atom = mol->atoms.atom;
2176
2177   ntype = sqrt(nbp->nr);
2178
2179   for (i=0; i<MAXATOMLIST; i++)
2180     param.a[i] = NOTSET;
2181   for (i=0; i<MAXFORCEPARAM; i++)
2182     param.c[i] = NOTSET;
2183
2184   /* Add a pair interaction for all non-excluded atom pairs */
2185   excl = &mol->excls;
2186   for(i=0; i<n; i++) {
2187     for(j=i+1; j<n; j++) {
2188       bExcl = FALSE;
2189       for(k=excl->index[i]; k<excl->index[i+1]; k++) {
2190         if (excl->a[k] == j)
2191           bExcl = TRUE;
2192       }
2193       if (!bExcl) {
2194         if (nb_funct != F_LJ)
2195           gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2196         param.a[0] = i;
2197         param.a[1] = j;
2198         param.c[0] = atom[i].q;
2199         param.c[1] = atom[j].q;
2200         param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2201         param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2202         add_param_to_list(&mol->plist[F_LJC_PAIRS_NB],&param);
2203       }
2204     }
2205   }
2206 }
2207
2208 static void set_excl_all(t_blocka *excl)
2209 {
2210   int nat,i,j,k;
2211
2212   /* Get rid of the current exclusions and exclude all atom pairs */
2213   nat = excl->nr;
2214   excl->nra = nat*nat;
2215   srenew(excl->a,excl->nra);
2216   k = 0;
2217   for(i=0; i<nat; i++) {
2218     excl->index[i] = k;
2219     for(j=0; j<nat; j++)
2220       excl->a[k++] = j;
2221   }
2222   excl->index[nat] = k;
2223 }
2224
2225 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
2226                            int couple_lam0,int couple_lam1)
2227 {
2228   int i;
2229
2230   for(i=0; i<atoms->nr; i++) {
2231     if (couple_lam0 != ecouplamVDWQ)
2232       atoms->atom[i].q     = 0.0;
2233     if (couple_lam0 == ecouplamNONE)
2234       atoms->atom[i].type  = atomtype_decouple;
2235     if (couple_lam1 != ecouplamVDWQ)
2236       atoms->atom[i].qB    = 0.0;
2237     if (couple_lam1 == ecouplamNONE)
2238       atoms->atom[i].typeB = atomtype_decouple;
2239   }
2240 }
2241
2242 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
2243                             int couple_lam0,int couple_lam1,
2244                             bool bCoupleIntra,int nb_funct,t_params *nbp)
2245 {
2246   convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
2247   if (!bCoupleIntra) {
2248     generate_LJCpairsNB(mol,nb_funct,nbp);
2249     set_excl_all(&mol->excls);
2250   }
2251   decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
2252 }