src/gromacs/gmxpreprocess/toppush.c

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