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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_toppush_c = "$Id$";
  30
  31 #ifdef HAVE_IDENT
  32 #ident "@(#) toppush.c 1.72 9/30/97"
  33 #endif
  34
  35 #include <math.h>
  36
  37 #include "sysstuff.h"
  38 #include "smalloc.h"
  39 #include "macros.h"
  40 #include "assert.h"
  41 #include "string2.h"
  42 #include "names.h"
  43 #include "toputil.h"
  44 #include "toppush.h"
  45 #include "topdirs.h"
  46 #include "symtab.h"
  47 #include "fatal.h"
  48
  49 static char  errbuf[256];
  50
  51 void generate_nbparams(int comb,int ftype,t_params *plist,t_atomtype *atype)
  52 {
  53   int   i,j,k=-1,nf;
  54   int   nr,nrfp;
  55   real  c,sig6,sigma_ij,eps_ij,bi,bj;
  56
  57   /* Lean mean shortcuts */
  58   nr   = atype->nr;
  59   nrfp = NRFP(ftype);
  60   snew(plist->param,nr*nr);
  61   plist->nr=nr*nr;
  62
  63   /* Fill the matrix with force parameters */
  64   switch(ftype) {
  65   case F_LJ:
  66     switch (comb) {
  67     case eCOMB_ARITHMETIC:
  68       /* Gromos rules */
  69       for(i=k=0; (i<nr); i++)
  70         for(j=0; (j<nr); j++,k++)
  71           for(nf=0; (nf<nrfp); nf++) {
  72             c = sqrt(atype->nb[i].c[nf] * atype->nb[j].c[nf]);
  73             plist->param[k].c[nf]      = c;
  74           }
  75       break;
  76
  77     case eCOMB_GEOMETRIC:
  78       /* c0 and c1 are epsilon and sigma */
  79       for (i=k=0; (i < nr); i++)
  80         for (j=0; (j < nr); j++,k++) {
  81           sigma_ij = (atype->nb[i].C0+atype->nb[j].C0)*0.5;
  82           eps_ij   = sqrt(atype->nb[i].C1*atype->nb[j].C1);
  83           sig6     = pow(sigma_ij,6.0);
  84
  85           plist->param[k].c[0] = 4*eps_ij*sig6;
  86           plist->param[k].c[1] = 4*eps_ij*sig6*sig6;
  87         }
  88
  89       break;
  90     case eCOMB_ARITH_SIG_EPS:
  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           sigma_ij = sqrt(atype->nb[i].C0*atype->nb[j].C0);
  95           eps_ij   = sqrt(atype->nb[i].C1*atype->nb[j].C1);
  96           sig6     = pow(sigma_ij,6.0);
  97
  98           plist->param[k].c[0] = 4*eps_ij*sig6;
  99           plist->param[k].c[1] = 4*eps_ij*sig6*sig6;
 100         }
 101
 102       break;
 103     default:
 104       fatal_error(0,"No such combination rule %d",comb);
 105     }
 106     assert(plist->nr == k);
 107     break;
 108
 109   case F_BHAM:
 110     /* Buckingham rules */
 111     for(i=k=0; (i<nr); i++)
 112       for(j=0; (j<nr); j++,k++) {
 113         plist->param[k].c[0] = sqrt(atype->nb[i].c[0] * atype->nb[j].c[0]);
 114         bi = atype->nb[i].c[1];
 115         bj = atype->nb[j].c[1];
 116         if ((bi == 0) || (bj == 0))
 117           plist->param[k].c[1] = 0;
 118         else
 119           plist->param[k].c[1] = 2.0/(1/bi+1/bj);
 120         plist->param[k].c[2] = sqrt(atype->nb[i].c[2] * atype->nb[j].c[2]);
 121       }
 122
 123     break;
 124   default:
 125     sprintf(errbuf,"Invalid nonbonded type %s",
 126             interaction_function[ftype].longname);
 127     warning(errbuf);
 128   }
 129 }
 130
 131 void push_at (t_symtab *symtab, t_atomtype *at, char *line,int nb_funct)
 132 {
 133   typedef struct {
 134     char *entry;
 135     int  ptype;
 136   } t_xlate;
 137   t_xlate xl[eptNR] = {
 138     { "A",   eptAtom },
 139     { "N",   eptNucleus },
 140     { "S",   eptShell },
 141     { "B",   eptBond },
 142     { "D",   eptDummy },
 143   };
 144
 145   int    nr,i,j,pt,nfp0=-1;
 146   char   type[STRLEN],ptype[STRLEN];
 147   double m,q;
 148   double c[MAXFORCEPARAM];
 149
 150   nr = at->nr;
 151
 152   switch (nb_funct) {
 153   case F_LJ:
 154     nfp0 = 2;
 155     if (sscanf (line,"%s%lf%lf%s%lf%lf",type,&m,&q,ptype,&c[0],&c[1]) != 6) {
 156       too_few();
 157       return;
 158     }
 159     break;
 160   case F_BHAM:
 161     nfp0 = 3;
 162     if (sscanf (line,"%s%lf%lf%s%lf%lf%lf",
 163                 type,&m,&q,ptype,&c[0],&c[1],&c[2]) != 7) {
 164       too_few();
 165       return;
 166     }
 167     break;
 168   default:
 169     fatal_error(0,"Invalid function type %d in push_at %s %d",nb_funct,
 170                 __FILE__,__LINE__);
 171   }
 172   for(j=nfp0; (j<MAXFORCEPARAM); j++)
 173     c[j]=0.0;
 174
 175   for(j=0; (j<eptNR); j++)
 176     if (strcasecmp(ptype,xl[j].entry) == 0)
 177       break;
 178   if (j == eptNR)
 179     fatal_error(0,"Invalid particle type %s on line %s",
 180                 ptype,line);
 181   pt=xl[j].ptype;
 182   if (debug)
 183     fprintf(debug,"ptype: %s\n",ptype_str[pt]);
 184
 185   /* Test if this type overwrites another */
 186   for (i=0; ((i<nr) && (strcasecmp(*(at->atomname[i]),type) != 0)); i++)
 187     ;
 188
 189   if ((i==nr) || (nr==0)) {
 190     /* New atomtype, get new space for arrays */
 191     srenew(at->atom,nr+1);
 192     srenew(at->atomname,nr+1);
 193     srenew(at->nb,nr+1);
 194     at->nr++;
 195   }
 196   else {
 197     sprintf(errbuf,"Overriding atomtype %s",type);
 198     warning(errbuf);
 199     nr = i;
 200   }
 201   /* fill the arrays */
 202   at->atomname[nr] = put_symtab(symtab,type);
 203   at->atom[nr].ptype = pt;
 204   at->atom[nr].m = m;
 205   at->atom[nr].q = q;
 206
 207   for (i=0; (i<MAXFORCEPARAM); i++)
 208     at->nb[nr].c[i] = c[i];
 209 }
 210
 211 static void push_bondtype(t_params *bt,t_param *b,int nral,int ftype)
 212 {
 213   int  i,j;
 214   bool bTest,bFound;
 215   int  nr   = bt->nr;
 216   int  nrfp = NRFP(ftype);
 217
 218   /* Check if this entry overwrites another */
 219   bFound=FALSE;
 220   for (i=0; (i < nr); i++) {
 221     bTest = TRUE;
 222     for (j=0; (j < nral); j++)
 223       bTest=(bTest && (b->a[j] == bt->param[i].a[j]));
 224     if (!bTest) {
 225       bTest=TRUE;
 226       for(j=0; (j<nral); j++)
 227         bTest=(bTest && (b->a[nral-1-j] == bt->param[i].a[j]));
 228     }
 229     if (bTest) {
 230       /* Overwrite it! */
 231       for (j=0; (j < nrfp); j++)
 232         bt->param[i].c[j] = b->c[j];
 233       bFound=TRUE;
 234     }
 235   }
 236   if (!bFound) {
 237     /* alloc */
 238     pr_alloc (2,bt);
 239
 240     /* fill the arrays up and down */
 241     memcpy((char *) bt->param[bt->nr].c,  (char *) b->c,(size_t)sizeof(b->c));
 242     memcpy((char *) bt->param[bt->nr].a,  (char *) b->a,(size_t)sizeof(b->a));
 243     memcpy((char *) bt->param[bt->nr+1].c,(char *) b->c,(size_t)sizeof(b->c));
 244     for (j=0; (j < nral); j++)
 245       bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
 246
 247     bt->nr += 2;
 248   }
 249 }
 250
 251 void push_bt(directive d,t_params bt[],int nral,t_atomtype *at,char *line)
 252 {
 253   static   char *formal[MAXATOMLIST+1] = {
 254     "%s",
 255     "%s%s",
 256     "%s%s%s",
 257     "%s%s%s%s",
 258     "%s%s%s%s%s",
 259     "%s%s%s%s%s%s"
 260   };
 261   static   char *formnl[MAXATOMLIST+1] = {
 262     "%*s",
 263     "%*s%*s",
 264     "%*s%*s%*s",
 265     "%*s%*s%*s%*s",
 266     "%*s%*s%*s%*s%*s",
 267     "%*s%*s%*s%*s%*s%*s"
 268   };
 269   static   char *formlf[MAXFORCEPARAM] = {
 270     "%lf",
 271     "%lf%lf",
 272     "%lf%lf%lf",
 273     "%lf%lf%lf%lf",
 274     "%lf%lf%lf%lf%lf",
 275     "%lf%lf%lf%lf%lf%lf",
 276   };
 277   int      i,ft,ftype,nn,nrfp;
 278   char     f1[STRLEN];
 279   char     alc[MAXATOMLIST+1][20];
 280   double   c[MAXFORCEPARAM];
 281   t_param  p;
 282
 283   /* Make format string (nral ints+functype) */
 284   if ((nn=sscanf(line,formal[nral],
 285                  alc[0],alc[1],alc[2],alc[3],alc[4],alc[5])) != nral+1) {
 286     sprintf(errbuf,"Not enough atomtypes (%d instead of %d)",nn-1,nral);
 287     warning(errbuf);
 288     return;
 289   }
 290
 291   ft    = atoi(alc[nral]);
 292   ftype = ifunc_index(d,ft);
 293   nrfp  = NRFP(ftype);
 294   strcpy(f1,formnl[nral]);
 295   strcat(f1,formlf[nrfp-1]);
 296   if ((nn=sscanf(line,f1,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5]))
 297       != nrfp) {
 298     for( ; (nn<nrfp); nn++)
 299       c[nn] = 0.0;
 300   }
 301   for(i=0; (i<nral); i++)
 302     p.a[i]=at2type(alc[i],at);
 303   for(i=0; (i<nrfp); i++)
 304     p.c[i]=c[i];
 305   push_bondtype (&(bt[ftype]),&p,nral,ftype);
 306 }
 307
 308 void push_nbt(directive d,t_params nbt[],t_atomtype *atype,
 309               char *pline,int nb_funct)
 310 {
 311   /* swap the atoms */
 312   static char *form2="%*s%*s%*s%lf%lf";
 313   static char *form3="%*s%*s%*s%lf%lf%lf";
 314   char    a0[80],a1[80];
 315   int     i,f,k,ftype,atnr,nrfp;
 316   double  c[3];
 317   atom_id ai,aj;
 318
 319   if (sscanf (pline,"%s%s%d",a0,a1,&f) != 3) {
 320     too_few();
 321     return;
 322   }
 323
 324   ftype=ifunc_index(d,f);
 325
 326   if (ftype != nb_funct) {
 327     sprintf(errbuf,"Trying to add %s while the default nonbond type is %s",
 328             interaction_function[ftype].longname,
 329             interaction_function[nb_funct].longname);
 330     warning(errbuf);
 331     return;
 332   }
 333
 334   /* Get the force parameters */
 335   if (NRFP(ftype) == 2) {
 336     if (sscanf(pline,form2,&c[0],&c[1]) != 2) {
 337       too_few();
 338       return;
 339     }
 340   }
 341   else {
 342     if (sscanf(pline,form3,&c[0],&c[1],&c[2]) != 3) {
 343       too_few();
 344       return;
 345     }
 346   }
 347   /* Put the parameters in the matrix */
 348   ai = at2type (a0,atype);
 349   aj = at2type (a1,atype);
 350
 351   atnr = atype->nr;
 352   nrfp = NRFP(ftype);
 353   for (i=0; (i < nrfp); i++) {
 354     k = atnr*ai+aj;
 355     nbt[ftype].param[k].c[i] = c[i];
 356     k = atnr*aj+ai;
 357     nbt[ftype].param[k].c[i] = c[i];
 358   }
 359 }
 360
 361 static void push_atom_now(t_symtab *symtab,t_atoms *at,int atomnr,
 362                           int type,int ptype,int resnumber,int cgnumber,
 363                           char *resname,char *name,real m0,real q0,
 364                           int typeB,real mB,real qB)
 365 {
 366   int j,resnr_diff=-1;
 367   int nr = at->nr;
 368
 369   if (((nr==0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
 370     fatal_error(0,"Atoms in the .top are not numbered consecutively from 1\n");
 371   if (nr)
 372     resnr_diff = resnumber - (at->atom[at->nr-1].resnr+1);
 373   if (((nr==0) && (resnumber != 1)) ||
 374       (nr && (resnr_diff != 0) && (resnr_diff != 1)))
 375     fatal_error(0,"Residue numbers in the .top are not numbered consecutively from 1\n");
 376
 377   /* New atom instance
 378    * get new space for arrays
 379    */
 380   srenew(at->atom,nr+1);
 381   srenew(at->atomname,nr+1);
 382
 383   if (resnumber > at->nres) {
 384     at->nres=resnumber;
 385     srenew(at->resname,resnumber);
 386     at->resname[resnumber-1] = put_symtab(symtab,resname);
 387   }
 388   /* fill the list */
 389   at->atom[nr].type  = type;
 390   at->atom[nr].ptype = ptype;
 391   at->atom[nr].q     = q0;
 392   at->atom[nr].m     = m0;
 393   at->atom[nr].typeB = typeB;
 394   at->atom[nr].qB    = qB;
 395   at->atom[nr].mB    = mB;
 396
 397   for(j=0; (j<egcNR); j++)
 398     at->atom[nr].grpnr[j] = -1;
 399   at->atom[nr].resnr = resnumber-1;
 400   at->atomname[nr] = put_symtab(symtab,name);
 401   at->nr++;
 402 }
 403
 404 void push_cg(t_block *block, int *lastindex, int index, int a)
 405 {
 406   if (debug)
 407     fprintf (debug,"Index %d, Atom %d\n",index,a);
 408
 409   if (((block->nr) && (*lastindex != index)) || (!block->nr)) {
 410     /* add a new block */
 411     block->nr++;
 412     srenew(block->index,block->nr+1);
 413   }
 414   srenew(block->a,block->nra+1);
 415   block->a[block->nra++]=a;
 416   block->index[block->nr]=block->nra;
 417   *lastindex = index;
 418 }
 419
 420 void push_atom(t_symtab *symtab,t_block *cgs,
 421                t_atoms *at,t_atomtype *atype,char *line)
 422 {
 423   int           nr,ptype;
 424   static int    lastcg;
 425   int           resnumber,cgnumber,atomnr,type,typeB,nscan;
 426   char          id[STRLEN],ctype[STRLEN],
 427                 resname[STRLEN],name[STRLEN];
 428   double        m,q,mb,qb;
 429   real          m0,q0,mB,qB;
 430
 431   /* Make a shortcut for writing in this molecule  */
 432   nr = at->nr;
 433
 434   /* Fixed parameters */
 435   if (sscanf(line,"%s%s%d%s%s%d",
 436              id,ctype,&resnumber,resname,name,&cgnumber) != 6) {
 437     too_few();
 438     return;
 439   }
 440   sscanf(id,"%d",&atomnr);
 441   type  = at2type(ctype,atype);
 442   ptype = atype->atom[type].ptype;
 443
 444   /* Set default from type */
 445   q0    = atype->atom[type].q;
 446   m0    = atype->atom[type].m;
 447   typeB = type;
 448   qB    = q0;
 449   mB    = m0;
 450
 451   /* Optional parameters */
 452   nscan=sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf",
 453                &q,&m,ctype,&qb,&mb);
 454
 455   /* Nasty switch that falls thru all the way down! */
 456   if (nscan > 0) {
 457     q0 = qB = q;
 458     if (nscan > 1) {
 459       m0 = mB = m;
 460       if (nscan > 2) {
 461         typeB=at2type(ctype,atype);
 462         qB = atype->atom[typeB].q;
 463         mB = atype->atom[typeB].m;
 464         if (nscan > 3) {
 465           qB = qb;
 466           if (nscan > 4)
 467             mB = mb;
 468         }
 469       }
 470     }
 471   }
 472   if (debug)
 473     fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
 474
 475   push_cg(cgs,&lastcg,cgnumber,nr);
 476
 477   push_atom_now(symtab,at,atomnr,type,ptype,resnumber,cgnumber,
 478                 resname,name,m0,q0,typeB,mB,qB);
 479 }
 480
 481 void push_molt(t_symtab *symtab,t_molinfo *newmol,char *line)
 482 {
 483   char type[STRLEN];
 484   int nrexcl;
 485
 486   if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
 487     too_few();
 488     return;
 489   }
 490
 491   /* Fill in the values */
 492   newmol->name     = put_symtab(symtab,type);
 493   newmol->nrexcl   = nrexcl;
 494   newmol->excl_set = FALSE;
 495 }
 496
 497 static bool default_params(int ftype,t_params bt[],t_atoms *at,t_param *p,
 498                            bool bB)
 499 {
 500   int      i,j;
 501   bool     bFound;
 502   t_param  *pi=NULL;
 503   int      nr   = bt[ftype].nr;
 504   int      nral = NRAL(ftype);
 505   int      nrfp = interaction_function[ftype].nrfpA;
 506   int      nrfpB= interaction_function[ftype].nrfpB;
 507
 508   bFound=FALSE;
 509   for (i=0; ((i < nr) && !bFound); i++) {
 510     pi=&(bt[ftype].param[i]);
 511     if ((ftype == F_PDIHS)  || (ftype == F_RBDIHS)) {
 512       /* The j and k atoms are decisive about which dihedral type
 513        * we should take.
 514        */
 515       if (bB)
 516         bFound=((at->atom[p->AJ].typeB==pi->AI) &&
 517                 (at->atom[p->AK].typeB==pi->AJ));
 518       else
 519         bFound=((at->atom[p->AJ].type==pi->AI) &&
 520                 (at->atom[p->AK].type==pi->AJ));
 521     }
 522     else if (ftype == F_IDIHS) {
 523       /* The i and l atoms are decisive about which dihedral type
 524        * we should take.
 525        */
 526       if (bB)
 527         bFound=((at->atom[p->AI].typeB==pi->AI) &&
 528                 (at->atom[p->AL].typeB==pi->AJ));
 529       else
 530         bFound=((at->atom[p->AI].type==pi->AI) &&
 531                 (at->atom[p->AL].type==pi->AJ));
 532     }
 533     else {
 534       if (bB)
 535         for (j=0; ((j < nral) &&
 536                    (at->atom[p->a[j]].typeB == pi->a[j])); j++);
 537       else
 538         for (j=0; ((j < nral) &&
 539                    (at->atom[p->a[j]].type == pi->a[j])); j++);
 540       bFound=(j==nral);
 541     }
 542   }
 543   if (bFound) {
 544     if (bB) {
 545       assert(nrfp+nrfpB <= MAXFORCEPARAM);
 546       for (j=0; (j < nrfpB); j++)
 547         p->c[nrfp+j] = pi->c[j];
 548     }
 549     else
 550       for (j=0; (j < nrfp); j++)
 551         p->c[j] = pi->c[j];
 552   }
 553   else {
 554     for (j=0; (j < nrfp); j++)
 555       p->c[j] = 0.0;
 556   }
 557   return bFound;
 558 }
 559
 560 void push_bondnow(t_params *bond, t_param *b)
 561 {
 562   if (debug) {
 563     fprintf(debug,"push_bondnow: nr = %d\n",bond->nr);
 564     fflush(debug);
 565   }
 566   /* allocate one position extra */
 567   pr_alloc (1,bond);
 568
 569   /* fill the arrays */
 570   memcpy ((char *)&(bond->param[bond->nr]),(char *)b,(size_t)sizeof(t_param));
 571
 572   bond->nr++;
 573 }
 574
 575 void push_bond(directive d,t_params bondtype[],t_params bond[],
 576                t_atoms *at,char *line)
 577 {
 578   static char *aaformat[MAXATOMLIST]= {
 579     "%d%d",
 580     "%d%d%d",
 581     "%d%d%d%d",
 582     "%d%d%d%d%d",
 583     "%d%d%d%d%d%d"
 584   };
 585   static char *asformat[MAXATOMLIST]= {
 586     "%*s%*s",
 587     "%*s%*s%*s",
 588     "%*s%*s%*s%*s",
 589     "%*s%*s%*s%*s%*s",
 590     "%*s%*s%*s%*s%*s%*s"
 591   };
 592   static char *ccformat[MAXFORCEPARAM]= {
 593     "%lf",
 594     "%lf%lf",
 595     "%lf%lf%lf",
 596     "%lf%lf%lf%lf",
 597     "%lf%lf%lf%lf%lf",
 598     "%lf%lf%lf%lf%lf%lf"
 599   };
 600   int      nr,i,j,nrfp,nrfpA,nral,nread,ftype;
 601   char     format[STRLEN];
 602   double   cc[MAXFORCEPARAM];
 603   int      aa[MAXATOMLIST+1];
 604   t_param  param,paramB;
 605   bool     bFoundA,bFoundB,bDef,bPert,bSwapParity=FALSE;
 606
 607   ftype = ifunc_index(d,1);
 608   nral  = NRAL(ftype);
 609   for(j=0; j<MAXATOMLIST; j++)
 610     aa[j]=NOTSET;
 611   nread = sscanf(line,aaformat[nral-1],
 612                  &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
 613   if (nread < nral) {
 614     too_few();
 615     return;
 616   } else if (nread == nral)
 617     ftype = ifunc_index(d,1);
 618   else {
 619     /* this is a hack to allow for dummy atoms with swapped parity */
 620     bSwapParity = (aa[nral]<0);
 621     if (bSwapParity)
 622       aa[nral] = -aa[nral];
 623     ftype = ifunc_index(d,aa[nral]);
 624     if (bSwapParity)
 625       switch(ftype) {
 626       case F_DUMMY3FAD:
 627       case F_DUMMY3OUT:
 628         break;
 629       default:
 630         fatal_error(0,"Negative function types only allowed for %s and %s",
 631                     interaction_function[F_DUMMY3FAD].longname,
 632                     interaction_function[F_DUMMY3OUT].longname);
 633       }
 634   }
 635
 636   /* Check for double atoms */
 637   for(i=0; (i<nral); i++)
 638     for(j=i+1; (j<nral); j++)
 639       if (aa[i] == aa[j]) {
 640         sprintf(errbuf,"Double atoms in %s (%d)",dir2str(d),aa[i]);
 641         warning(errbuf);
 642       }
 643
 644   /* default force parameters  */
 645   for(j=0; (j<MAXATOMLIST); j++)
 646     param.a[j] = aa[j]-1;
 647   for(j=0; (j<MAXFORCEPARAM); j++)
 648     param.c[j] = 0.0;
 649
 650   /* Get force params for normal and free energy perturbation
 651    * studies, as determined by types!
 652    */
 653   bFoundA = default_params(ftype,bondtype,at,&param,FALSE);
 654   bFoundB = default_params(ftype,bondtype,at,&param,TRUE);
 655   bDef    = TRUE;
 656
 657   nrfp  = NRFP(ftype);
 658   if (nread > nral) {
 659     strcpy(format,asformat[nral-1]);
 660     strcat(format,ccformat[nrfp-1]);
 661
 662     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5]);
 663     if (nread > nrfp) {
 664       warning("Too many parameters");
 665       nread = nrfp;
 666     }
 667
 668     for(j=0; (j<nread); j++)
 669       param.c[j]=cc[j];
 670
 671     /* Check whether we have to use the defaults */
 672     if (nread == nrfp)
 673       bDef = FALSE;
 674   } else
 675     nread = 0;
 676   /* nread now holds the number of force parameters read! */
 677
 678   if (bDef) {
 679     /* Use defaults */
 680     nrfpA=interaction_function[ftype].nrfpA;
 681
 682     if (nread < nrfpA) {
 683       if (!bFoundA) {
 684         if (interaction_function[ftype].flags & IF_DUMMY) {
 685           /* set them to NOTSET, will be calculated later */
 686           for(j=0; (j<MAXFORCEPARAM); j++)
 687             param.c[j] = NOTSET;
 688           if (bSwapParity)
 689             param.C1 = -1; /* flag to swap parity of dummy construction */
 690         } else {
 691           sprintf(errbuf,"No default %s types, using zeroes",
 692                   interaction_function[ftype].longname);
 693           warning(errbuf);
 694         }
 695       } else
 696         if (bSwapParity)
 697           switch(ftype) {
 698           case F_DUMMY3FAD:
 699             param.C0 = 360-param.C0;
 700             break;
 701           case F_DUMMY3OUT:
 702             param.C2 = -param.C2;
 703             break;
 704           }
 705     } else {
 706       /* This implies nread < nrfp, and so perturbed values have not
 707        * been read */
 708       if (!bFoundB) {
 709         /* We only have to issue a warning if these atoms are perturbed! */
 710         bPert = FALSE;
 711         for(j=0; (j<nral); j++)
 712           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
 713
 714         if (bPert) {
 715           sprintf(errbuf,"No default %s types for perturbed atoms, "
 716                   "using normal values",interaction_function[ftype].longname);
 717           warning(errbuf);
 718         }
 719         for(j=nrfpA; (j<nrfp); j++)
 720           param.c[j]=param.c[j-nrfpA];
 721       }
 722     }
 723   }
 724   /* Put the values in the appropriate arrays */
 725   push_bondnow (&bond[ftype],&param);
 726 }
 727
 728 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
 729                   int *nrcopies)
 730 {
 731   int  i,copies;
 732   char type[STRLEN];
 733
 734   *nrcopies=0;
 735   if (sscanf(pline,"%s%d",type,&copies) != 2) {
 736     too_few();
 737     return;
 738   }
 739
 740   /* search moleculename */
 741   for (i=0; ((i<nrmols) && strcasecmp(type,*(mols[i].name))); i++)
 742     ;
 743
 744   if (i<nrmols) {
 745     *nrcopies        = copies;
 746     *whichmol        = i;
 747   }
 748   else {
 749     sprintf(errbuf,"No such moleculetype %s",type);
 750     warning(errbuf);
 751   }
 752 }
 753
 754 void init_block2(t_block2 *b2, int natom)
 755 {
 756   int i;
 757
 758   b2->nr=natom;
 759   snew(b2->nra,b2->nr);
 760   snew(b2->a,b2->nr);
 761   for(i=0; (i<b2->nr); i++)
 762     b2->a[i]=NULL;
 763 }
 764
 765 void done_block2(t_block2 *b2)
 766 {
 767   int i;
 768
 769   if (b2->nr) {
 770     for(i=0; (i<b2->nr); i++)
 771       sfree(b2->a[i]);
 772     sfree(b2->a);
 773     sfree(b2->nra);
 774     b2->nr=0;
 775   }
 776 }
 777
 778 void push_excl(char *line, t_block2 *b2)
 779 {
 780   int  i,j;
 781   int  n;
 782   char base[STRLEN],format[STRLEN];
 783
 784   if (sscanf(line,"%d",&i)==0)
 785     return;
 786
 787   if ((1 <= i) && (i <= b2->nr))
 788     i--;
 789   else {
 790     if (debug)
 791       fprintf(debug,"Unbound atom %d\n",i-1);
 792     return;
 793   }
 794   strcpy(base,"%*d");
 795   do {
 796     strcpy(format,base);
 797     strcat(format,"%d");
 798     n=sscanf(line,format,&j);
 799     if (n == 1) {
 800       if ((1 <= j) && (j <= b2->nr)) {
 801         j--;
 802         srenew(b2->a[i],++(b2->nra[i]));
 803         b2->a[i][b2->nra[i]-1]=j;
 804         /* also add the reverse exclusion! */
 805         srenew(b2->a[j],++(b2->nra[j]));
 806         b2->a[j][b2->nra[j]-1]=i;
 807         strcat(base,"%*d");
 808       }
 809       else
 810         fatal_error(0,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
 811     }
 812   } while (n == 1);
 813 }
 814
 815 void b_to_b2(t_block *b, t_block2 *b2)
 816 {
 817   int     i;
 818   atom_id j,a;
 819
 820   for(i=0; (i<b->nr); i++)
 821     for(j=b->index[i]; (j<b->index[i+1]); j++) {
 822       a=b->a[j];
 823       srenew(b2->a[i],++b2->nra[i]);
 824       b2->a[i][b2->nra[i]-1]=a;
 825     }
 826 }
 827
 828 void b2_to_b(t_block2 *b2, t_block *b)
 829 {
 830   int     i,nra;
 831   atom_id j;
 832
 833   nra=0;
 834   for(i=0; (i<b2->nr); i++) {
 835     b->index[i]=nra;
 836     for(j=0; (j<b2->nra[i]); j++)
 837       b->a[nra+j]=b2->a[i][j];
 838     nra+=b2->nra[i];
 839   }
 840   /* terminate list */
 841   b->index[i]=nra;
 842 }
 843
 844 static int icomp(const void *v1, const void *v2)
 845 {
 846   return (*((atom_id *) v1))-(*((atom_id *) v2));
 847 }
 848
 849 void merge_excl(t_block *excl, t_block2 *b2)
 850 {
 851   int     i,k;
 852   atom_id j;
 853   int     nra;
 854
 855   if (!b2->nr)
 856     return;
 857   else if (b2->nr != excl->nr) {
 858     fatal_error(0,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
 859                 b2->nr,excl->nr);
 860   }
 861   else if (debug)
 862     fprintf(debug,"Entering merge_excl\n");
 863
 864   /* First copy all entries from excl to b2 */
 865   b_to_b2(excl,b2);
 866
 867   /* Count and sort the exclusions */
 868   nra=0;
 869   for(i=0; (i<b2->nr); i++) {
 870     if (b2->nra[i] > 0) {
 871       /* remove double entries */
 872       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
 873       k=1;
 874       for(j=1; (j<b2->nra[i]); j++)
 875         if (b2->a[i][j]!=b2->a[i][k-1]) {
 876           b2->a[i][k]=b2->a[i][j];
 877           k++;
 878         }
 879       b2->nra[i]=k;
 880       nra+=b2->nra[i];
 881     }
 882   }
 883   excl->nra=nra;
 884   srenew(excl->a,excl->nra);
 885
 886   b2_to_b(b2,excl);
 887 }
 888