src/kernel/toppush.c

   1 /*
   2  * $Id$
   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 "gpp_atomtype.h"
  55 #include "gpp_bond_atomtype.h"
  56
  57 void generate_nbparams(int comb,int ftype,t_params *plist,t_atomtype atype)
  58 {
  59   int   i,j,k=-1,nf;
  60   int   nr,nrfp;
  61   real  c,bi,bj,ci,cj,ci0,ci1,ci2,cj0,cj1,cj2;
  62   char  errbuf[256];
  63
  64   /* Lean mean shortcuts */
  65   nr   = get_atomtype_ntypes(atype);
  66   nrfp = NRFP(ftype);
  67   snew(plist->param,nr*nr);
  68   plist->nr=nr*nr;
  69
  70   /* Fill the matrix with force parameters */
  71   switch(ftype) {
  72   case F_LJ:
  73     switch (comb) {
  74     case eCOMB_GEOMETRIC:
  75       /* Gromos rules */
  76       for(i=k=0; (i<nr); i++) {
  77         for(j=0; (j<nr); j++,k++) {
  78           for(nf=0; (nf<nrfp); nf++) {
  79             ci = get_atomtype_nbparam(i,nf,atype);
  80             cj = get_atomtype_nbparam(j,nf,atype);
  81             c = sqrt(ci * cj);
  82             plist->param[k].c[nf]      = c;
  83           }
  84         }
  85       }
  86       break;
  87
  88     case eCOMB_ARITHMETIC:
  89       /* c0 and c1 are epsilon and sigma */
  90       for (i=k=0; (i < nr); i++)
  91         for (j=0; (j < nr); j++,k++) {
  92           ci0 = get_atomtype_nbparam(i,0,atype);
  93           cj0 = get_atomtype_nbparam(j,0,atype);
  94           ci1 = get_atomtype_nbparam(i,1,atype);
  95           cj1 = get_atomtype_nbparam(j,1,atype);
  96           plist->param[k].c[0] = (ci0+cj0)*0.5;
  97           plist->param[k].c[1] = sqrt(ci1*cj1);
  98         }
  99
 100       break;
 101     case eCOMB_GEOM_SIG_EPS:
 102       /* c0 and c1 are epsilon and sigma */
 103       for (i=k=0; (i < nr); i++)
 104         for (j=0; (j < nr); j++,k++) {
 105           ci0 = get_atomtype_nbparam(i,0,atype);
 106           cj0 = get_atomtype_nbparam(j,0,atype);
 107           ci1 = get_atomtype_nbparam(i,1,atype);
 108           cj1 = get_atomtype_nbparam(j,1,atype);
 109           plist->param[k].c[0] = sqrt(ci0*cj0);
 110           plist->param[k].c[1] = sqrt(ci1*cj1);
 111         }
 112
 113       break;
 114     default:
 115       gmx_fatal(FARGS,"No such combination rule %d",comb);
 116     }
 117     if (plist->nr != k)
 118       gmx_incons("Topology processing, generate nb parameters");
 119     break;
 120
 121   case F_BHAM:
 122     /* Buckingham rules */
 123     for(i=k=0; (i<nr); i++)
 124       for(j=0; (j<nr); j++,k++) {
 125         ci0 = get_atomtype_nbparam(i,0,atype);
 126         cj0 = get_atomtype_nbparam(j,0,atype);
 127         ci2 = get_atomtype_nbparam(i,2,atype);
 128         cj2 = get_atomtype_nbparam(j,2,atype);
 129         bi = get_atomtype_nbparam(i,1,atype);
 130         bj = get_atomtype_nbparam(j,1,atype);
 131         plist->param[k].c[0] = sqrt(ci0 * cj0);
 132         if ((bi == 0) || (bj == 0))
 133           plist->param[k].c[1] = 0;
 134         else
 135           plist->param[k].c[1] = 2.0/(1/bi+1/bj);
 136         plist->param[k].c[2] = sqrt(ci2 * cj2);
 137       }
 138
 139     break;
 140   default:
 141     sprintf(errbuf,"Invalid nonbonded type %s",
 142             interaction_function[ftype].longname);
 143     warning_error(errbuf);
 144   }
 145 }
 146
 147 static void realloc_nb_params(t_atomtype at,
 148                               t_nbparam ***nbparam,t_nbparam ***pair)
 149 {
 150   /* Add space in the non-bonded parameters matrix */
 151   int atnr = get_atomtype_ntypes(at);
 152   srenew(*nbparam,atnr);
 153   snew((*nbparam)[atnr-1],atnr);
 154   if (pair) {
 155     srenew(*pair,atnr);
 156     snew((*pair)[atnr-1],atnr);
 157   }
 158 }
 159
 160 static void copy_B_from_A(int ftype,double *c)
 161 {
 162   int nrfpA,nrfpB,i;
 163
 164   nrfpA = NRFPA(ftype);
 165   nrfpB = NRFPB(ftype);
 166
 167   /* Copy the B-state from the A-state */
 168   if (interaction_function[ftype].flags & IF_TABULATED) {
 169     /* The only case where the first B parameter does not correspond
 170      * to the first A parameter.
 171      */
 172     c[2] = c[1];
 173   } else {
 174     /* Copy the B parameters from the first nrfpB A parameters */
 175     for(i=0; (i<nrfpB); i++) {
 176       c[nrfpA+i] = c[i];
 177     }
 178   }
 179 }
 180
 181 void push_at (t_symtab *symtab, t_atomtype at, t_bond_atomtype bat,
 182               char *line,int nb_funct,
 183               t_nbparam ***nbparam,t_nbparam ***pair)
 184 {
 185   typedef struct {
 186     char *entry;
 187     int  ptype;
 188   } t_xlate;
 189   t_xlate xl[eptNR] = {
 190     { "A",   eptAtom },
 191     { "N",   eptNucleus },
 192     { "S",   eptShell },
 193     { "B",   eptBond },
 194     { "V",   eptVSite },
 195   };
 196
 197   int    nr,i,nfields,j,pt,nfp0=-1;
 198   int    batype_nr,nread;
 199   char   type[STRLEN],btype[STRLEN],ptype[STRLEN];
 200   double m,q;
 201   double c[MAXFORCEPARAM];
 202   double radius,vol,surftens;
 203   char   tmpfield[12][100]; /* Max 12 fields of width 100 */
 204   char   errbuf[256];
 205   t_atom  *atom;
 206   t_param *param;
 207   int    atomnr;
 208   bool   have_atomic_number;
 209   bool   have_bonded_type;
 210
 211   snew(atom,1);
 212   snew(param,1);
 213
 214   /* First assign input line to temporary array */
 215   nfields=sscanf(line,"%s%s%s%s%s%s%s%s%s%s%s%s",
 216                  tmpfield[0],tmpfield[1],tmpfield[2],tmpfield[3],tmpfield[4],tmpfield[5],
 217                  tmpfield[6],tmpfield[7],tmpfield[8],tmpfield[9],tmpfield[10],tmpfield[11]);
 218
 219   /* Comments on optional fields in the atomtypes section:
 220    *
 221    * The force field format is getting a bit old. For OPLS-AA we needed
 222    * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
 223    * we also needed the atomic numbers.
 224    * To avoid making all old or user-generated force fields unusable we
 225    * have introduced both these quantities as optional columns, and do some
 226    * acrobatics to check whether they are present or not.
 227    * This will all look much nicer when we switch to XML... sigh.
 228    *
 229    * Field 0 (mandatory) is the nonbonded type name. (string)
 230    * Field 1 (optional)  is the bonded type (string)
 231    * Field 2 (optional)  is the atomic number (int)
 232    * Field 3 (mandatory) is the mass (numerical)
 233    * Field 4 (mandatory) is the charge (numerical)
 234    * Field 5 (mandatory) is the particle type (single character)
 235    * This is followed by a number of nonbonded parameters.
 236    *
 237    * The safest way to identify the format is the particle type field.
 238    *
 239    * So, here is what we do:
 240    *
 241    * A. Read in the first six fields as strings
 242    * B. If field 3 (starting from 0) is a single char, we have neither
 243    *    bonded_type or atomic numbers.
 244    * C. If field 5 is a single char we have both.
 245    * D. If field 4 is a single char we check field 1. If this begins with
 246    *    an alphabetical character we have bonded types, otherwise atomic numbers.
 247    */
 248
 249   if(nfields<6)
 250   {
 251       too_few();
 252       return;
 253   }
 254
 255   if( (strlen(tmpfield[5])==1) && isalpha(tmpfield[5][0]) )
 256   {
 257       have_bonded_type   = TRUE;
 258       have_atomic_number = TRUE;
 259   }
 260   else if( (strlen(tmpfield[3])==1) && isalpha(tmpfield[3][0]) )
 261   {
 262       have_bonded_type   = FALSE;
 263       have_atomic_number = FALSE;
 264   }
 265   else
 266   {
 267       have_bonded_type   = ( isalpha(tmpfield[1][0]) != 0 );
 268       have_atomic_number = !have_bonded_type;
 269   }
 270
 271   /* optional fields */
 272   surftens = -1;
 273   vol      =  0;
 274   radius   =  0;
 275   atomnr   = -1;
 276
 277   switch (nb_funct) {
 278
 279   case F_LJ:
 280     nfp0 = 2;
 281
 282     if( have_atomic_number )
 283     {
 284         if ( have_bonded_type )
 285         {
 286             nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf",
 287                            type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],
 288                            &radius,&vol,&surftens);
 289             if(nread < 8)
 290             {
 291                 too_few();
 292                 return;
 293             }
 294         }
 295         else
 296         {
 297             /* have_atomic_number && !have_bonded_type */
 298             nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf",
 299                            type,&atomnr,&m,&q,ptype,&c[0],&c[1],
 300                            &radius,&vol,&surftens);
 301             if(nread < 7)
 302             {
 303                 too_few();
 304                 return;
 305             }
 306         }
 307     }
 308     else
 309     {
 310         if ( have_bonded_type )
 311         {
 312             /* !have_atomic_number && have_bonded_type */
 313             nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf",
 314                            type,btype,&m,&q,ptype,&c[0],&c[1],
 315                            &radius,&vol,&surftens);
 316             if(nread < 7)
 317             {
 318                 too_few();
 319                 return;
 320             }
 321         }
 322         else
 323         {
 324             /* !have_atomic_number && !have_bonded_type */
 325             nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf",
 326                            type,&m,&q,ptype,&c[0],&c[1],
 327                            &radius,&vol,&surftens);
 328             if(nread < 6)
 329             {
 330                 too_few();
 331                 return;
 332             }
 333         }
 334     }
 335
 336     if( !have_bonded_type )
 337     {
 338         strcpy(btype,type);
 339     }
 340
 341     if( !have_atomic_number )
 342     {
 343         atomnr = -1;
 344     }
 345
 346     break;
 347
 348   case F_BHAM:
 349     nfp0 = 3;
 350
 351       if( have_atomic_number )
 352       {
 353           if ( have_bonded_type )
 354           {
 355               nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
 356                              type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
 357                              &radius,&vol,&surftens);
 358               if(nread < 9)
 359               {
 360                   too_few();
 361                   return;
 362               }
 363           }
 364           else
 365           {
 366               /* have_atomic_number && !have_bonded_type */
 367               nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
 368                              type,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
 369                              &radius,&vol,&surftens);
 370               if(nread < 8)
 371               {
 372                   too_few();
 373                   return;
 374               }
 375           }
 376       }
 377       else
 378       {
 379           if ( have_bonded_type )
 380           {
 381               /* !have_atomic_number && have_bonded_type */
 382               nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
 383                              type,btype,&m,&q,ptype,&c[0],&c[1],&c[2],
 384                              &radius,&vol,&surftens);
 385               if(nread < 8)
 386               {
 387                   too_few();
 388                   return;
 389               }
 390           }
 391           else
 392           {
 393               /* !have_atomic_number && !have_bonded_type */
 394               nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
 395                              type,&m,&q,ptype,&c[0],&c[1],&c[2],
 396                              &radius,&vol,&surftens);
 397               if(nread < 7)
 398               {
 399                   too_few();
 400                   return;
 401               }
 402           }
 403       }
 404
 405       if( !have_bonded_type )
 406       {
 407           strcpy(btype,type);
 408       }
 409
 410       if( !have_atomic_number )
 411       {
 412           atomnr = -1;
 413       }
 414
 415       break;
 416
 417       default:
 418           gmx_fatal(FARGS,"Invalid function type %d in push_at %s %d",nb_funct,
 419                     __FILE__,__LINE__);
 420   }
 421   for(j=nfp0; (j<MAXFORCEPARAM); j++)
 422     c[j]=0.0;
 423
 424   if(strlen(type)==1 && isdigit(type[0]))
 425      gmx_fatal(FARGS,"Atom type names can't be single digits.");
 426
 427   if(strlen(btype)==1 && isdigit(btype[0]))
 428      gmx_fatal(FARGS,"Bond atom type names can't be single digits.");
 429
 430   /* Hack to read old topologies */
 431   if (strcasecmp(ptype,"D") == 0)
 432     sprintf(ptype,"V");
 433   for(j=0; (j<eptNR); j++)
 434     if (strcasecmp(ptype,xl[j].entry) == 0)
 435       break;
 436   if (j == eptNR)
 437     gmx_fatal(FARGS,"Invalid particle type %s on line %s",
 438                 ptype,line);
 439   pt=xl[j].ptype;
 440   if (debug)
 441     fprintf(debug,"ptype: %s\n",ptype_str[pt]);
 442
 443   atom->q = q;
 444   atom->m = m;
 445   atom->ptype = pt;
 446   for (i=0; (i<MAXFORCEPARAM); i++)
 447     param->c[i] = c[i];
 448
 449   if ((batype_nr = get_bond_atomtype_type(btype,bat)) == NOTSET)
 450     add_bond_atomtype(bat,symtab,btype);
 451   batype_nr = get_bond_atomtype_type(btype,bat);
 452
 453   if ((nr = get_atomtype_type(type,at)) != NOTSET) {
 454     sprintf(errbuf,"Overriding atomtype %s",type);
 455     warning(errbuf);
 456     if ((nr = set_atomtype(nr,at,symtab,atom,type,param,batype_nr,
 457                            radius,vol,surftens,atomnr)) == NOTSET)
 458       gmx_fatal(FARGS,"Replacing atomtype %s failed",type);
 459   }
 460   else if ((nr = add_atomtype(at,symtab,atom,type,param,
 461                               batype_nr,radius,vol,
 462                               surftens,atomnr)) == NOTSET)
 463     gmx_fatal(FARGS,"Adding atomtype %s failed",type);
 464   else {
 465     /* Add space in the non-bonded parameters matrix */
 466     realloc_nb_params(at,nbparam,pair);
 467   }
 468   sfree(atom);
 469   sfree(param);
 470 }
 471
 472 static void push_bondtype(t_params *bt,t_param *b,int nral,int ftype,
 473                           char *line)
 474 {
 475   int  i,j;
 476   bool bTest,bFound,bId;
 477   int  nr   = bt->nr;
 478   int  nrfp = NRFP(ftype);
 479   char  errbuf[256];
 480
 481   /* Check if this entry overwrites another */
 482   bFound=FALSE;
 483   for (i=0; (i < nr); i++) {
 484     bTest = TRUE;
 485     for (j=0; (j < nral); j++)
 486       bTest=(bTest && (b->a[j] == bt->param[i].a[j]));
 487     if (!bTest) {
 488       bTest=TRUE;
 489       for(j=0; (j<nral); j++)
 490         bTest=(bTest && (b->a[nral-1-j] == bt->param[i].a[j]));
 491     }
 492     if (bTest) {
 493       if (!bFound) {
 494         bId = TRUE;
 495         for (j=0; (j < nrfp); j++)
 496           bId = bId && (bt->param[i].c[j] == b->c[j]);
 497         if (!bId) {
 498           sprintf(errbuf,"Overriding %s parameters,",
 499                   interaction_function[ftype].longname);
 500           warning(errbuf);
 501           fprintf(stderr,"  old:");
 502           for (j=0; (j < nrfp); j++)
 503             fprintf(stderr," %g",bt->param[i].c[j]);
 504           fprintf(stderr," \n  new: %s\n\n",line);
 505         }
 506       }
 507       /* Overwrite it! */
 508       for (j=0; (j < nrfp); j++)
 509         bt->param[i].c[j] = b->c[j];
 510       bFound=TRUE;
 511     }
 512   }
 513   if (!bFound) {
 514     /* alloc */
 515     pr_alloc (2,bt);
 516
 517     /* fill the arrays up and down */
 518     for (j=0; (j < nrfp); j++) {
 519       bt->param[bt->nr].c[j]   = b->c[j];
 520       bt->param[bt->nr+1].c[j] = b->c[j];
 521     }
 522     for (j=0; (j<nral); j++) {
 523       bt->param[bt->nr].a[j]   = b->a[j];
 524       bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
 525     }
 526     bt->nr += 2;
 527   }
 528 }
 529
 530 void push_bt(directive d,t_params bt[],int nral,
 531              t_atomtype at,
 532              t_bond_atomtype bat,char *line)
 533 {
 534   const char *formal[MAXATOMLIST+1] = {
 535     "%s",
 536     "%s%s",
 537     "%s%s%s",
 538     "%s%s%s%s",
 539     "%s%s%s%s%s",
 540     "%s%s%s%s%s%s"
 541   };
 542   const char *formnl[MAXATOMLIST+1] = {
 543     "%*s",
 544     "%*s%*s",
 545     "%*s%*s%*s",
 546     "%*s%*s%*s%*s",
 547     "%*s%*s%*s%*s%*s",
 548     "%*s%*s%*s%*s%*s%*s"
 549   };
 550   const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
 551   int      i,ft,ftype,nn,nrfp,nrfpA,nrfpB;
 552   char     f1[STRLEN];
 553   char     alc[MAXATOMLIST+1][20];
 554   /* One force parameter more, so we can check if we read too many */
 555   double   c[MAXFORCEPARAM+1];
 556   t_param  p;
 557   char  errbuf[256];
 558
 559   if ((bat && at) || (!bat && !at))
 560     gmx_incons("You should pass either bat or at to push_bt");
 561
 562   /* Make format string (nral ints+functype) */
 563   if ((nn=sscanf(line,formal[nral],
 564                  alc[0],alc[1],alc[2],alc[3],alc[4],alc[5])) != nral+1) {
 565     sprintf(errbuf,"Not enough atomtypes (%d instead of %d)",nn-1,nral);
 566     warning_error(errbuf);
 567     return;
 568   }
 569
 570   ft    = atoi(alc[nral]);
 571   ftype = ifunc_index(d,ft);
 572   nrfp  = NRFP(ftype);
 573   nrfpA = interaction_function[ftype].nrfpA;
 574   nrfpB = interaction_function[ftype].nrfpB;
 575   strcpy(f1,formnl[nral]);
 576   strcat(f1,formlf);
 577   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]))
 578       != nrfp) {
 579     if (nn == nrfpA) {
 580       /* Copy the B-state from the A-state */
 581       copy_B_from_A(ftype,c);
 582     } else {
 583       if (nn < nrfpA) {
 584       warning_error("Not enough parameters");
 585       } else if (nn > nrfpA && nn < nrfp) {
 586         warning_error("Too many parameters or not enough parameters for topology B");
 587       } else if (nn > nrfp) {
 588         warning_error("Too many parameters");
 589       }
 590       for(i=nn; (i<nrfp); i++)
 591         c[i] = 0.0;
 592     }
 593   }
 594   for(i=0; (i<nral); i++) {
 595     if (at && ((p.a[i]=get_atomtype_type(alc[i],at)) == NOTSET))
 596       gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
 597     else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
 598       gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
 599   }
 600   for(i=0; (i<nrfp); i++)
 601     p.c[i]=c[i];
 602   push_bondtype (&(bt[ftype]),&p,nral,ftype,line);
 603 }
 604
 605
 606 void push_dihedraltype(directive d,t_params bt[],
 607                        t_bond_atomtype bat,char *line)
 608 {
 609   const char *formal[MAXATOMLIST+1] = {
 610     "%s",
 611     "%s%s",
 612     "%s%s%s",
 613     "%s%s%s%s",
 614     "%s%s%s%s%s",
 615     "%s%s%s%s%s%s"
 616   };
 617   const char *formnl[MAXATOMLIST+1] = {
 618     "%*s",
 619     "%*s%*s",
 620     "%*s%*s%*s",
 621     "%*s%*s%*s%*s",
 622     "%*s%*s%*s%*s%*s",
 623     "%*s%*s%*s%*s%*s%*s"
 624   };
 625   const char *formlf[MAXFORCEPARAM] = {
 626     "%lf",
 627     "%lf%lf",
 628     "%lf%lf%lf",
 629     "%lf%lf%lf%lf",
 630     "%lf%lf%lf%lf%lf",
 631     "%lf%lf%lf%lf%lf%lf",
 632     "%lf%lf%lf%lf%lf%lf%lf",
 633     "%lf%lf%lf%lf%lf%lf%lf%lf",
 634     "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 635     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 636     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 637     "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 638   };
 639   int      i,ft,ftype,nn,nrfp,nrfpA,nrfpB,nral;
 640   char     f1[STRLEN];
 641   char     alc[MAXATOMLIST+1][20];
 642   double   c[MAXFORCEPARAM];
 643   t_param  p;
 644   char  errbuf[256];
 645
 646   /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
 647    *
 648    * We first check for 2 atoms with the 3th column being an integer
 649    * defining the type. If this isn't the case, we try it with 4 atoms
 650    * and the 5th column defining the dihedral type.
 651    */
 652   nn=sscanf(line,formal[4],alc[0],alc[1],alc[2],alc[3],alc[4]);
 653   if(nn>=3 && strlen(alc[2])==1 && isdigit(alc[2][0])) {
 654     nral=2;
 655     ft    = atoi(alc[nral]);
 656     /* Move atom types around a bit and use 'X' for wildcard atoms
 657      * to create a 4-atom dihedral definition with arbitrary atoms in
 658      * position 1 and 4.
 659      */
 660     if(alc[2][0]=='2') {
 661       /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
 662       strcpy(alc[3],alc[1]);
 663       sprintf(alc[2],"X");
 664       sprintf(alc[1],"X");
 665       /* alc[0] stays put */
 666     } else {
 667       /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
 668       sprintf(alc[3],"X");
 669       strcpy(alc[2],alc[1]);
 670       strcpy(alc[1],alc[0]);
 671       sprintf(alc[0],"X");
 672     }
 673   } else if(nn==5 && strlen(alc[4])==1 && isdigit(alc[4][0])) {
 674     nral=4;
 675     ft    = atoi(alc[nral]);
 676   } else {
 677     sprintf(errbuf,"Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)",nn-1);
 678     warning_error(errbuf);
 679     return;
 680   }
 681
 682   ftype = ifunc_index(d,ft);
 683   nrfp  = NRFP(ftype);
 684   nrfpA = interaction_function[ftype].nrfpA;
 685   nrfpB = interaction_function[ftype].nrfpB;
 686
 687   strcpy(f1,formnl[nral]);
 688   strcat(f1,formlf[nrfp-1]);
 689
 690   /* Check number of parameters given */
 691   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]))
 692       != nrfp) {
 693     if (nn == nrfpA) {
 694       /* Copy the B-state from the A-state */
 695       copy_B_from_A(ftype,c);
 696     } else {
 697       if (nn < nrfpA) {
 698       warning_error("Not enough parameters");
 699       } else if (nn > nrfpA && nn < nrfp) {
 700         warning_error("Too many parameters or not enough parameters for topology B");
 701       } else if (nn > nrfp) {
 702         warning_error("Too many parameters");
 703       }
 704       for(i=nn; (i<nrfp); i++)
 705         c[i] = 0.0;
 706     }
 707   }
 708
 709   for(i=0; (i<4); i++) {
 710     if(!strcmp(alc[i],"X"))
 711       p.a[i]=-1;
 712     else {
 713       if ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET)
 714         gmx_fatal(FARGS,"Unknown bond_atomtype %s",alc[i]);
 715     }
 716   }
 717   for(i=0; (i<nrfp); i++)
 718     p.c[i]=c[i];
 719   /* Always use 4 atoms here, since we created two wildcard atoms
 720    * if there wasn't of them 4 already.
 721    */
 722   push_bondtype (&(bt[ftype]),&p,4,ftype,line);
 723 }
 724
 725
 726 void push_nbt(directive d,t_nbparam **nbt,t_atomtype atype,
 727               char *pline,int nb_funct)
 728 {
 729   /* swap the atoms */
 730   const char *form2="%*s%*s%*s%lf%lf";
 731   const char *form3="%*s%*s%*s%lf%lf%lf";
 732   const char *form4="%*s%*s%*s%lf%lf%lf%lf";
 733   char    a0[80],a1[80];
 734   int     i,f,n,ftype,atnr,nrfp;
 735   double  c[4];
 736   real    cr[4],sig6;
 737   atom_id ai,aj;
 738   t_nbparam *nbp;
 739   bool    bId;
 740   char  errbuf[256];
 741
 742   if (sscanf (pline,"%s%s%d",a0,a1,&f) != 3) {
 743     too_few();
 744     return;
 745   }
 746
 747   ftype=ifunc_index(d,f);
 748
 749   if (ftype != nb_funct) {
 750     sprintf(errbuf,"Trying to add %s while the default nonbond type is %s",
 751             interaction_function[ftype].longname,
 752             interaction_function[nb_funct].longname);
 753     warning_error(errbuf);
 754     return;
 755   }
 756
 757   /* Get the force parameters */
 758   nrfp = NRFP(ftype);
 759   if (ftype == F_LJ14) {
 760     n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
 761     if (n < 2) {
 762       too_few();
 763       return;
 764     }
 765     /* When the B topology parameters are not set,
 766      * copy them from topology A
 767      */
 768     for(i=n; i<nrfp; i++)
 769       c[i] = c[i-2];
 770   }
 771   else if (ftype == F_LJC14_Q) {
 772     n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
 773     if (n < 4) {
 774       too_few();
 775       return;
 776     }
 777   }
 778   else if (nrfp == 2) {
 779     if (sscanf(pline,form2,&c[0],&c[1]) != 2) {
 780       too_few();
 781       return;
 782     }
 783   }
 784   else if (nrfp == 3) {
 785     if (sscanf(pline,form3,&c[0],&c[1],&c[2]) != 3) {
 786       too_few();
 787       return;
 788     }
 789   }
 790   else {
 791     gmx_fatal(FARGS,"Number of force parameters for nonbonded interactions is %d"
 792                 " in file %s, line %d",nrfp,__FILE__,__LINE__);
 793   }
 794   for(i=0; (i<nrfp); i++)
 795     cr[i] = c[i];
 796
 797   /* Put the parameters in the matrix */
 798   if ((ai = get_atomtype_type (a0,atype)) == NOTSET)
 799     gmx_fatal(FARGS,"Atomtype %s not found",a0);
 800   if ((aj = get_atomtype_type (a1,atype)) == NOTSET)
 801     gmx_fatal(FARGS,"Atomtype %s not found",a1);
 802   nbp = &(nbt[max(ai,aj)][min(ai,aj)]);
 803
 804   if (nbp->bSet) {
 805     bId = TRUE;
 806     for (i=0; i<nrfp; i++)
 807       bId = bId && (nbp->c[i] == cr[i]);
 808     if (!bId) {
 809       sprintf(errbuf,"Overriding non-bonded parameters,");
 810       warning(errbuf);
 811       fprintf(stderr,"  old:");
 812       for (i=0; i<nrfp; i++)
 813         fprintf(stderr," %g",nbp->c[i]);
 814       fprintf(stderr," new\n%s\n",pline);
 815     }
 816   }
 817   nbp->bSet = TRUE;
 818   for (i=0; i<nrfp; i++)
 819     nbp->c[i] = cr[i];
 820 }
 821
 822 static void push_atom_now(t_symtab *symtab,t_atoms *at,int atomnr,
 823                           int atomicnumber,
 824                           int type,char *ctype,int ptype,
 825                           int resnumber,int cgnumber,
 826                           char *resname,char *name,real m0,real q0,
 827                           int typeB,char *ctypeB,real mB,real qB)
 828 {
 829   int j,resnr_diff=-1;
 830   int nr = at->nr;
 831
 832   if (((nr==0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
 833     gmx_fatal(FARGS,"Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)",atomnr,at->nr);
 834   if (nr)
 835     resnr_diff = resnumber - (at->atom[at->nr-1].resnr+1);
 836   if (((nr==0) && (resnumber != 1)) ||
 837       (nr && (resnr_diff != 0) && (resnr_diff != 1)))
 838     gmx_fatal(FARGS,"Residue numbers in the .top are not numbered consecutively from 1 (rather, resnumber = %d and resnr_diff = %d)",resnumber,resnr_diff);
 839
 840   /* New atom instance
 841    * get new space for arrays
 842    */
 843   srenew(at->atom,nr+1);
 844   srenew(at->atomname,nr+1);
 845   srenew(at->atomtype,nr+1);
 846   srenew(at->atomtypeB,nr+1);
 847
 848   if (resnumber > at->nres) {
 849     at->nres=resnumber;
 850     srenew(at->resname,resnumber);
 851     at->resname[resnumber-1] = put_symtab(symtab,resname);
 852   }
 853   /* fill the list */
 854   at->atom[nr].type  = type;
 855   at->atom[nr].ptype = ptype;
 856   at->atom[nr].q     = q0;
 857   at->atom[nr].m     = m0;
 858   at->atom[nr].typeB = typeB;
 859   at->atom[nr].qB    = qB;
 860   at->atom[nr].mB    = mB;
 861
 862   at->atom[nr].resnr = resnumber-1;
 863   at->atom[nr].atomnumber = atomicnumber;
 864   at->atomname[nr] = put_symtab(symtab,name);
 865   at->atomtype[nr] = put_symtab(symtab,ctype);
 866   at->atomtypeB[nr] = put_symtab(symtab,ctypeB);
 867   at->nr++;
 868 }
 869
 870 void push_cg(t_block *block, int *lastindex, int index, int a)
 871 {
 872   if (debug)
 873     fprintf (debug,"Index %d, Atom %d\n",index,a);
 874
 875   if (((block->nr) && (*lastindex != index)) || (!block->nr)) {
 876     /* add a new block */
 877     block->nr++;
 878     srenew(block->index,block->nr+1);
 879   }
 880   block->index[block->nr] = a + 1;
 881   *lastindex = index;
 882 }
 883
 884 void push_atom(t_symtab *symtab,t_block *cgs,
 885                t_atoms *at,t_atomtype atype,char *line,int *lastcg)
 886 {
 887   int           nr,ptype;
 888   int           resnumber,cgnumber,atomnr,type,typeB,nscan;
 889   char          id[STRLEN],ctype[STRLEN],ctypeB[STRLEN],
 890                 resname[STRLEN],name[STRLEN],check[STRLEN];
 891   double        m,q,mb,qb;
 892   real          m0,q0,mB,qB;
 893
 894   /* Make a shortcut for writing in this molecule  */
 895   nr = at->nr;
 896
 897   /* Fixed parameters */
 898   if (sscanf(line,"%s%s%d%s%s%d",
 899              id,ctype,&resnumber,resname,name,&cgnumber) != 6) {
 900     too_few();
 901     return;
 902   }
 903   sscanf(id,"%d",&atomnr);
 904   if ((type  = get_atomtype_type(ctype,atype)) == NOTSET)
 905     gmx_fatal(FARGS,"Atomtype %s not found",ctype);
 906   ptype = get_atomtype_ptype(type,atype);
 907
 908   /* Set default from type */
 909   q0    = get_atomtype_qA(type,atype);
 910   m0    = get_atomtype_massA(type,atype);
 911   typeB = type;
 912   qB    = q0;
 913   mB    = m0;
 914
 915   /* Optional parameters */
 916   nscan = sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
 917                  &q,&m,ctypeB,&qb,&mb,check);
 918
 919   /* Nasty switch that falls thru all the way down! */
 920   if (nscan > 0) {
 921     q0 = qB = q;
 922     if (nscan > 1) {
 923       m0 = mB = m;
 924       if (nscan > 2) {
 925         typeB=get_atomtype_type(ctypeB,atype);
 926         qB = get_atomtype_qA(typeB,atype);
 927         mB = get_atomtype_massA(typeB,atype);
 928         if (nscan > 3) {
 929           qB = qb;
 930           if (nscan > 4) {
 931             mB = mb;
 932             if (nscan > 5)
 933               warning_error("Too many parameters");
 934           }
 935         }
 936       }
 937     }
 938   }
 939   if (debug)
 940     fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
 941
 942   push_cg(cgs,lastcg,cgnumber,nr);
 943
 944   push_atom_now(symtab,at,atomnr,get_atomtype_atomnumber(type,atype),
 945                 type,ctype,ptype,resnumber,cgnumber,
 946                 resname,name,m0,q0,typeB,
 947                 typeB==type ? ctype : ctypeB,mB,qB);
 948 }
 949
 950 void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line)
 951 {
 952   char type[STRLEN];
 953   int nrexcl,i;
 954   t_molinfo *newmol;
 955
 956   if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
 957     warning_error("Expected a molecule type name and nrexcl");
 958   }
 959
 960   /* Test if this atomtype overwrites another */
 961   i = 0;
 962   while (i < *nmol) {
 963     if (strcasecmp(*((*mol)[i].name),type) == 0)
 964       gmx_fatal(FARGS,"moleculetype %s is redefined",type);
 965     i++;
 966   }
 967
 968   (*nmol)++;
 969   srenew(*mol,*nmol);
 970   newmol = &((*mol)[*nmol-1]);
 971   init_molinfo(newmol);
 972
 973   /* Fill in the values */
 974   newmol->name     = put_symtab(symtab,type);
 975   newmol->nrexcl   = nrexcl;
 976   newmol->excl_set = FALSE;
 977 }
 978
 979 static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
 980                               t_param *p,int c_start,bool bB,bool bGenPairs)
 981 {
 982   int      i,j,ti,tj,ntype;
 983   bool     bFound;
 984   t_param  *pi=NULL;
 985   int      nr   = bt[ftype].nr;
 986   int      nral = NRAL(ftype);
 987   int      nrfp = interaction_function[ftype].nrfpA;
 988   int      nrfpB= interaction_function[ftype].nrfpB;
 989
 990   if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
 991     return TRUE;
 992
 993   bFound=FALSE;
 994   if(bGenPairs) {
 995     /* First test the generated-pair position to save
 996      * time when we have 1000*1000 entries for e.g. OPLS...
 997      */
 998     ntype=sqrt(nr);
 999     if (bB) {
1000       ti=at->atom[p->a[0]].typeB;
1001       tj=at->atom[p->a[1]].typeB;
1002     } else {
1003       ti=at->atom[p->a[0]].type;
1004       tj=at->atom[p->a[1]].type;
1005     }
1006     pi=&(bt[ftype].param[ntype*ti+tj]);
1007     bFound=((ti == pi->a[0]) && (tj == pi->a[1]));
1008   }
1009
1010   /* Search explicitly if we didnt find it */
1011   if(!bFound) {
1012     for (i=0; ((i < nr) && !bFound); i++) {
1013       pi=&(bt[ftype].param[i]);
1014       if (bB)
1015         for (j=0; ((j < nral) &&
1016                    (at->atom[p->a[j]].typeB == pi->a[j])); j++);
1017       else
1018         for (j=0; ((j < nral) &&
1019                    (at->atom[p->a[j]].type == pi->a[j])); j++);
1020       bFound=(j==nral);
1021     }
1022   }
1023
1024   if (bFound) {
1025     if (bB) {
1026       if (nrfp+nrfpB > MAXFORCEPARAM) {
1027         gmx_incons("Too many force parameters");
1028       }
1029       for (j=c_start; (j < nrfpB); j++)
1030         p->c[nrfp+j] = pi->c[j];
1031     }
1032     else
1033       for (j=c_start; (j < nrfp); j++)
1034         p->c[j] = pi->c[j];
1035   }
1036   else {
1037     for (j=c_start; (j < nrfp); j++)
1038       p->c[j] = 0.0;
1039   }
1040   return bFound;
1041 }
1042
1043
1044 static bool default_params(int ftype,t_params bt[],
1045                            t_atoms *at,t_atomtype atype,
1046                            t_param *p,bool bB,
1047                            t_param **param_def)
1048 {
1049   int      i,j;
1050   bool     bFound;
1051   t_param  *pi=NULL;
1052   int      nr   = bt[ftype].nr;
1053   int      nral = NRAL(ftype);
1054   int      nrfpA= interaction_function[ftype].nrfpA;
1055   int      nrfpB= interaction_function[ftype].nrfpB;
1056
1057   if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1058     return TRUE;
1059
1060
1061   /* We allow wildcards now. The first type (with or without wildcards) that
1062    * fits is used, so you should probably put the wildcarded bondtypes
1063    * at the end of each section.
1064    */
1065   bFound=FALSE;
1066   /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1067    * special case for this. Check for B state outside loop to speed it up.
1068    */
1069   if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS) {
1070     if (bB) {
1071       for (i=0; ((i < nr) && !bFound); i++) {
1072         pi=&(bt[ftype].param[i]);
1073         bFound=
1074           (((pi->AI==-1) ||
1075             (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1076            ((pi->AJ==-1) ||
1077             (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1078            ((pi->AK==-1) ||
1079             (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1080            ((pi->AL==-1) ||
1081             (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1082            );
1083       }
1084     } else {
1085       for (i=0; ((i < nr) && !bFound); i++) {
1086         pi=&(bt[ftype].param[i]);
1087         bFound=
1088           (((pi->AI==-1) ||
1089             (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1090            ((pi->AJ==-1) ||
1091             (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1092            ((pi->AK==-1) ||
1093             (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1094            ((pi->AL==-1) ||
1095             (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1096            );
1097       }
1098     }
1099   } else { /* Not a dihedral */
1100     for (i=0; ((i < nr) && !bFound); i++) {
1101       pi=&(bt[ftype].param[i]);
1102       if (bB)
1103         for (j=0; ((j < nral) &&
1104                    (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
1105           ;
1106       else
1107         for (j=0; ((j < nral) &&
1108                    (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
1109           ;
1110       bFound=(j==nral);
1111     }
1112   }
1113
1114   *param_def = pi;
1115
1116   return bFound;
1117 }
1118
1119
1120
1121 void push_bondnow(t_params *bond, t_param *b)
1122 {
1123   int j;
1124
1125   if (debug) {
1126     fprintf(debug,"push_bondnow: nr = %d\n",bond->nr);
1127     fflush(debug);
1128   }
1129   /* allocate one position extra */
1130   pr_alloc (1,bond);
1131
1132   /* fill the arrays */
1133   for (j=0; (j < MAXFORCEPARAM); j++)
1134     bond->param[bond->nr].c[j]   = b->c[j];
1135   for (j=0; (j < MAXATOMLIST); j++)
1136     bond->param[bond->nr].a[j]   = b->a[j];
1137
1138   bond->nr++;
1139 }
1140
1141 void push_bond(directive d,t_params bondtype[],t_params bond[],
1142                t_atoms *at,t_atomtype atype,char *line,
1143                bool bBonded,bool bGenPairs,real fudgeQQ,
1144                bool bZero,bool *bWarn_copy_A_B)
1145 {
1146   const char *aaformat[MAXATOMLIST]= {
1147     "%d%d",
1148     "%d%d%d",
1149     "%d%d%d%d",
1150     "%d%d%d%d%d",
1151     "%d%d%d%d%d%d"
1152   };
1153   const char *asformat[MAXATOMLIST]= {
1154     "%*s%*s",
1155     "%*s%*s%*s",
1156     "%*s%*s%*s%*s",
1157     "%*s%*s%*s%*s%*s",
1158     "%*s%*s%*s%*s%*s%*s"
1159   };
1160   const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1161   int      nr,i,j,nral,nread,ftype;
1162   char     format[STRLEN];
1163   /* One force parameter more, so we can check if we read too many */
1164   double   cc[MAXFORCEPARAM+1];
1165   int      aa[MAXATOMLIST+1];
1166   t_param  param,paramB,*param_def;
1167   bool     bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
1168   char  errbuf[256];
1169
1170   ftype = ifunc_index(d,1);
1171   nral  = NRAL(ftype);
1172   for(j=0; j<MAXATOMLIST; j++)
1173     aa[j]=NOTSET;
1174   bDef = (NRFP(ftype) > 0);
1175
1176   nread = sscanf(line,aaformat[nral-1],
1177                  &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1178   if (nread < nral) {
1179     too_few();
1180     return;
1181   } else if (nread == nral)
1182     ftype = ifunc_index(d,1);
1183   else {
1184     /* this is a hack to allow for virtual sites with swapped parity */
1185     bSwapParity = (aa[nral]<0);
1186     if (bSwapParity)
1187       aa[nral] = -aa[nral];
1188     ftype = ifunc_index(d,aa[nral]);
1189     if (bSwapParity)
1190       switch(ftype) {
1191       case F_VSITE3FAD:
1192       case F_VSITE3OUT:
1193         break;
1194       default:
1195         gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1196                     interaction_function[F_VSITE3FAD].longname,
1197                     interaction_function[F_VSITE3OUT].longname);
1198       }
1199   }
1200
1201   /* Check for double atoms and atoms out of bounds */
1202   for(i=0; (i<nral); i++) {
1203     if ( aa[i] < 1 || aa[i] > at->nr )
1204       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1205                 "Atom index (%d) in %s out of bounds (1-%d).\n"
1206                 "This probably means that you have inserted topology section \"%s\"\n"
1207                 "in a part belonging to a different molecule than you intended to.\n"
1208                 "In that case move the \"%s\" section to the right molecule.",
1209                 get_warning_file(),get_warning_line(),
1210                 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1211     for(j=i+1; (j<nral); j++)
1212       if (aa[i] == aa[j]) {
1213         sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1214         warning(errbuf);
1215       }
1216   }
1217   if (ftype == F_SETTLE)
1218     if (aa[0]+2 > at->nr)
1219       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1220                   "             Atom index (%d) in %s out of bounds (1-%d)\n"
1221                   "             Settle works on atoms %d, %d and %d",
1222                   get_warning_file(),get_warning_line(),
1223                   aa[0],dir2str(d),at->nr,
1224                   aa[0],aa[0]+1,aa[0]+2);
1225
1226   /* default force parameters  */
1227   for(j=0; (j<MAXATOMLIST); j++)
1228     param.a[j] = aa[j]-1;
1229   for(j=0; (j<MAXFORCEPARAM); j++)
1230     param.c[j] = 0.0;
1231
1232   /* Get force params for normal and free energy perturbation
1233    * studies, as determined by types!
1234    */
1235   if (bBonded) {
1236     bFoundA = default_params(ftype,bondtype,at,atype,&param,FALSE,&param_def);
1237     if (bFoundA) {
1238       /* Copy the A-state and B-state default parameters */
1239       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1240         param.c[j] = param_def->c[j];
1241     }
1242     bFoundB = default_params(ftype,bondtype,at,atype,&param,TRUE,&param_def);
1243     if (bFoundB) {
1244       /* Copy only the B-state default parameters */
1245       for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1246         param.c[j] = param_def->c[j];
1247     }
1248   } else if (ftype == F_LJ14) {
1249     bFoundA = default_nb_params(ftype, bondtype,at,&param,0,FALSE,bGenPairs);
1250     bFoundB = default_nb_params(ftype, bondtype,at,&param,0,TRUE, bGenPairs);
1251   } else if (ftype == F_LJC14_Q) {
1252     param.c[0] = fudgeQQ;
1253     /* Fill in the A-state charges as default parameters */
1254     param.c[1] = at->atom[param.a[0]].q;
1255     param.c[2] = at->atom[param.a[1]].q;
1256     /* The default LJ parameters are the standard 1-4 parameters */
1257     bFoundA = default_nb_params(F_LJ14,bondtype,at,&param,3,FALSE,bGenPairs);
1258     bFoundB = TRUE;
1259   } else if (ftype == F_LJC_PAIRS_NB) {
1260     /* Defaults are not supported here */
1261     bFoundA = FALSE;
1262     bFoundB = TRUE;
1263   } else {
1264     gmx_incons("Unknown function type in push_bond");
1265   }
1266
1267   if (nread > nral) {
1268     strcpy(format,asformat[nral-1]);
1269     strcat(format,ccformat);
1270
1271     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1272                    &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1273
1274     if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1275       {
1276         /* We only have to issue a warning if these atoms are perturbed! */
1277         bPert = FALSE;
1278         for(j=0; (j<nral); j++)
1279           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1280
1281         if (bPert && *bWarn_copy_A_B)
1282           {
1283             sprintf(errbuf,
1284                     "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1285             warning(errbuf);
1286             *bWarn_copy_A_B = FALSE;
1287           }
1288
1289         /* If only the A parameters were specified, copy them to the B state */
1290         /* The B-state parameters correspond to the first nrfpB
1291          * A-state parameters.
1292          */
1293         for(j=0; (j<NRFPB(ftype)); j++)
1294           cc[nread++] = cc[j];
1295       }
1296
1297     /* If nread was 0 or EOF, no parameters were read => use defaults.
1298      * If nread was nrfpA we copied above so nread=nrfp.
1299      * If nread was nrfp we are cool.
1300      * For F_LJC14_Q we allow supplying fudgeQQ only.
1301      * Anything else is an error!
1302      */
1303     if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1304         !(ftype == F_LJC14_Q && nread == 1))
1305       {
1306         gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1307                   nread,NRFPA(ftype),NRFP(ftype),
1308                   interaction_function[ftype].longname);
1309       }
1310
1311     for(j=0; (j<nread); j++)
1312       param.c[j]=cc[j];
1313
1314     /* Check whether we have to use the defaults */
1315     if (nread == NRFP(ftype))
1316       bDef = FALSE;
1317   }
1318   else
1319   {
1320     nread = 0;
1321   }
1322
1323   /* nread now holds the number of force parameters read! */
1324
1325   if (bDef) {
1326     /* Use defaults */
1327
1328     if (nread > 0 && nread < NRFPA(ftype)) {
1329       /* Issue an error, do not use defaults */
1330       sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1331       warning_error(errbuf);
1332     } else {
1333       if (nread == 0 || nread == EOF) {
1334         if (!bFoundA) {
1335           if (interaction_function[ftype].flags & IF_VSITE) {
1336             /* set them to NOTSET, will be calculated later */
1337             for(j=0; (j<MAXFORCEPARAM); j++)
1338               param.c[j] = NOTSET;
1339             if (bSwapParity)
1340               param.C1 = -1; /* flag to swap parity of vsite construction */
1341           } else {
1342             if (bZero) {
1343               fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1344                       interaction_function[ftype].longname);
1345             } else {
1346               sprintf(errbuf,"No default %s types",
1347                       interaction_function[ftype].longname);
1348               warning_error(errbuf);
1349             }
1350           }
1351         } else
1352           if (bSwapParity)
1353             switch(ftype) {
1354             case F_VSITE3FAD:
1355               param.C0 = 360 - param.C0;
1356               break;
1357             case F_VSITE3OUT:
1358               param.C2 = -param.C2;
1359               break;
1360             }
1361       }
1362       if (!bFoundB) {
1363         /* We only have to issue a warning if these atoms are perturbed! */
1364         bPert = FALSE;
1365         for(j=0; (j<nral); j++)
1366           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1367
1368         if (bPert) {
1369           sprintf(errbuf,"No default %s types for perturbed atoms, "
1370                   "using normal values",interaction_function[ftype].longname);
1371           warning(errbuf);
1372         }
1373       }
1374     }
1375   }
1376
1377   if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1378       && param.c[5]!=param.c[2])
1379     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1380                 "             %s multiplicity can not be perturbed %f!=%f",
1381                 get_warning_file(),get_warning_line(),
1382                 interaction_function[ftype].longname,
1383                 param.c[2],param.c[5]);
1384
1385   if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1386     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1387                 "             %s table number can not be perturbed %d!=%d",
1388                 get_warning_file(),get_warning_line(),
1389                 interaction_function[ftype].longname,
1390                 (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1391
1392   /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1393   if (ftype==F_RBDIHS) {
1394     nr=0;
1395     for(i=0;i<NRFP(ftype);i++) {
1396       if(param.c[i]!=0)
1397         nr++;
1398     }
1399     if(nr==0)
1400       return;
1401   }
1402
1403   /* Put the values in the appropriate arrays */
1404   push_bondnow (&bond[ftype],&param);
1405 }
1406
1407 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1408                   t_atoms *at,t_atomtype atype,char *line)
1409 {
1410   char *ptr;
1411   int  type,ftype,j,n,ret,nj,a;
1412   int  *atc=NULL;
1413   double *weight=NULL,weight_tot;
1414   t_param param;
1415
1416   /* default force parameters  */
1417   for(j=0; (j<MAXATOMLIST); j++)
1418     param.a[j] = NOTSET;
1419   for(j=0; (j<MAXFORCEPARAM); j++)
1420     param.c[j] = 0.0;
1421
1422   ptr = line;
1423   ret = sscanf(ptr,"%d%n",&a,&n);
1424   ptr += n;
1425   if (ret == 0)
1426     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1427               "             Expected an atom index in section \"%s\"",
1428                 get_warning_file(),get_warning_line(),
1429               dir2str(d));
1430
1431   param.a[0] = a - 1;
1432
1433   ret = sscanf(ptr,"%d%n",&type,&n);
1434   ptr += n;
1435   ftype = ifunc_index(d,type);
1436
1437   weight_tot = 0;
1438   nj = 0;
1439   do {
1440     ret = sscanf(ptr,"%d%n",&a,&n);
1441     ptr += n;
1442     if (ret > 0) {
1443       if (nj % 20 == 0) {
1444         srenew(atc,nj+20);
1445         srenew(weight,nj+20);
1446       }
1447       atc[nj] = a - 1;
1448       switch (type) {
1449       case 1:
1450         weight[nj] = 1;
1451         break;
1452       case 2:
1453         /* Here we use the A-state mass as a parameter.
1454          * Note that the B-state mass has no influence.
1455          */
1456         weight[nj] = at->atom[atc[nj]].m;
1457         break;
1458       case 3:
1459         weight[nj] = -1;
1460         ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1461         ptr += n;
1462         if (weight[nj] < 0)
1463           gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1464                     "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1465                     get_warning_file(),get_warning_line(),
1466                     nj+1,atc[nj]+1);
1467         break;
1468       default:
1469         gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1470       }
1471       weight_tot += weight[nj];
1472       nj++;
1473     }
1474   } while (ret > 0);
1475
1476   if (nj == 0)
1477     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1478               "             Expected more than one atom index in section \"%s\"",
1479               get_warning_file(),get_warning_line(),
1480               dir2str(d));
1481
1482   if (weight_tot == 0)
1483      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1484                "             The total mass of the construting atoms is zero",
1485                get_warning_file(),get_warning_line());
1486
1487   for(j=0; j<nj; j++) {
1488     param.a[1] = atc[j];
1489     param.c[0] = nj;
1490     param.c[1] = weight[j]/weight_tot;
1491     /* Put the values in the appropriate arrays */
1492     push_bondnow (&bond[ftype],&param);
1493   }
1494
1495   sfree(atc);
1496   sfree(weight);
1497 }
1498
1499 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1500               int *nrcopies)
1501 {
1502   int  i,copies;
1503   char type[STRLEN];
1504
1505   *nrcopies=0;
1506   if (sscanf(pline,"%s%d",type,&copies) != 2) {
1507     too_few();
1508     return;
1509   }
1510
1511   /* search moleculename */
1512   for (i=0; ((i<nrmols) && strcasecmp(type,*(mols[i].name))); i++)
1513     ;
1514
1515   if (i<nrmols) {
1516     *nrcopies        = copies;
1517     *whichmol        = i;
1518   } else
1519     gmx_fatal(FARGS,"No such moleculetype %s",type);
1520 }
1521
1522 void init_block2(t_block2 *b2, int natom)
1523 {
1524   int i;
1525
1526   b2->nr=natom;
1527   snew(b2->nra,b2->nr);
1528   snew(b2->a,b2->nr);
1529   for(i=0; (i<b2->nr); i++)
1530     b2->a[i]=NULL;
1531 }
1532
1533 void done_block2(t_block2 *b2)
1534 {
1535   int i;
1536
1537   if (b2->nr) {
1538     for(i=0; (i<b2->nr); i++)
1539       sfree(b2->a[i]);
1540     sfree(b2->a);
1541     sfree(b2->nra);
1542     b2->nr=0;
1543   }
1544 }
1545
1546 void push_excl(char *line, t_block2 *b2)
1547 {
1548   int  i,j;
1549   int  n;
1550   char base[STRLEN],format[STRLEN];
1551
1552   if (sscanf(line,"%d",&i)==0)
1553     return;
1554
1555   if ((1 <= i) && (i <= b2->nr))
1556     i--;
1557   else {
1558     if (debug)
1559       fprintf(debug,"Unbound atom %d\n",i-1);
1560     return;
1561   }
1562   strcpy(base,"%*d");
1563   do {
1564     strcpy(format,base);
1565     strcat(format,"%d");
1566     n=sscanf(line,format,&j);
1567     if (n == 1) {
1568       if ((1 <= j) && (j <= b2->nr)) {
1569         j--;
1570         srenew(b2->a[i],++(b2->nra[i]));
1571         b2->a[i][b2->nra[i]-1]=j;
1572         /* also add the reverse exclusion! */
1573         srenew(b2->a[j],++(b2->nra[j]));
1574         b2->a[j][b2->nra[j]-1]=i;
1575         strcat(base,"%*d");
1576       }
1577       else
1578         gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
1579     }
1580   } while (n == 1);
1581 }
1582
1583 void b_to_b2(t_blocka *b, t_block2 *b2)
1584 {
1585   int     i;
1586   atom_id j,a;
1587
1588   for(i=0; (i<b->nr); i++)
1589     for(j=b->index[i]; (j<b->index[i+1]); j++) {
1590       a=b->a[j];
1591       srenew(b2->a[i],++b2->nra[i]);
1592       b2->a[i][b2->nra[i]-1]=a;
1593     }
1594 }
1595
1596 void b2_to_b(t_block2 *b2, t_blocka *b)
1597 {
1598   int     i,nra;
1599   atom_id j;
1600
1601   nra=0;
1602   for(i=0; (i<b2->nr); i++) {
1603     b->index[i]=nra;
1604     for(j=0; (j<b2->nra[i]); j++)
1605       b->a[nra+j]=b2->a[i][j];
1606     nra+=b2->nra[i];
1607   }
1608   /* terminate list */
1609   b->index[i]=nra;
1610 }
1611
1612 static int icomp(const void *v1, const void *v2)
1613 {
1614   return (*((atom_id *) v1))-(*((atom_id *) v2));
1615 }
1616
1617 void merge_excl(t_blocka *excl, t_block2 *b2)
1618 {
1619   int     i,k;
1620   atom_id j;
1621   int     nra;
1622
1623   if (!b2->nr)
1624     return;
1625   else if (b2->nr != excl->nr) {
1626     gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
1627                 b2->nr,excl->nr);
1628   }
1629   else if (debug)
1630     fprintf(debug,"Entering merge_excl\n");
1631
1632   /* First copy all entries from excl to b2 */
1633   b_to_b2(excl,b2);
1634
1635   /* Count and sort the exclusions */
1636   nra=0;
1637   for(i=0; (i<b2->nr); i++) {
1638     if (b2->nra[i] > 0) {
1639       /* remove double entries */
1640       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
1641       k=1;
1642       for(j=1; (j<b2->nra[i]); j++)
1643         if (b2->a[i][j]!=b2->a[i][k-1]) {
1644           b2->a[i][k]=b2->a[i][j];
1645           k++;
1646         }
1647       b2->nra[i]=k;
1648       nra+=b2->nra[i];
1649     }
1650   }
1651   excl->nra=nra;
1652   srenew(excl->a,excl->nra);
1653
1654   b2_to_b(b2,excl);
1655 }
1656
1657 int add_atomtype_decoupled(t_symtab *symtab,t_atomtype at,
1658                            t_nbparam ***nbparam,t_nbparam ***pair)
1659 {
1660   t_atom  atom;
1661   t_param param;
1662   int     i,nr;
1663
1664   /* Define an atom type with all parameters set to zero (no interactions) */
1665   atom.q = 0.0;
1666   atom.m = 0.0;
1667   /* Type for decoupled atoms could be anything,
1668    * this should be changed automatically later when required.
1669    */
1670   atom.ptype = eptAtom;
1671   for (i=0; (i<MAXFORCEPARAM); i++)
1672     param.c[i] = 0.0;
1673
1674   nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0);
1675
1676   /* Add space in the non-bonded parameters matrix */
1677   realloc_nb_params(at,nbparam,pair);
1678
1679   return nr;
1680 }
1681
1682 static void convert_pairs_to_pairsQ(t_params *plist,
1683                                     real fudgeQQ,t_atoms *atoms)
1684 {
1685   t_param *param;
1686   int  i;
1687   real v,w;
1688
1689   /* Copy the pair list to the pairQ list */
1690   plist[F_LJC14_Q] = plist[F_LJ14];
1691   /* Empty the pair list */
1692   plist[F_LJ14].nr    = 0;
1693   plist[F_LJ14].param = NULL;
1694   param = plist[F_LJC14_Q].param;
1695   for(i=0; i<plist[F_LJC14_Q].nr; i++) {
1696     v = param[i].c[0];
1697     w = param[i].c[1];
1698     param[i].c[0] = fudgeQQ;
1699     param[i].c[1] = atoms->atom[param[i].a[0]].q;
1700     param[i].c[2] = atoms->atom[param[i].a[1]].q;
1701     param[i].c[3] = v;
1702     param[i].c[4] = w;
1703   }
1704 }
1705
1706 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
1707 {
1708   int n,ntype,i,j,k;
1709   t_atom *atom;
1710   t_blocka *excl;
1711   bool bExcl;
1712   t_param param;
1713
1714   n = mol->atoms.nr;
1715   atom = mol->atoms.atom;
1716
1717   ntype = sqrt(nbp->nr);
1718
1719   for (i=0; i<MAXATOMLIST; i++)
1720     param.a[i] = NOTSET;
1721   for (i=0; i<MAXFORCEPARAM; i++)
1722     param.c[i] = NOTSET;
1723
1724   /* Add a pair interaction for all non-excluded atom pairs */
1725   excl = &mol->excls;
1726   for(i=0; i<n; i++) {
1727     for(j=i+1; j<n; j++) {
1728       bExcl = FALSE;
1729       for(k=excl->index[i]; k<excl->index[i+1]; k++) {
1730         if (excl->a[k] == j)
1731           bExcl = TRUE;
1732       }
1733       if (!bExcl) {
1734         if (nb_funct != F_LJ)
1735           gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
1736         param.a[0] = i;
1737         param.a[1] = j;
1738         param.c[0] = atom[i].q;
1739         param.c[1] = atom[j].q;
1740         param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
1741         param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
1742         push_bondnow(&mol->plist[F_LJC_PAIRS_NB],&param);
1743       }
1744     }
1745   }
1746 }
1747
1748 static void set_excl_all(t_blocka *excl)
1749 {
1750   int nat,i,j,k;
1751
1752   /* Get rid of the current exclusions and exclude all atom pairs */
1753   nat = excl->nr;
1754   excl->nra = nat*nat;
1755   srenew(excl->a,excl->nra);
1756   k = 0;
1757   for(i=0; i<nat; i++) {
1758     excl->index[i] = k;
1759     for(j=0; j<nat; j++)
1760       excl->a[k++] = j;
1761   }
1762   excl->index[nat] = k;
1763 }
1764
1765 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
1766                            int couple_lam0,int couple_lam1)
1767 {
1768   int i;
1769
1770   for(i=0; i<atoms->nr; i++) {
1771     if (couple_lam0 != ecouplamVDWQ)
1772       atoms->atom[i].q     = 0.0;
1773     if (couple_lam0 == ecouplamNONE)
1774       atoms->atom[i].type  = atomtype_decouple;
1775     if (couple_lam1 != ecouplamVDWQ)
1776       atoms->atom[i].qB    = 0.0;
1777     if (couple_lam1 == ecouplamNONE)
1778       atoms->atom[i].typeB = atomtype_decouple;
1779   }
1780 }
1781
1782 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
1783                             int couple_lam0,int couple_lam1,
1784                             bool bCoupleIntra,int nb_funct,t_params *nbp)
1785 {
1786   convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
1787   if (!bCoupleIntra) {
1788     generate_LJCpairsNB(mol,nb_funct,nbp);
1789     set_excl_all(&mol->excls);
1790   }
1791   decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
1792 }