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