src/kernel/toppush.c

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