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,nral_fmt,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     {
1506         aa[j] = NOTSET;
1507     }
1508     bDef = (NRFP(ftype) > 0);
1509
1510     if (ftype == F_SETTLE)
1511     {
1512         /* SETTLE acts on 3 atoms, but the topology format only specifies
1513          * the first atom (for historical reasons).
1514          */
1515         nral_fmt = 1;
1516     }
1517     else
1518     {
1519         nral_fmt = nral;
1520     }
1521
1522     nread = sscanf(line,aaformat[nral_fmt-1],
1523                    &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1524
1525     if (ftype == F_SETTLE)
1526     {
1527         aa[3] = aa[1];
1528         aa[1] = aa[0] + 1;
1529         aa[2] = aa[0] + 2;
1530     }
1531
1532     if (nread < nral_fmt)
1533     {
1534         too_few(wi);
1535         return;
1536     }
1537     else if (nread > nral_fmt)
1538     {
1539         /* this is a hack to allow for virtual sites with swapped parity */
1540         bSwapParity = (aa[nral]<0);
1541         if (bSwapParity)
1542         {
1543             aa[nral] = -aa[nral];
1544         }
1545         ftype = ifunc_index(d,aa[nral]);
1546         if (bSwapParity)
1547         {
1548             switch(ftype)
1549             {
1550             case F_VSITE3FAD:
1551             case F_VSITE3OUT:
1552                 break;
1553             default:
1554                 gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1555                           interaction_function[F_VSITE3FAD].longname,
1556                           interaction_function[F_VSITE3OUT].longname);
1557             }
1558         }
1559     }
1560
1561
1562   /* Check for double atoms and atoms out of bounds */
1563   for(i=0; (i<nral); i++) {
1564     if ( aa[i] < 1 || aa[i] > at->nr )
1565       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1566                 "Atom index (%d) in %s out of bounds (1-%d).\n"
1567                 "This probably means that you have inserted topology section \"%s\"\n"
1568                 "in a part belonging to a different molecule than you intended to.\n"
1569                 "In that case move the \"%s\" section to the right molecule.",
1570                 get_warning_file(wi),get_warning_line(wi),
1571                 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1572     for(j=i+1; (j<nral); j++)
1573       if (aa[i] == aa[j]) {
1574         sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1575         warning(wi,errbuf);
1576       }
1577   }
1578
1579   /* default force parameters  */
1580   for(j=0; (j<MAXATOMLIST); j++)
1581     param.a[j] = aa[j]-1;
1582   for(j=0; (j<MAXFORCEPARAM); j++)
1583     param.c[j] = 0.0;
1584
1585   /* Get force params for normal and free energy perturbation
1586    * studies, as determined by types!
1587    */
1588
1589   if (bBonded) {
1590     bFoundA = default_params(ftype,bondtype,at,atype,&param,FALSE,&param_defA,&nparam_defA);
1591     if (bFoundA) {
1592       /* Copy the A-state and B-state default parameters */
1593       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1594         param.c[j] = param_defA->c[j];
1595     }
1596     bFoundB = default_params(ftype,bondtype,at,atype,&param,TRUE,&param_defB,&nparam_defB);
1597     if (bFoundB) {
1598       /* Copy only the B-state default parameters */
1599       for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1600         param.c[j] = param_defB->c[j];
1601     }
1602   } else if (ftype == F_LJ14) {
1603     bFoundA = default_nb_params(ftype, bondtype,at,&param,0,FALSE,bGenPairs);
1604     bFoundB = default_nb_params(ftype, bondtype,at,&param,0,TRUE, bGenPairs);
1605   } else if (ftype == F_LJC14_Q) {
1606     param.c[0] = fudgeQQ;
1607     /* Fill in the A-state charges as default parameters */
1608     param.c[1] = at->atom[param.a[0]].q;
1609     param.c[2] = at->atom[param.a[1]].q;
1610     /* The default LJ parameters are the standard 1-4 parameters */
1611     bFoundA = default_nb_params(F_LJ14,bondtype,at,&param,3,FALSE,bGenPairs);
1612     bFoundB = TRUE;
1613   } else if (ftype == F_LJC_PAIRS_NB) {
1614     /* Defaults are not supported here */
1615     bFoundA = FALSE;
1616     bFoundB = TRUE;
1617   } else {
1618     gmx_incons("Unknown function type in push_bond");
1619   }
1620
1621   if (nread > nral_fmt) {
1622           /* Manually specified parameters - in this case we discard multiple torsion info! */
1623
1624     strcpy(format,asformat[nral_fmt-1]);
1625     strcat(format,ccformat);
1626
1627     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1628                    &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1629
1630     if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1631       {
1632         /* We only have to issue a warning if these atoms are perturbed! */
1633         bPert = FALSE;
1634         for(j=0; (j<nral); j++)
1635           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1636
1637         if (bPert && *bWarn_copy_A_B)
1638           {
1639             sprintf(errbuf,
1640                     "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1641             warning(wi,errbuf);
1642             *bWarn_copy_A_B = FALSE;
1643           }
1644
1645         /* If only the A parameters were specified, copy them to the B state */
1646         /* The B-state parameters correspond to the first nrfpB
1647          * A-state parameters.
1648          */
1649         for(j=0; (j<NRFPB(ftype)); j++)
1650           cc[nread++] = cc[j];
1651       }
1652
1653     /* If nread was 0 or EOF, no parameters were read => use defaults.
1654      * If nread was nrfpA we copied above so nread=nrfp.
1655      * If nread was nrfp we are cool.
1656      * For F_LJC14_Q we allow supplying fudgeQQ only.
1657      * Anything else is an error!
1658      */
1659     if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1660         !(ftype == F_LJC14_Q && nread == 1))
1661       {
1662         gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1663                   nread,NRFPA(ftype),NRFP(ftype),
1664                   interaction_function[ftype].longname);
1665       }
1666
1667     for(j=0; (j<nread); j++)
1668       param.c[j]=cc[j];
1669
1670     /* Check whether we have to use the defaults */
1671     if (nread == NRFP(ftype))
1672       bDef = FALSE;
1673   }
1674   else
1675   {
1676     nread = 0;
1677   }
1678   /* nread now holds the number of force parameters read! */
1679
1680         if (bDef)
1681         {
1682           /* Use defaults */
1683                 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1684                 if(ftype==F_PDIHS)
1685                 {
1686                         if ((nparam_defA != nparam_defB) || ((nparam_defA>1 || nparam_defB>1) && (param_defA!=param_defB)))
1687                         {
1688                                 sprintf(errbuf,
1689                                                 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1690                                                 "Please specify perturbed parameters manually for this torsion in your topology!");
1691                                 warning_error(wi,errbuf);
1692                         }
1693                 }
1694
1695           if (nread > 0 && nread < NRFPA(ftype))
1696           {
1697                   /* Issue an error, do not use defaults */
1698                   sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1699                   warning_error(wi,errbuf);
1700           }
1701
1702       if (nread == 0 || nread == EOF)
1703           {
1704                   if (!bFoundA)
1705                   {
1706                           if (interaction_function[ftype].flags & IF_VSITE)
1707                           {
1708                                   /* set them to NOTSET, will be calculated later */
1709                                   for(j=0; (j<MAXFORCEPARAM); j++)
1710                                           param.c[j] = NOTSET;
1711
1712                                   if (bSwapParity)
1713                                           param.C1 = -1; /* flag to swap parity of vsite construction */
1714                           }
1715                           else
1716                           {
1717                                   if (bZero)
1718                                   {
1719                                           fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1720                                                           interaction_function[ftype].longname);
1721                                   }
1722                                   else
1723                                   {
1724                                           sprintf(errbuf,"No default %s types",interaction_function[ftype].longname);
1725                                           warning_error(wi,errbuf);
1726                                   }
1727                           }
1728                   }
1729                   else
1730                   {
1731                           if (bSwapParity)
1732                           {
1733                                   switch(ftype)
1734                                   {
1735                                           case F_VSITE3FAD:
1736                                                   param.C0 = 360 - param.C0;
1737                                                   break;
1738                                           case F_VSITE3OUT:
1739                                                   param.C2 = -param.C2;
1740                                                   break;
1741                                   }
1742                           }
1743                   }
1744                   if (!bFoundB)
1745                   {
1746                           /* We only have to issue a warning if these atoms are perturbed! */
1747                           bPert = FALSE;
1748                           for(j=0; (j<nral); j++)
1749                                   bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1750
1751                           if (bPert)
1752                           {
1753                                   sprintf(errbuf,"No default %s types for perturbed atoms, "
1754                                                   "using normal values",interaction_function[ftype].longname);
1755                                   warning(wi,errbuf);
1756                           }
1757                   }
1758           }
1759   }
1760
1761         if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1762                 && param.c[5]!=param.c[2])
1763                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1764                                   "             %s multiplicity can not be perturbed %f!=%f",
1765                                   get_warning_file(wi),get_warning_line(wi),
1766                                   interaction_function[ftype].longname,
1767                                   param.c[2],param.c[5]);
1768
1769         if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1770                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1771                                   "             %s table number can not be perturbed %d!=%d",
1772                                   get_warning_file(wi),get_warning_line(wi),
1773                                   interaction_function[ftype].longname,
1774                                   (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1775
1776         /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1777         if (ftype==F_RBDIHS) {
1778                 nr=0;
1779                 for(i=0;i<NRFP(ftype);i++) {
1780                         if(param.c[i]!=0)
1781                                 nr++;
1782                 }
1783                 if(nr==0)
1784                         return;
1785         }
1786
1787         /* Put the values in the appropriate arrays */
1788         add_param_to_list (&bond[ftype],&param);
1789
1790         /* Push additional torsions from FF for ftype==9 if we have them.
1791          * We have already checked that the A/B states do not differ in this case,
1792          * so we do not have to double-check that again, or the vsite stuff.
1793          * In addition, those torsions cannot be automatically perturbed.
1794          */
1795         if(bDef && ftype==F_PDIHS)
1796         {
1797                 for(i=1;i<nparam_defA;i++)
1798                 {
1799                         /* Advance pointer! */
1800                         param_defA+=2;
1801                         for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1802                                 param.c[j] = param_defA->c[j];
1803                         /* And push the next term for this torsion */
1804                         add_param_to_list (&bond[ftype],&param);
1805                 }
1806         }
1807 }
1808
1809 void push_cmap(directive d, t_params bondtype[], t_params bond[],
1810                            t_atoms *at, gpp_atomtype_t atype, char *line,
1811                            gmx_bool *bWarn_copy_A_B,
1812                warninp_t wi)
1813 {
1814         const char *aaformat[MAXATOMLIST+1]=
1815         {
1816                 "%d",
1817                 "%d%d",
1818                 "%d%d%d",
1819                 "%d%d%d%d",
1820                 "%d%d%d%d%d",
1821                 "%d%d%d%d%d%d",
1822                 "%d%d%d%d%d%d%d"
1823         };
1824
1825         int i,j,nr,ftype,nral,nread,ncmap_params;
1826         int cmap_type;
1827         int aa[MAXATOMLIST+1];
1828         char  errbuf[256];
1829         gmx_bool bFound;
1830         t_param param,paramB,*param_defA,*param_defB;
1831
1832         ftype        = ifunc_index(d,1);
1833         nral         = NRAL(ftype);
1834         nr           = bondtype[ftype].nr;
1835         ncmap_params = 0;
1836
1837         nread = sscanf(line,aaformat[nral-1],
1838                                    &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1839
1840         if (nread < nral)
1841         {
1842                 too_few(wi);
1843                 return;
1844         }
1845         else if (nread == nral)
1846         {
1847                 ftype = ifunc_index(d,1);
1848         }
1849
1850         /* Check for double atoms and atoms out of bounds */
1851         for(i=0;i<nral;i++)
1852         {
1853                 if(aa[i] < 1 || aa[i] > at->nr)
1854                 {
1855                         gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1856                                           "Atom index (%d) in %s out of bounds (1-%d).\n"
1857                                           "This probably means that you have inserted topology section \"%s\"\n"
1858                                           "in a part belonging to a different molecule than you intended to.\n"
1859                                           "In that case move the \"%s\" section to the right molecule.",
1860                                           get_warning_file(wi),get_warning_line(wi),
1861                                           aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1862                 }
1863
1864                 for(j=i+1; (j<nral); j++)
1865                 {
1866                         if (aa[i] == aa[j])
1867                         {
1868                                 sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1869                                 warning(wi,errbuf);
1870                         }
1871                 }
1872         }
1873
1874         /* default force parameters  */
1875         for(j=0; (j<MAXATOMLIST); j++)
1876         {
1877                 param.a[j] = aa[j]-1;
1878         }
1879         for(j=0; (j<MAXFORCEPARAM); j++)
1880         {
1881                 param.c[j] = 0.0;
1882         }
1883
1884         /* Get the cmap type for this cmap angle */
1885         bFound = default_cmap_params(ftype,bondtype,at,atype,&param,FALSE,&cmap_type,&ncmap_params);
1886
1887         /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
1888         if(bFound && ncmap_params==1)
1889         {
1890                 /* Put the values in the appropriate arrays */
1891                 param.c[0]=cmap_type;
1892                 add_param_to_list(&bond[ftype],&param);
1893         }
1894         else
1895         {
1896                 /* This is essentially the same check as in default_cmap_params() done one more time */
1897                 gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
1898                                   param.a[0]+1,param.a[1]+1,param.a[2]+1,param.a[3]+1,param.a[4]+1);
1899         }
1900 }
1901
1902
1903
1904 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1905                   t_atoms *at,gpp_atomtype_t atype,char *line,
1906                   warninp_t wi)
1907 {
1908   char *ptr;
1909   int  type,ftype,j,n,ret,nj,a;
1910   int  *atc=NULL;
1911   double *weight=NULL,weight_tot;
1912   t_param param;
1913
1914   /* default force parameters  */
1915   for(j=0; (j<MAXATOMLIST); j++)
1916     param.a[j] = NOTSET;
1917   for(j=0; (j<MAXFORCEPARAM); j++)
1918     param.c[j] = 0.0;
1919
1920   ptr = line;
1921   ret = sscanf(ptr,"%d%n",&a,&n);
1922   ptr += n;
1923   if (ret == 0)
1924     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1925               "             Expected an atom index in section \"%s\"",
1926               get_warning_file(wi),get_warning_line(wi),
1927               dir2str(d));
1928
1929   param.a[0] = a - 1;
1930
1931   ret = sscanf(ptr,"%d%n",&type,&n);
1932   ptr += n;
1933   ftype = ifunc_index(d,type);
1934
1935   weight_tot = 0;
1936   nj = 0;
1937   do {
1938     ret = sscanf(ptr,"%d%n",&a,&n);
1939     ptr += n;
1940     if (ret > 0) {
1941       if (nj % 20 == 0) {
1942         srenew(atc,nj+20);
1943         srenew(weight,nj+20);
1944       }
1945       atc[nj] = a - 1;
1946       switch (type) {
1947       case 1:
1948         weight[nj] = 1;
1949         break;
1950       case 2:
1951         /* Here we use the A-state mass as a parameter.
1952          * Note that the B-state mass has no influence.
1953          */
1954         weight[nj] = at->atom[atc[nj]].m;
1955         break;
1956       case 3:
1957         weight[nj] = -1;
1958         ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1959         ptr += n;
1960         if (weight[nj] < 0)
1961           gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1962                 "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1963                 get_warning_file(wi),get_warning_line(wi),
1964                 nj+1,atc[nj]+1);
1965         break;
1966       default:
1967         gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1968       }
1969       weight_tot += weight[nj];
1970       nj++;
1971     }
1972   } while (ret > 0);
1973
1974   if (nj == 0)
1975     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1976               "             Expected more than one atom index in section \"%s\"",
1977               get_warning_file(wi),get_warning_line(wi),
1978               dir2str(d));
1979
1980   if (weight_tot == 0)
1981      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1982                "             The total mass of the construting atoms is zero",
1983                get_warning_file(wi),get_warning_line(wi));
1984
1985   for(j=0; j<nj; j++) {
1986     param.a[1] = atc[j];
1987     param.c[0] = nj;
1988     param.c[1] = weight[j]/weight_tot;
1989     /* Put the values in the appropriate arrays */
1990     add_param_to_list (&bond[ftype],&param);
1991   }
1992
1993   sfree(atc);
1994   sfree(weight);
1995 }
1996
1997 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1998               int *nrcopies,
1999               warninp_t wi)
2000 {
2001   int  i,copies;
2002   char type[STRLEN];
2003
2004   *nrcopies=0;
2005   if (sscanf(pline,"%s%d",type,&copies) != 2) {
2006     too_few(wi);
2007     return;
2008   }
2009
2010   /* search moleculename */
2011   for (i=0; ((i<nrmols) && gmx_strcasecmp(type,*(mols[i].name))); i++)
2012     ;
2013
2014   if (i<nrmols) {
2015     *nrcopies        = copies;
2016     *whichmol        = i;
2017   } else
2018     gmx_fatal(FARGS,"No such moleculetype %s",type);
2019 }
2020
2021 void init_block2(t_block2 *b2, int natom)
2022 {
2023   int i;
2024
2025   b2->nr=natom;
2026   snew(b2->nra,b2->nr);
2027   snew(b2->a,b2->nr);
2028   for(i=0; (i<b2->nr); i++)
2029     b2->a[i]=NULL;
2030 }
2031
2032 void done_block2(t_block2 *b2)
2033 {
2034   int i;
2035
2036   if (b2->nr) {
2037     for(i=0; (i<b2->nr); i++)
2038       sfree(b2->a[i]);
2039     sfree(b2->a);
2040     sfree(b2->nra);
2041     b2->nr=0;
2042   }
2043 }
2044
2045 void push_excl(char *line, t_block2 *b2)
2046 {
2047   int  i,j;
2048   int  n;
2049   char base[STRLEN],format[STRLEN];
2050
2051   if (sscanf(line,"%d",&i)==0)
2052     return;
2053
2054   if ((1 <= i) && (i <= b2->nr))
2055     i--;
2056   else {
2057     if (debug)
2058       fprintf(debug,"Unbound atom %d\n",i-1);
2059     return;
2060   }
2061   strcpy(base,"%*d");
2062   do {
2063     strcpy(format,base);
2064     strcat(format,"%d");
2065     n=sscanf(line,format,&j);
2066     if (n == 1) {
2067       if ((1 <= j) && (j <= b2->nr)) {
2068         j--;
2069         srenew(b2->a[i],++(b2->nra[i]));
2070         b2->a[i][b2->nra[i]-1]=j;
2071         /* also add the reverse exclusion! */
2072         srenew(b2->a[j],++(b2->nra[j]));
2073         b2->a[j][b2->nra[j]-1]=i;
2074         strcat(base,"%*d");
2075       }
2076       else
2077         gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
2078     }
2079   } while (n == 1);
2080 }
2081
2082 void b_to_b2(t_blocka *b, t_block2 *b2)
2083 {
2084   int     i;
2085   atom_id j,a;
2086
2087   for(i=0; (i<b->nr); i++)
2088     for(j=b->index[i]; (j<b->index[i+1]); j++) {
2089       a=b->a[j];
2090       srenew(b2->a[i],++b2->nra[i]);
2091       b2->a[i][b2->nra[i]-1]=a;
2092     }
2093 }
2094
2095 void b2_to_b(t_block2 *b2, t_blocka *b)
2096 {
2097   int     i,nra;
2098   atom_id j;
2099
2100   nra=0;
2101   for(i=0; (i<b2->nr); i++) {
2102     b->index[i]=nra;
2103     for(j=0; (j<b2->nra[i]); j++)
2104       b->a[nra+j]=b2->a[i][j];
2105     nra+=b2->nra[i];
2106   }
2107   /* terminate list */
2108   b->index[i]=nra;
2109 }
2110
2111 static int icomp(const void *v1, const void *v2)
2112 {
2113   return (*((atom_id *) v1))-(*((atom_id *) v2));
2114 }
2115
2116 void merge_excl(t_blocka *excl, t_block2 *b2)
2117 {
2118   int     i,k;
2119   atom_id j;
2120   int     nra;
2121
2122   if (!b2->nr)
2123     return;
2124   else if (b2->nr != excl->nr) {
2125     gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2126                 b2->nr,excl->nr);
2127   }
2128   else if (debug)
2129     fprintf(debug,"Entering merge_excl\n");
2130
2131   /* First copy all entries from excl to b2 */
2132   b_to_b2(excl,b2);
2133
2134   /* Count and sort the exclusions */
2135   nra=0;
2136   for(i=0; (i<b2->nr); i++) {
2137     if (b2->nra[i] > 0) {
2138       /* remove double entries */
2139       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
2140       k=1;
2141       for(j=1; (j<b2->nra[i]); j++)
2142         if (b2->a[i][j]!=b2->a[i][k-1]) {
2143           b2->a[i][k]=b2->a[i][j];
2144           k++;
2145         }
2146       b2->nra[i]=k;
2147       nra+=b2->nra[i];
2148     }
2149   }
2150   excl->nra=nra;
2151   srenew(excl->a,excl->nra);
2152
2153   b2_to_b(b2,excl);
2154 }
2155
2156 int add_atomtype_decoupled(t_symtab *symtab,gpp_atomtype_t at,
2157                            t_nbparam ***nbparam,t_nbparam ***pair)
2158 {
2159   t_atom  atom;
2160   t_param param;
2161   int     i,nr;
2162
2163   /* Define an atom type with all parameters set to zero (no interactions) */
2164   atom.q = 0.0;
2165   atom.m = 0.0;
2166   /* Type for decoupled atoms could be anything,
2167    * this should be changed automatically later when required.
2168    */
2169   atom.ptype = eptAtom;
2170   for (i=0; (i<MAXFORCEPARAM); i++)
2171     param.c[i] = 0.0;
2172
2173   nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0,0,0);
2174
2175   /* Add space in the non-bonded parameters matrix */
2176   realloc_nb_params(at,nbparam,pair);
2177
2178   return nr;
2179 }
2180
2181 static void convert_pairs_to_pairsQ(t_params *plist,
2182                                     real fudgeQQ,t_atoms *atoms)
2183 {
2184   t_param *param;
2185   int  i;
2186   real v,w;
2187
2188   /* Copy the pair list to the pairQ list */
2189   plist[F_LJC14_Q] = plist[F_LJ14];
2190   /* Empty the pair list */
2191   plist[F_LJ14].nr    = 0;
2192   plist[F_LJ14].param = NULL;
2193   param = plist[F_LJC14_Q].param;
2194   for(i=0; i<plist[F_LJC14_Q].nr; i++) {
2195     v = param[i].c[0];
2196     w = param[i].c[1];
2197     param[i].c[0] = fudgeQQ;
2198     param[i].c[1] = atoms->atom[param[i].a[0]].q;
2199     param[i].c[2] = atoms->atom[param[i].a[1]].q;
2200     param[i].c[3] = v;
2201     param[i].c[4] = w;
2202   }
2203 }
2204
2205 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
2206 {
2207   int n,ntype,i,j,k;
2208   t_atom *atom;
2209   t_blocka *excl;
2210   gmx_bool bExcl;
2211   t_param param;
2212
2213   n = mol->atoms.nr;
2214   atom = mol->atoms.atom;
2215
2216   ntype = sqrt(nbp->nr);
2217
2218   for (i=0; i<MAXATOMLIST; i++)
2219     param.a[i] = NOTSET;
2220   for (i=0; i<MAXFORCEPARAM; i++)
2221     param.c[i] = NOTSET;
2222
2223   /* Add a pair interaction for all non-excluded atom pairs */
2224   excl = &mol->excls;
2225   for(i=0; i<n; i++) {
2226     for(j=i+1; j<n; j++) {
2227       bExcl = FALSE;
2228       for(k=excl->index[i]; k<excl->index[i+1]; k++) {
2229         if (excl->a[k] == j)
2230           bExcl = TRUE;
2231       }
2232       if (!bExcl) {
2233         if (nb_funct != F_LJ)
2234           gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2235         param.a[0] = i;
2236         param.a[1] = j;
2237         param.c[0] = atom[i].q;
2238         param.c[1] = atom[j].q;
2239         param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2240         param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2241         add_param_to_list(&mol->plist[F_LJC_PAIRS_NB],&param);
2242       }
2243     }
2244   }
2245 }
2246
2247 static void set_excl_all(t_blocka *excl)
2248 {
2249   int nat,i,j,k;
2250
2251   /* Get rid of the current exclusions and exclude all atom pairs */
2252   nat = excl->nr;
2253   excl->nra = nat*nat;
2254   srenew(excl->a,excl->nra);
2255   k = 0;
2256   for(i=0; i<nat; i++) {
2257     excl->index[i] = k;
2258     for(j=0; j<nat; j++)
2259       excl->a[k++] = j;
2260   }
2261   excl->index[nat] = k;
2262 }
2263
2264 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
2265                            int couple_lam0,int couple_lam1)
2266 {
2267   int i;
2268
2269   for(i=0; i<atoms->nr; i++) {
2270     if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW) {
2271       atoms->atom[i].q     = 0.0;
2272     }
2273     if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ) {
2274       atoms->atom[i].type  = atomtype_decouple;
2275     }
2276     if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW) {
2277       atoms->atom[i].qB    = 0.0;
2278     }
2279     if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ) {
2280       atoms->atom[i].typeB = atomtype_decouple;
2281     }
2282   }
2283 }
2284
2285 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
2286                             int couple_lam0,int couple_lam1,
2287                             gmx_bool bCoupleIntra,int nb_funct,t_params *nbp)
2288 {
2289   convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
2290   if (!bCoupleIntra) {
2291     generate_LJCpairsNB(mol,nb_funct,nbp);
2292     set_excl_all(&mol->excls);
2293   }
2294   decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
2295 }