resolved merge conflicts
[alexxy/gromacs.git] / 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 i,n,k,found,gfound;
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     found = 0;
898     gfound = -1;
899     for(i=0;i<get_atomtype_ntypes(at) && !found;i++)
900     {
901         if(gmx_strncasecmp(atypename,get_atomtype_name(i,at),STRLEN-1)==0)
902         {
903             found = i;
904                         gfound = i;
905                 }
906     }
907         
908         if (gfound==-1)
909     {
910                 printf("Couldn't find topology match for atomtype %s\n",atypename);
911                 abort();
912     }
913     
914         set_atomtype_gbparam(at,found,radius,vol,surftens,gb_radius,S_hct);
915 }
916
917 void 
918 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at, 
919               t_bond_atomtype bat, char *line,
920               warninp_t wi)
921 {
922         const char *formal="%s%s%s%s%s%s%s%s";
923         
924         int      i,j,ft,ftype,nn,nrfp,nrfpA,nrfpB;
925         int      start;
926         int      nxcmap,nycmap,ncmap,read_cmap,sl,nct;
927         char     s[20],alc[MAXATOMLIST+1][20];
928         t_param  p;
929         bool     bAllowRepeat;
930         char     errbuf[256];
931         
932         /* Keep the compiler happy */
933         read_cmap = 0;
934         start     = 0;
935         
936         if((nn=sscanf(line,formal,alc[0],alc[1],alc[2],alc[3],alc[4],alc[5],alc[6],alc[7])) != nral+3)
937         {
938                 sprintf(errbuf,"Incorrect number of atomtypes for cmap (%d instead of 5)",nn-1);
939                 warning_error(wi,errbuf);
940                 return;
941         }
942         
943         /* Compute an offset for each line where the cmap parameters start
944          * ie. where the atom types and grid spacing information ends
945          */
946         for(i=0;i<nn;i++)
947         {
948                 start += (int)strlen(alc[i]);
949         }
950         
951         /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */
952         /* start is the position on the line where we start to read the actual cmap grid data from the itp file */
953         start = start + nn -1; 
954         
955         ft     = strtol(alc[nral],NULL,10);
956         nxcmap = strtol(alc[nral+1],NULL,10);
957         nycmap = strtol(alc[nral+2],NULL,10);
958         
959         /* Check for equal grid spacing in x and y dims */
960         if(nxcmap!=nycmap)
961         {
962                 gmx_fatal(FARGS,"Not the same grid spacing in x and y for cmap grid: x=%d, y=%d",nxcmap,nycmap);
963         }
964         
965         ncmap  = nxcmap*nycmap;
966         ftype = ifunc_index(d,ft);
967         nrfpA = strtol(alc[6],NULL,10)*strtol(alc[6],NULL,10);
968         nrfpB = strtol(alc[7],NULL,10)*strtol(alc[7],NULL,10);
969         nrfp  = nrfpA+nrfpB;
970         
971         /* Allocate memory for the CMAP grid */
972         bt->ncmap+=nrfp;
973         srenew(bt->cmap,bt->ncmap);
974         
975         /* Read in CMAP parameters */
976         sl = 0;
977         for(i=0;i<ncmap;i++)
978         {
979                 nn=sscanf(line+start+sl,"%s",s);
980                 sl+=strlen(s)+1;
981                 bt->cmap[i+(bt->ncmap)-nrfp]=strtod(s,NULL);
982                 
983                 if(nn==1)
984                 {
985                         read_cmap++;
986                 }
987                 else
988                 {
989                         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]);
990                 }
991                 
992         }
993         
994         /* Check do that we got the number of parameters we expected */
995         if(read_cmap==nrfpA)
996         {
997                 for(i=0;i<ncmap;i++)
998                 {
999                         bt->cmap[i+ncmap]=bt->cmap[i];
1000                 }
1001         }
1002         else
1003         {
1004                 if(read_cmap<nrfpA)
1005                 {
1006                         warning_error(wi,"Not enough cmap parameters");
1007                 }
1008                 else if(read_cmap > nrfpA && read_cmap < nrfp)
1009                 {
1010                         warning_error(wi,"Too many cmap parameters or not enough parameters for topology B");
1011                 }
1012                 else if(read_cmap>nrfp)
1013                 {
1014                         warning_error(wi,"Too many cmap parameters");
1015                 }
1016         }
1017         
1018         
1019         /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1020          * so we can safely assign them each time 
1021          */
1022         bt->grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
1023         bt->nc           = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1024         nct              = (nral+1) * bt->nc; 
1025         
1026         /* Allocate memory for the cmap_types information */
1027         srenew(bt->cmap_types,nct);
1028         
1029         for(i=0; (i<nral); i++) 
1030         {
1031                 if(at && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1032                 {
1033                         gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
1034                 }
1035                 else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1036                 {
1037                         gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
1038                 }
1039                 
1040                 /* Assign a grid number to each cmap_type */
1041                 bt->cmap_types[bt->nct++]=get_bond_atomtype_type(alc[i],bat);
1042         }
1043         
1044         /* Assign a type number to this cmap */
1045         bt->cmap_types[bt->nct++]=bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */ 
1046         
1047         /* Check for the correct number of atoms (again) */
1048         if(bt->nct!=nct)
1049         {
1050                 gmx_fatal(FARGS,"Incorrect number of atom types (%d) in cmap type %d\n",nct,bt->nc);
1051         }
1052         
1053         /* Is this correct?? */
1054         for(i=0;i<MAXFORCEPARAM;i++)
1055         {
1056                 p.c[i]=NOTSET;
1057         }
1058         
1059         /* Push the bond to the bondlist */
1060         push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line,wi);
1061 }
1062
1063
1064 static void push_atom_now(t_symtab *symtab,t_atoms *at,int atomnr,
1065                           int atomicnumber,
1066                           int type,char *ctype,int ptype,
1067                           char *resnumberic,int cgnumber,
1068                           char *resname,char *name,real m0,real q0,
1069                           int typeB,char *ctypeB,real mB,real qB)
1070 {
1071   int j,resind=0,resnr;
1072   unsigned char ric;
1073   int nr = at->nr;
1074
1075   if (((nr==0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1076     gmx_fatal(FARGS,"Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)",atomnr,at->nr);
1077
1078   j = strlen(resnumberic) - 1;
1079   if (isdigit(resnumberic[j])) {
1080     ric = ' ';
1081   } else {
1082     ric = resnumberic[j];
1083     if (j == 0 || !isdigit(resnumberic[j-1])) {
1084       gmx_fatal(FARGS,"Invalid residue number '%s' for atom %d",
1085                 resnumberic,atomnr);
1086     }
1087   }
1088   resnr = strtol(resnumberic,NULL,10);
1089
1090   if (nr > 0) {
1091     resind = at->atom[nr-1].resind;
1092   }
1093   if (nr == 0 || strcmp(resname,*at->resinfo[resind].name) != 0 ||
1094       resnr != at->resinfo[resind].nr ||
1095       ric   != at->resinfo[resind].ic) {
1096     if (nr == 0) {
1097       resind = 0;
1098     } else {
1099       resind++;
1100     }
1101     at->nres = resind + 1;
1102     srenew(at->resinfo,at->nres);
1103     at->resinfo[resind].name = put_symtab(symtab,resname);
1104     at->resinfo[resind].nr   = resnr;
1105     at->resinfo[resind].ic   = ric;
1106   } else {
1107     resind = at->atom[at->nr-1].resind;
1108   }
1109
1110   /* New atom instance
1111    * get new space for arrays 
1112    */
1113   srenew(at->atom,nr+1);
1114   srenew(at->atomname,nr+1);
1115   srenew(at->atomtype,nr+1);
1116   srenew(at->atomtypeB,nr+1);
1117
1118   /* fill the list */
1119   at->atom[nr].type  = type;
1120   at->atom[nr].ptype = ptype;
1121   at->atom[nr].q     = q0;
1122   at->atom[nr].m     = m0;
1123   at->atom[nr].typeB = typeB;
1124   at->atom[nr].qB    = qB;
1125   at->atom[nr].mB    = mB;
1126   
1127   at->atom[nr].resind = resind;
1128   at->atom[nr].atomnumber = atomicnumber;
1129   at->atomname[nr] = put_symtab(symtab,name);
1130   at->atomtype[nr] = put_symtab(symtab,ctype);
1131   at->atomtypeB[nr] = put_symtab(symtab,ctypeB);
1132   at->nr++;
1133 }
1134
1135 void push_cg(t_block *block, int *lastindex, int index, int a)
1136 {
1137   if (debug)
1138     fprintf (debug,"Index %d, Atom %d\n",index,a);
1139
1140   if (((block->nr) && (*lastindex != index)) || (!block->nr)) {
1141     /* add a new block */
1142     block->nr++;
1143     srenew(block->index,block->nr+1);
1144   }
1145   block->index[block->nr] = a + 1;
1146   *lastindex = index;
1147 }
1148
1149 void push_atom(t_symtab *symtab,t_block *cgs,
1150                t_atoms *at,gpp_atomtype_t atype,char *line,int *lastcg,
1151                warninp_t wi)
1152 {
1153   int           nr,ptype;
1154   int           cgnumber,atomnr,type,typeB,nscan;
1155   char          id[STRLEN],ctype[STRLEN],ctypeB[STRLEN],
1156                 resnumberic[STRLEN],resname[STRLEN],name[STRLEN],check[STRLEN];
1157   double        m,q,mb,qb;
1158   real          m0,q0,mB,qB;
1159
1160   /* Make a shortcut for writing in this molecule  */
1161   nr = at->nr;
1162
1163   /* Fixed parameters */
1164   if (sscanf(line,"%s%s%s%s%s%d",
1165              id,ctype,resnumberic,resname,name,&cgnumber) != 6) {
1166     too_few(wi);
1167     return;
1168   }
1169   sscanf(id,"%d",&atomnr);
1170   if ((type  = get_atomtype_type(ctype,atype)) == NOTSET)
1171     gmx_fatal(FARGS,"Atomtype %s not found",ctype);
1172   ptype = get_atomtype_ptype(type,atype);
1173
1174   /* Set default from type */  
1175   q0    = get_atomtype_qA(type,atype);
1176   m0    = get_atomtype_massA(type,atype);
1177   typeB = type;
1178   qB    = q0;
1179   mB    = m0;
1180   
1181   /* Optional parameters */
1182   nscan = sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1183                  &q,&m,ctypeB,&qb,&mb,check);
1184   
1185   /* Nasty switch that falls thru all the way down! */
1186   if (nscan > 0) {
1187     q0 = qB = q;
1188     if (nscan > 1) {
1189       m0 = mB = m;
1190       if (nscan > 2) {
1191         if ((typeB = get_atomtype_type(ctypeB,atype)) == NOTSET) {
1192           gmx_fatal(FARGS,"Atomtype %s not found",ctypeB);
1193         }
1194         qB = get_atomtype_qA(typeB,atype);
1195         mB = get_atomtype_massA(typeB,atype);
1196         if (nscan > 3) {
1197           qB = qb;
1198           if (nscan > 4) {
1199             mB = mb;
1200             if (nscan > 5)
1201             warning_error(wi,"Too many parameters");
1202           }
1203         }
1204       }
1205     }
1206   }
1207   if (debug) 
1208     fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
1209   
1210   push_cg(cgs,lastcg,cgnumber,nr);
1211
1212   push_atom_now(symtab,at,atomnr,get_atomtype_atomnumber(type,atype),
1213                 type,ctype,ptype,resnumberic,cgnumber,
1214                 resname,name,m0,q0,typeB,
1215                 typeB==type ? ctype : ctypeB,mB,qB);
1216 }
1217
1218 void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line,
1219                warninp_t wi)
1220 {
1221   char type[STRLEN];
1222   int nrexcl,i;
1223   t_molinfo *newmol;
1224
1225   if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
1226       warning_error(wi,"Expected a molecule type name and nrexcl");
1227   }
1228   
1229   /* Test if this atomtype overwrites another */
1230   i = 0;
1231   while (i < *nmol) {
1232     if (strcasecmp(*((*mol)[i].name),type) == 0)
1233       gmx_fatal(FARGS,"moleculetype %s is redefined",type);
1234     i++;
1235   }
1236   
1237   (*nmol)++;
1238   srenew(*mol,*nmol);
1239   newmol = &((*mol)[*nmol-1]);
1240   init_molinfo(newmol);
1241   
1242   /* Fill in the values */
1243   newmol->name     = put_symtab(symtab,type);
1244   newmol->nrexcl   = nrexcl;
1245   newmol->excl_set = FALSE;
1246 }
1247
1248 static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
1249                               t_param *p,int c_start,bool bB,bool bGenPairs)
1250 {
1251   int      i,j,ti,tj,ntype;
1252   bool     bFound;
1253   t_param  *pi=NULL;
1254   int      nr   = bt[ftype].nr;
1255   int      nral = NRAL(ftype);
1256   int      nrfp = interaction_function[ftype].nrfpA;
1257   int      nrfpB= interaction_function[ftype].nrfpB;
1258
1259   if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1260     return TRUE;
1261
1262   bFound=FALSE;
1263   if(bGenPairs) {
1264     /* First test the generated-pair position to save
1265      * time when we have 1000*1000 entries for e.g. OPLS...
1266      */
1267     ntype=sqrt(nr);
1268     if (bB) {
1269       ti=at->atom[p->a[0]].typeB;
1270       tj=at->atom[p->a[1]].typeB;
1271     } else {
1272       ti=at->atom[p->a[0]].type;
1273       tj=at->atom[p->a[1]].type;
1274     } 
1275     pi=&(bt[ftype].param[ntype*ti+tj]);
1276     bFound=((ti == pi->a[0]) && (tj == pi->a[1]));
1277   }
1278
1279   /* Search explicitly if we didnt find it */
1280   if(!bFound) {
1281     for (i=0; ((i < nr) && !bFound); i++) {
1282       pi=&(bt[ftype].param[i]);
1283       if (bB)
1284         for (j=0; ((j < nral) && 
1285                    (at->atom[p->a[j]].typeB == pi->a[j])); j++);
1286       else
1287         for (j=0; ((j < nral) && 
1288                    (at->atom[p->a[j]].type == pi->a[j])); j++);
1289       bFound=(j==nral);
1290     }
1291   }
1292   
1293   if (bFound) {
1294     if (bB) {
1295       if (nrfp+nrfpB > MAXFORCEPARAM) {
1296         gmx_incons("Too many force parameters");
1297       }
1298       for (j=c_start; (j < nrfpB); j++)
1299         p->c[nrfp+j] = pi->c[j];
1300     }
1301     else
1302       for (j=c_start; (j < nrfp); j++)
1303         p->c[j] = pi->c[j];
1304   }
1305   else {
1306     for (j=c_start; (j < nrfp); j++)
1307       p->c[j] = 0.0;
1308   }
1309   return bFound;
1310 }
1311
1312 static bool default_cmap_params(int ftype, t_params bondtype[],
1313                                 t_atoms *at, gpp_atomtype_t atype,
1314                                 t_param *p, bool bB,
1315                                 int *cmap_type, int *nparam_def)
1316 {
1317         int i,j,nparam_found;
1318         int ct;
1319         bool bFound=FALSE;
1320         
1321         nparam_found = 0;
1322         ct           = 0;
1323         
1324         /* Match the current cmap angle against the list of cmap_types */
1325         for(i=0;i<bondtype->nct && !bFound;i+=6)
1326         {
1327                 if(bB)
1328                 {
1329                         
1330                 }
1331                 else
1332                 {
1333                         if( 
1334                            (get_atomtype_batype(at->atom[p->a[0]].type,atype)==bondtype->cmap_types[i])   &&
1335                            (get_atomtype_batype(at->atom[p->a[1]].type,atype)==bondtype->cmap_types[i+1]) &&
1336                            (get_atomtype_batype(at->atom[p->a[2]].type,atype)==bondtype->cmap_types[i+2]) &&
1337                            (get_atomtype_batype(at->atom[p->a[3]].type,atype)==bondtype->cmap_types[i+3]) &&
1338                            (get_atomtype_batype(at->atom[p->a[4]].type,atype)==bondtype->cmap_types[i+4]))
1339                         {
1340                                 /* Found cmap torsion */
1341                                 bFound=TRUE;
1342                                 ct = bondtype->cmap_types[i+5];
1343                                 nparam_found=1;
1344                         }
1345                 }
1346         }
1347         
1348         /* If we did not find a matching type for this cmap torsion */
1349         if(!bFound)
1350         {
1351                 gmx_fatal(FARGS,"Unknown cmap torsion between atoms %d %d %d %d %d\n",
1352                                   p->a[0]+1,p->a[1]+1,p->a[2]+1,p->a[3]+1,p->a[4]+1);
1353         }
1354         
1355         *nparam_def = nparam_found;
1356         *cmap_type  = ct;
1357         
1358         return bFound;
1359 }
1360
1361 static bool default_params(int ftype,t_params bt[],
1362                            t_atoms *at,gpp_atomtype_t atype,
1363                            t_param *p,bool bB,
1364                            t_param **param_def,
1365                int *nparam_def)
1366 {
1367   int      i,j,nparam_found;
1368   bool     bFound,bSame;
1369   t_param  *pi=NULL;
1370   t_param  *pj=NULL;
1371   int      nr   = bt[ftype].nr;
1372   int      nral = NRAL(ftype);
1373   int      nrfpA= interaction_function[ftype].nrfpA;
1374   int      nrfpB= interaction_function[ftype].nrfpB;
1375
1376   if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1377     return TRUE;
1378
1379         
1380   /* We allow wildcards now. The first type (with or without wildcards) that
1381    * fits is used, so you should probably put the wildcarded bondtypes
1382    * at the end of each section.
1383    */
1384   bFound=FALSE;
1385   nparam_found=0;
1386   /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a 
1387    * special case for this. Check for B state outside loop to speed it up.
1388    */
1389   if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS) 
1390   {
1391           if (bB) 
1392           {
1393                   for (i=0; ((i < nr) && !bFound); i++) 
1394                   {
1395                           pi=&(bt[ftype].param[i]);
1396                           bFound=
1397                           (
1398                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1399                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1400                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1401                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1402                            );
1403                   }
1404           } 
1405           else 
1406           {
1407                   /* State A */
1408                   for (i=0; ((i < nr) && !bFound); i++) 
1409                   {
1410                           pi=&(bt[ftype].param[i]);
1411                           bFound=
1412                           (
1413                            ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].type,atype)==pi->AI)) &&         
1414                            ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].type,atype)==pi->AJ)) &&
1415                            ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].type,atype)==pi->AK)) &&
1416                            ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].type,atype)==pi->AL))
1417                            );
1418                   }
1419           }
1420           /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1421            * The rules in that case is that additional matches HAVE to be on adjacent lines!
1422            */
1423           if(bFound==TRUE)
1424           {
1425                   nparam_found++;
1426                   bSame = TRUE;
1427                   /* Continue from current i value */
1428                   for(j=i+1 ; j<nr && bSame ; j+=2)
1429                   {
1430                           pj=&(bt[ftype].param[j]);
1431                           bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1432                           if(bSame)
1433                                   nparam_found++;
1434                           /* nparam_found will be increased as long as the numbers match */
1435                   }
1436           }
1437   } else { /* Not a dihedral */
1438           for (i=0; ((i < nr) && !bFound); i++) {
1439                   pi=&(bt[ftype].param[i]);
1440                   if (bB)
1441                           for (j=0; ((j < nral) && 
1442                                                  (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
1443                                   ;
1444                   else
1445                           for (j=0; ((j < nral) && 
1446                                                  (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
1447                                   ;
1448                   bFound=(j==nral);
1449           }
1450           if(bFound)
1451                   nparam_found=1;
1452   }
1453
1454   *param_def = pi;
1455   *nparam_def = nparam_found;
1456           
1457   return bFound;
1458 }
1459
1460
1461
1462 void push_bond(directive d,t_params bondtype[],t_params bond[],
1463                t_atoms *at,gpp_atomtype_t atype,char *line,
1464                bool bBonded,bool bGenPairs,real fudgeQQ,
1465                bool bZero,bool *bWarn_copy_A_B,
1466                warninp_t wi)
1467 {
1468   const char *aaformat[MAXATOMLIST]= {
1469     "%d%d",
1470     "%d%d%d",
1471     "%d%d%d%d",
1472     "%d%d%d%d%d",
1473     "%d%d%d%d%d%d",
1474     "%d%d%d%d%d%d%d"
1475   };
1476   const char *asformat[MAXATOMLIST]= {
1477     "%*s%*s",
1478     "%*s%*s%*s",
1479     "%*s%*s%*s%*s",
1480     "%*s%*s%*s%*s%*s",
1481     "%*s%*s%*s%*s%*s%*s",
1482     "%*s%*s%*s%*s%*s%*s%*s"
1483   };
1484   const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1485   int      nr,i,j,nral,nread,ftype;
1486   char     format[STRLEN];
1487   /* One force parameter more, so we can check if we read too many */
1488   double   cc[MAXFORCEPARAM+1];
1489   int      aa[MAXATOMLIST+1];
1490   t_param  param,paramB,*param_defA,*param_defB;
1491   bool     bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
1492   int      nparam_defA,nparam_defB;
1493   char  errbuf[256];
1494
1495   nparam_defA=nparam_defB=0;
1496         
1497   ftype = ifunc_index(d,1);
1498   nral  = NRAL(ftype);
1499   for(j=0; j<MAXATOMLIST; j++)
1500     aa[j]=NOTSET;
1501   bDef = (NRFP(ftype) > 0);
1502   
1503   nread = sscanf(line,aaformat[nral-1],
1504                  &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1505   if (nread < nral) {
1506     too_few(wi);
1507     return;
1508   } else if (nread == nral) 
1509     ftype = ifunc_index(d,1);
1510   else {
1511     /* this is a hack to allow for virtual sites with swapped parity */
1512     bSwapParity = (aa[nral]<0);
1513     if (bSwapParity)
1514       aa[nral] = -aa[nral];
1515     ftype = ifunc_index(d,aa[nral]);
1516     if (bSwapParity)
1517       switch(ftype) {
1518       case F_VSITE3FAD:
1519       case F_VSITE3OUT:
1520         break;
1521       default:
1522         gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1523                     interaction_function[F_VSITE3FAD].longname,
1524                     interaction_function[F_VSITE3OUT].longname);
1525       }
1526   }
1527   
1528   /* Check for double atoms and atoms out of bounds */
1529   for(i=0; (i<nral); i++) {
1530     if ( aa[i] < 1 || aa[i] > at->nr )
1531       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1532                 "Atom index (%d) in %s out of bounds (1-%d).\n"
1533                 "This probably means that you have inserted topology section \"%s\"\n"
1534                 "in a part belonging to a different molecule than you intended to.\n" 
1535                 "In that case move the \"%s\" section to the right molecule.",
1536                 get_warning_file(wi),get_warning_line(wi),
1537                 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1538     for(j=i+1; (j<nral); j++)
1539       if (aa[i] == aa[j]) {
1540         sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1541         warning(wi,errbuf);
1542       }
1543   }
1544   if (ftype == F_SETTLE)
1545     if (aa[0]+2 > at->nr)
1546       gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1547                   "             Atom index (%d) in %s out of bounds (1-%d)\n"
1548                   "             Settle works on atoms %d, %d and %d",
1549                   get_warning_file(wi),get_warning_line(wi),
1550                   aa[0],dir2str(d),at->nr,
1551                   aa[0],aa[0]+1,aa[0]+2);
1552   
1553   /* default force parameters  */
1554   for(j=0; (j<MAXATOMLIST); j++)
1555     param.a[j] = aa[j]-1;
1556   for(j=0; (j<MAXFORCEPARAM); j++)
1557     param.c[j] = 0.0;
1558   
1559   /* Get force params for normal and free energy perturbation
1560    * studies, as determined by types!
1561    */
1562   if (bBonded) {
1563     bFoundA = default_params(ftype,bondtype,at,atype,&param,FALSE,&param_defA,&nparam_defA);
1564     if (bFoundA) {
1565       /* Copy the A-state and B-state default parameters */
1566       for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1567         param.c[j] = param_defA->c[j];
1568     }
1569     bFoundB = default_params(ftype,bondtype,at,atype,&param,TRUE,&param_defB,&nparam_defB);
1570     if (bFoundB) {
1571       /* Copy only the B-state default parameters */
1572       for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1573         param.c[j] = param_defB->c[j];
1574     }
1575   } else if (ftype == F_LJ14) {
1576     bFoundA = default_nb_params(ftype, bondtype,at,&param,0,FALSE,bGenPairs);
1577     bFoundB = default_nb_params(ftype, bondtype,at,&param,0,TRUE, bGenPairs);
1578   } else if (ftype == F_LJC14_Q) {
1579     param.c[0] = fudgeQQ;
1580     /* Fill in the A-state charges as default parameters */
1581     param.c[1] = at->atom[param.a[0]].q;
1582     param.c[2] = at->atom[param.a[1]].q;
1583     /* The default LJ parameters are the standard 1-4 parameters */
1584     bFoundA = default_nb_params(F_LJ14,bondtype,at,&param,3,FALSE,bGenPairs);
1585     bFoundB = TRUE;
1586   } else if (ftype == F_LJC_PAIRS_NB) {
1587     /* Defaults are not supported here */
1588     bFoundA = FALSE;
1589     bFoundB = TRUE;
1590   } else {
1591     gmx_incons("Unknown function type in push_bond");
1592   }
1593         
1594   if (nread > nral) {  
1595           /* Manually specified parameters - in this case we discard multiple torsion info! */
1596           
1597     strcpy(format,asformat[nral-1]);
1598     strcat(format,ccformat);
1599     
1600     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1601                    &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1602
1603     if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1604       {
1605         /* We only have to issue a warning if these atoms are perturbed! */
1606         bPert = FALSE;
1607         for(j=0; (j<nral); j++)
1608           bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1609
1610         if (bPert && *bWarn_copy_A_B)
1611           {
1612             sprintf(errbuf,
1613                     "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1614             warning(wi,errbuf);
1615             *bWarn_copy_A_B = FALSE;
1616           }
1617         
1618         /* If only the A parameters were specified, copy them to the B state */
1619         /* The B-state parameters correspond to the first nrfpB
1620          * A-state parameters.
1621          */
1622         for(j=0; (j<NRFPB(ftype)); j++)
1623           cc[nread++] = cc[j];
1624       }
1625     
1626     /* If nread was 0 or EOF, no parameters were read => use defaults.
1627      * If nread was nrfpA we copied above so nread=nrfp.
1628      * If nread was nrfp we are cool.
1629      * For F_LJC14_Q we allow supplying fudgeQQ only.
1630      * Anything else is an error!
1631      */ 
1632     if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1633         !(ftype == F_LJC14_Q && nread == 1))
1634       {
1635         gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1636                   nread,NRFPA(ftype),NRFP(ftype),
1637                   interaction_function[ftype].longname);        
1638       }
1639     
1640     for(j=0; (j<nread); j++)
1641       param.c[j]=cc[j];
1642       
1643     /* Check whether we have to use the defaults */
1644     if (nread == NRFP(ftype))
1645       bDef = FALSE;
1646   } 
1647   else
1648   {
1649     nread = 0;
1650   }
1651   /* nread now holds the number of force parameters read! */
1652  
1653         
1654         
1655         
1656         if (bDef) 
1657         {
1658           /* Use defaults */
1659                 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1660                 if(ftype==F_PDIHS)
1661                 {
1662                         if ((nparam_defA != nparam_defB) || ((nparam_defA>1 || nparam_defB>1) && (param_defA!=param_defB)))
1663                         {
1664                                 sprintf(errbuf,
1665                                                 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1666                                                 "Please specify perturbed parameters manually for this torsion in your topology!");
1667                                 warning_error(wi,errbuf);
1668                         }
1669                 }
1670                 
1671           if (nread > 0 && nread < NRFPA(ftype)) 
1672           {
1673                   /* Issue an error, do not use defaults */
1674                   sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1675                   warning_error(wi,errbuf);
1676           } 
1677                   
1678       if (nread == 0 || nread == EOF) 
1679           {
1680                   if (!bFoundA) 
1681                   {
1682                           if (interaction_function[ftype].flags & IF_VSITE) 
1683                           {
1684                                   /* set them to NOTSET, will be calculated later */
1685                                   for(j=0; (j<MAXFORCEPARAM); j++)
1686                                           param.c[j] = NOTSET;
1687                                   
1688                                   if (bSwapParity)
1689                                           param.C1 = -1; /* flag to swap parity of vsite construction */
1690                           } 
1691                           else 
1692                           {
1693                                   if (bZero) 
1694                                   {
1695                                           fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1696                                                           interaction_function[ftype].longname);
1697                                   }
1698                                   else 
1699                                   {
1700                                           sprintf(errbuf,"No default %s types",interaction_function[ftype].longname);
1701                                           warning_error(wi,errbuf);
1702                                   }
1703                           }
1704                   }
1705                   else
1706                   {
1707                           if (bSwapParity)
1708                           {
1709                                   switch(ftype) 
1710                                   {
1711                                           case F_VSITE3FAD:
1712                                                   param.C0 = 360 - param.C0;
1713                                                   break;
1714                                           case F_VSITE3OUT:
1715                                                   param.C2 = -param.C2;
1716                                                   break;
1717                                   }
1718                           }
1719                   }
1720                   if (!bFoundB) 
1721                   {
1722                           /* We only have to issue a warning if these atoms are perturbed! */
1723                           bPert = FALSE;
1724                           for(j=0; (j<nral); j++)
1725                                   bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1726                           
1727                           if (bPert)
1728                           {
1729                                   sprintf(errbuf,"No default %s types for perturbed atoms, "
1730                                                   "using normal values",interaction_function[ftype].longname);
1731                                   warning(wi,errbuf);
1732                           }
1733                   }
1734           }
1735   }
1736   
1737         if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1738                 && param.c[5]!=param.c[2])
1739                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1740                                   "             %s multiplicity can not be perturbed %f!=%f",
1741                                   get_warning_file(wi),get_warning_line(wi),
1742                                   interaction_function[ftype].longname,
1743                                   param.c[2],param.c[5]);
1744         
1745         if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1746                 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1747                                   "             %s table number can not be perturbed %d!=%d",
1748                                   get_warning_file(wi),get_warning_line(wi),
1749                                   interaction_function[ftype].longname,
1750                                   (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1751         
1752         /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1753         if (ftype==F_RBDIHS) {
1754                 nr=0;
1755                 for(i=0;i<NRFP(ftype);i++) {
1756                         if(param.c[i]!=0)
1757                                 nr++;
1758                 }
1759                 if(nr==0)
1760                         return;
1761         }
1762         
1763         /* Put the values in the appropriate arrays */
1764         add_param_to_list (&bond[ftype],&param);
1765
1766         /* Push additional torsions from FF for ftype==9 if we have them.
1767          * We have already checked that the A/B states do not differ in this case,
1768          * so we do not have to double-check that again, or the vsite stuff.
1769          * In addition, those torsions cannot be automatically perturbed.
1770          */
1771         if(bDef && ftype==F_PDIHS)
1772         {
1773                 for(i=1;i<nparam_defA;i++)
1774                 {
1775                         /* Advance pointer! */
1776                         param_defA+=2; 
1777                         for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1778                                 param.c[j] = param_defA->c[j];
1779                         /* And push the next term for this torsion */
1780                         add_param_to_list (&bond[ftype],&param);                
1781                 }
1782         }
1783 }
1784
1785 void push_cmap(directive d, t_params bondtype[], t_params bond[],
1786                            t_atoms *at, gpp_atomtype_t atype, char *line,
1787                            bool *bWarn_copy_A_B,
1788                warninp_t wi)
1789 {
1790         const char *aaformat[MAXATOMLIST+1]= 
1791         {
1792                 "%d",
1793                 "%d%d",
1794                 "%d%d%d",
1795                 "%d%d%d%d",
1796                 "%d%d%d%d%d",
1797                 "%d%d%d%d%d%d",
1798                 "%d%d%d%d%d%d%d"
1799         };
1800         
1801         int i,j,nr,ftype,nral,nread,ncmap_params;
1802         int cmap_type;
1803         int aa[MAXATOMLIST+1];
1804         char  errbuf[256];
1805         bool bFound;
1806         t_param param,paramB,*param_defA,*param_defB;
1807         
1808         ftype        = ifunc_index(d,1);
1809         nral         = NRAL(ftype);
1810         nr           = bondtype[ftype].nr;
1811         ncmap_params = 0;
1812         
1813         nread = sscanf(line,aaformat[nral-1],
1814                                    &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1815         
1816         if (nread < nral) 
1817         {
1818                 too_few(wi);
1819                 return;
1820         } 
1821         else if (nread == nral) 
1822         {
1823                 ftype = ifunc_index(d,1);
1824         }
1825         
1826         /* Check for double atoms and atoms out of bounds */
1827         for(i=0;i<nral;i++)
1828         {
1829                 if(aa[i] < 1 || aa[i] > at->nr)
1830                 {
1831                         gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1832                                           "Atom index (%d) in %s out of bounds (1-%d).\n"
1833                                           "This probably means that you have inserted topology section \"%s\"\n"
1834                                           "in a part belonging to a different molecule than you intended to.\n" 
1835                                           "In that case move the \"%s\" section to the right molecule.",
1836                                           get_warning_file(wi),get_warning_line(wi),
1837                                           aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1838                 }
1839                 
1840                 for(j=i+1; (j<nral); j++)
1841                 {
1842                         if (aa[i] == aa[j]) 
1843                         {
1844                                 sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1845                                 warning(wi,errbuf);
1846                         }
1847                 }
1848         }
1849         
1850         /* default force parameters  */
1851         for(j=0; (j<MAXATOMLIST); j++)
1852         {
1853                 param.a[j] = aa[j]-1;
1854         }
1855         for(j=0; (j<MAXFORCEPARAM); j++)
1856         {
1857                 param.c[j] = 0.0;
1858         }       
1859         
1860         /* Get the cmap type for this cmap angle */
1861         bFound = default_cmap_params(ftype,bondtype,at,atype,&param,FALSE,&cmap_type,&ncmap_params);
1862         
1863         /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
1864         if(bFound && ncmap_params==1)
1865         {
1866                 /* Put the values in the appropriate arrays */
1867                 param.c[0]=cmap_type;
1868                 add_param_to_list(&bond[ftype],&param);
1869         }
1870         else
1871         {
1872                 /* This is essentially the same check as in default_cmap_params() done one more time */
1873                 gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
1874                                   param.a[0]+1,param.a[1]+1,param.a[2]+1,param.a[3]+1,param.a[4]+1);
1875         }
1876 }
1877
1878
1879         
1880 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1881                   t_atoms *at,gpp_atomtype_t atype,char *line,
1882                   warninp_t wi)
1883 {
1884   char *ptr;
1885   int  type,ftype,j,n,ret,nj,a;
1886   int  *atc=NULL;
1887   double *weight=NULL,weight_tot;
1888   t_param param;
1889
1890   /* default force parameters  */
1891   for(j=0; (j<MAXATOMLIST); j++)
1892     param.a[j] = NOTSET;
1893   for(j=0; (j<MAXFORCEPARAM); j++)
1894     param.c[j] = 0.0;
1895
1896   ptr = line;
1897   ret = sscanf(ptr,"%d%n",&a,&n);
1898   ptr += n;
1899   if (ret == 0)
1900     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1901               "             Expected an atom index in section \"%s\"",
1902               get_warning_file(wi),get_warning_line(wi),
1903               dir2str(d));
1904   
1905   param.a[0] = a - 1;
1906
1907   ret = sscanf(ptr,"%d%n",&type,&n);
1908   ptr += n;
1909   ftype = ifunc_index(d,type);
1910
1911   weight_tot = 0;
1912   nj = 0;
1913   do {
1914     ret = sscanf(ptr,"%d%n",&a,&n);
1915     ptr += n;
1916     if (ret > 0) {
1917       if (nj % 20 == 0) {
1918         srenew(atc,nj+20);
1919         srenew(weight,nj+20);
1920       }
1921       atc[nj] = a - 1;
1922       switch (type) {
1923       case 1:
1924         weight[nj] = 1;
1925         break;
1926       case 2:
1927         /* Here we use the A-state mass as a parameter.
1928          * Note that the B-state mass has no influence.
1929          */
1930         weight[nj] = at->atom[atc[nj]].m;
1931         break;
1932       case 3:
1933         weight[nj] = -1;
1934         ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1935         ptr += n;
1936         if (weight[nj] < 0)
1937           gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1938                 "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1939                 get_warning_file(wi),get_warning_line(wi),
1940                 nj+1,atc[nj]+1);
1941         break;
1942       default:
1943         gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1944       }
1945       weight_tot += weight[nj];
1946       nj++;
1947     }
1948   } while (ret > 0);
1949
1950   if (nj == 0)
1951     gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1952               "             Expected more than one atom index in section \"%s\"",
1953               get_warning_file(wi),get_warning_line(wi),
1954               dir2str(d));
1955
1956   if (weight_tot == 0)
1957      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1958                "             The total mass of the construting atoms is zero",
1959                get_warning_file(wi),get_warning_line(wi));
1960
1961   for(j=0; j<nj; j++) {
1962     param.a[1] = atc[j];
1963     param.c[0] = nj;
1964     param.c[1] = weight[j]/weight_tot;
1965     /* Put the values in the appropriate arrays */
1966     add_param_to_list (&bond[ftype],&param);
1967   }
1968
1969   sfree(atc);
1970   sfree(weight);
1971 }
1972
1973 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1974               int *nrcopies,
1975               warninp_t wi)
1976 {
1977   int  i,copies;
1978   char type[STRLEN];
1979
1980   *nrcopies=0;
1981   if (sscanf(pline,"%s%d",type,&copies) != 2) {
1982     too_few(wi);
1983     return;
1984   }
1985   
1986   /* search moleculename */
1987   for (i=0; ((i<nrmols) && strcasecmp(type,*(mols[i].name))); i++)
1988     ;
1989
1990   if (i<nrmols) {
1991     *nrcopies        = copies;
1992     *whichmol        = i;
1993   } else
1994     gmx_fatal(FARGS,"No such moleculetype %s",type);
1995 }
1996
1997 void init_block2(t_block2 *b2, int natom)
1998 {
1999   int i;
2000
2001   b2->nr=natom;
2002   snew(b2->nra,b2->nr);
2003   snew(b2->a,b2->nr);
2004   for(i=0; (i<b2->nr); i++)
2005     b2->a[i]=NULL;
2006 }
2007
2008 void done_block2(t_block2 *b2)
2009 {
2010   int i;
2011   
2012   if (b2->nr) {
2013     for(i=0; (i<b2->nr); i++)
2014       sfree(b2->a[i]);
2015     sfree(b2->a);
2016     sfree(b2->nra);
2017     b2->nr=0;
2018   }
2019 }
2020
2021 void push_excl(char *line, t_block2 *b2)
2022 {
2023   int  i,j;
2024   int  n;
2025   char base[STRLEN],format[STRLEN];
2026
2027   if (sscanf(line,"%d",&i)==0)
2028     return;
2029     
2030   if ((1 <= i) && (i <= b2->nr))
2031     i--;
2032   else {
2033     if (debug)
2034       fprintf(debug,"Unbound atom %d\n",i-1);
2035     return;
2036   }
2037   strcpy(base,"%*d");
2038   do {
2039     strcpy(format,base);
2040     strcat(format,"%d");
2041     n=sscanf(line,format,&j);
2042     if (n == 1) {
2043       if ((1 <= j) && (j <= b2->nr)) {
2044         j--;
2045         srenew(b2->a[i],++(b2->nra[i]));
2046         b2->a[i][b2->nra[i]-1]=j;
2047         /* also add the reverse exclusion! */
2048         srenew(b2->a[j],++(b2->nra[j]));
2049         b2->a[j][b2->nra[j]-1]=i;
2050         strcat(base,"%*d");
2051       }
2052       else 
2053         gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
2054     }
2055   } while (n == 1);
2056 }
2057
2058 void b_to_b2(t_blocka *b, t_block2 *b2)
2059 {
2060   int     i;
2061   atom_id j,a;
2062
2063   for(i=0; (i<b->nr); i++)
2064     for(j=b->index[i]; (j<b->index[i+1]); j++) {
2065       a=b->a[j];
2066       srenew(b2->a[i],++b2->nra[i]);
2067       b2->a[i][b2->nra[i]-1]=a;
2068     }
2069 }
2070
2071 void b2_to_b(t_block2 *b2, t_blocka *b)
2072 {
2073   int     i,nra;
2074   atom_id j;
2075
2076   nra=0;
2077   for(i=0; (i<b2->nr); i++) {
2078     b->index[i]=nra;
2079     for(j=0; (j<b2->nra[i]); j++)
2080       b->a[nra+j]=b2->a[i][j];
2081     nra+=b2->nra[i];
2082   }
2083   /* terminate list */
2084   b->index[i]=nra;
2085 }
2086
2087 static int icomp(const void *v1, const void *v2)
2088 {
2089   return (*((atom_id *) v1))-(*((atom_id *) v2));
2090 }
2091
2092 void merge_excl(t_blocka *excl, t_block2 *b2)
2093 {
2094   int     i,k;
2095   atom_id j;
2096   int     nra;
2097
2098   if (!b2->nr)
2099     return;
2100   else if (b2->nr != excl->nr) {
2101     gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2102                 b2->nr,excl->nr);
2103   }
2104   else if (debug)
2105     fprintf(debug,"Entering merge_excl\n");
2106
2107   /* First copy all entries from excl to b2 */
2108   b_to_b2(excl,b2);
2109
2110   /* Count and sort the exclusions */
2111   nra=0;
2112   for(i=0; (i<b2->nr); i++) {
2113     if (b2->nra[i] > 0) {
2114       /* remove double entries */
2115       qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
2116       k=1;
2117       for(j=1; (j<b2->nra[i]); j++)
2118         if (b2->a[i][j]!=b2->a[i][k-1]) {
2119           b2->a[i][k]=b2->a[i][j];
2120           k++;
2121         }
2122       b2->nra[i]=k;
2123       nra+=b2->nra[i];
2124     }
2125   }
2126   excl->nra=nra;
2127   srenew(excl->a,excl->nra);
2128
2129   b2_to_b(b2,excl);
2130 }
2131
2132 int add_atomtype_decoupled(t_symtab *symtab,gpp_atomtype_t at,
2133                            t_nbparam ***nbparam,t_nbparam ***pair)
2134 {
2135   t_atom  atom;
2136   t_param param;
2137   int     i,nr;
2138
2139   /* Define an atom type with all parameters set to zero (no interactions) */
2140   atom.q = 0.0;
2141   atom.m = 0.0;
2142   /* Type for decoupled atoms could be anything,
2143    * this should be changed automatically later when required.
2144    */
2145   atom.ptype = eptAtom;
2146   for (i=0; (i<MAXFORCEPARAM); i++)
2147     param.c[i] = 0.0;
2148
2149   nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0,0,0);
2150
2151   /* Add space in the non-bonded parameters matrix */
2152   realloc_nb_params(at,nbparam,pair);
2153
2154   return nr;
2155 }
2156
2157 static void convert_pairs_to_pairsQ(t_params *plist,
2158                                     real fudgeQQ,t_atoms *atoms)
2159 {
2160   t_param *param;
2161   int  i;
2162   real v,w;
2163
2164   /* Copy the pair list to the pairQ list */
2165   plist[F_LJC14_Q] = plist[F_LJ14];
2166   /* Empty the pair list */
2167   plist[F_LJ14].nr    = 0;
2168   plist[F_LJ14].param = NULL;
2169   param = plist[F_LJC14_Q].param;
2170   for(i=0; i<plist[F_LJC14_Q].nr; i++) {
2171     v = param[i].c[0];
2172     w = param[i].c[1];
2173     param[i].c[0] = fudgeQQ;
2174     param[i].c[1] = atoms->atom[param[i].a[0]].q;
2175     param[i].c[2] = atoms->atom[param[i].a[1]].q;
2176     param[i].c[3] = v;
2177     param[i].c[4] = w;
2178   }
2179 }
2180
2181 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
2182 {
2183   int n,ntype,i,j,k;
2184   t_atom *atom;
2185   t_blocka *excl;
2186   bool bExcl;
2187   t_param param;
2188
2189   n = mol->atoms.nr;
2190   atom = mol->atoms.atom;
2191   
2192   ntype = sqrt(nbp->nr);
2193
2194   for (i=0; i<MAXATOMLIST; i++) 
2195     param.a[i] = NOTSET;
2196   for (i=0; i<MAXFORCEPARAM; i++) 
2197     param.c[i] = NOTSET;
2198
2199   /* Add a pair interaction for all non-excluded atom pairs */
2200   excl = &mol->excls;
2201   for(i=0; i<n; i++) {
2202     for(j=i+1; j<n; j++) {
2203       bExcl = FALSE;
2204       for(k=excl->index[i]; k<excl->index[i+1]; k++) {
2205         if (excl->a[k] == j)
2206           bExcl = TRUE;
2207       }
2208       if (!bExcl) {
2209         if (nb_funct != F_LJ)
2210           gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2211         param.a[0] = i;
2212         param.a[1] = j;
2213         param.c[0] = atom[i].q;
2214         param.c[1] = atom[j].q;
2215         param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2216         param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2217         add_param_to_list(&mol->plist[F_LJC_PAIRS_NB],&param);
2218       }
2219     }
2220   }
2221 }
2222
2223 static void set_excl_all(t_blocka *excl)
2224 {
2225   int nat,i,j,k;
2226
2227   /* Get rid of the current exclusions and exclude all atom pairs */
2228   nat = excl->nr;
2229   excl->nra = nat*nat;
2230   srenew(excl->a,excl->nra);
2231   k = 0;
2232   for(i=0; i<nat; i++) {
2233     excl->index[i] = k;
2234     for(j=0; j<nat; j++)
2235       excl->a[k++] = j;
2236   }
2237   excl->index[nat] = k;
2238 }
2239
2240 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
2241                            int couple_lam0,int couple_lam1)
2242 {
2243   int i;
2244
2245   for(i=0; i<atoms->nr; i++) {
2246     if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW) {
2247       atoms->atom[i].q     = 0.0;
2248     }
2249     if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ) {
2250       atoms->atom[i].type  = atomtype_decouple;
2251     }
2252     if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW) {
2253       atoms->atom[i].qB    = 0.0;
2254     }
2255     if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ) {
2256       atoms->atom[i].typeB = atomtype_decouple;
2257     }
2258   }
2259 }
2260
2261 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
2262                             int couple_lam0,int couple_lam1,
2263                             bool bCoupleIntra,int nb_funct,t_params *nbp)
2264 {
2265   convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
2266   if (!bCoupleIntra) {
2267     generate_LJCpairsNB(mol,nb_funct,nbp);
2268     set_excl_all(&mol->excls);
2269   }
2270   decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
2271 }