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