a5eff4763fe379618778bdf22c8b254a50b2f85a
[namd.git] / src / Parameters.C
1 /**
2 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
3 ***  The Board of Trustees of the University of Illinois.
4 ***  All rights reserved.
5 **/
6
7 /*
8    The class Parameters is used to hold all of the parameters read
9    in from the parameter files.  The class provides a routine to read in
10    parameter files (as many parameter files as desired can be read in) and
11    a series of routines that allow the parameters that have been read in
12    to be queried.
13 */
14
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <string.h>
18 #include <vector>
19 #ifndef WIN32
20 #include <strings.h>
21 #endif
22 #include "InfoStream.h"
23 #include <charm++.h>
24 #include "Parameters.h"
25 #include "Communicate.h"
26 #include "ConfigList.h"
27 //****** BEGIN CHARMM/XPLOR type changes
28 #include "SimParameters.h"
29 //****** END CHARMM/XPLOR type changes
30
31 #define MIN_DEBUG_LEVEL 3
32 //#define DEBUGM
33 #include "Debug.h"
34
35 #define INDEX(ncols,i,j)  ((i)*ncols + (j))
36
37 #define ENABLETABLES
38
39 static char** table_types;
40
41 //  struct bond_params is used to form a binary tree of bond parameters.
42 //  The two atom names are used to determine the order of the nodes in the
43 //  tree.  atom1name should ALWAYS be lexically before atom2name
44
45 struct bond_params
46 {
47   char atom1name[11];
48   char atom2name[11];
49   Real forceconstant;
50   Real distance;
51   Index index;
52   struct bond_params *left;
53   struct bond_params *right;
54 };
55
56 //  struct angle_params is used to form a binary tree of bond parameters.
57 //  The three atom names are used to determine the order of the nodes in
58 //  the tree.  atom1name should ALWAYS be lexically before atom3name
59
60 struct angle_params
61 {
62   char atom1name[11];
63   char atom2name[11];
64   char atom3name[11];
65   Real forceconstant;
66   int normal;
67   Real angle;
68   Real k_ub;
69   Real r_ub;
70   Index index;
71   struct angle_params *left;
72   struct angle_params *right;
73 };
74
75 //  struct dihedral_params is used to form a linked list of the dihedral
76 //  parameters.  The linked list is arranged in such a way that any
77 //  bonds with wildcards are at the end of the list so that a linear
78 //  search can be done but we will still find exact matches before
79 //  wildcard matches
80
81 struct dihedral_params
82 {
83   char atom1name[11];
84   char atom2name[11];
85   char atom3name[11];
86   char atom4name[11];
87   char atom1wild;
88   char atom2wild;
89   char atom3wild;
90   char atom4wild;
91   int multiplicity;
92   FourBodyConsts values[MAX_MULTIPLICITY];
93   Index index;
94   struct dihedral_params *next;
95   dihedral_params() { memset(this, 0, sizeof(dihedral_params)); }
96 };
97
98 //  struct improper_params is used to form a linked list of the improper
99 //  parameters.  The linked list is arranged in such a way that any
100 //  bonds with wildcards are at the end of the list so that a linear
101 //  search can be done but we will still find exact matches before
102 //  wildcard matches
103
104 struct improper_params
105 {
106   char atom1name[11];
107   char atom2name[11];
108   char atom3name[11];
109   char atom4name[11];
110   int multiplicity;
111   FourBodyConsts values[MAX_MULTIPLICITY];
112   Index index;
113   struct improper_params *next;
114 };
115
116 struct crossterm_params
117 {
118   crossterm_params(int dim) : dimension(dim) {
119     values = new double[dimension*dimension];
120   }
121   ~crossterm_params() {
122     delete [] values;
123   }
124   char atom1name[11];
125   char atom2name[11];
126   char atom3name[11];
127   char atom4name[11];
128   char atom5name[11];
129   char atom6name[11];
130   char atom7name[11];
131   char atom8name[11];
132   int dimension;  // usually 24
133   double *values;  // dimension * dimension data
134   Index index;
135   struct crossterm_params *next;
136 };
137
138 //  struct vdw_params is used to form a binary serach tree of the
139 //  vdw paramters for a single atom.
140
141 struct vdw_params
142 {
143   char atomname[11];
144   Real sigma;
145   Real epsilon;
146   Real sigma14;
147   Real epsilon14;
148   Index index;
149   struct vdw_params *left;
150   struct vdw_params *right;
151 };
152
153 //  struct vdw_pair_params is used to form a linked list of the
154 //  vdw parameters for a pair of atoms
155
156 struct vdw_pair_params
157 {
158   char atom1name[11];
159   char atom2name[11];
160   Real A;
161   Real B;
162   Real A14;
163   Real B14;
164   struct vdw_pair_params *next;
165 };
166
167 struct table_pair_params
168 {
169   char atom1name[11];
170   char atom2name[11];
171   int K;
172   struct table_pair_params *next;
173 };
174
175 struct nbthole_pair_params
176 {
177   char atom1name[11];
178   char atom2name[11];
179   Real alphai;
180   Real alphaj;
181   Real tholeij;
182   Index index;
183   struct nbthole_pair_params *next;
184 };
185
186 Parameters::Parameters() {
187   initialize();
188 }
189
190 void Parameters::initialize() {
191
192   paramType = -1;
193
194   /*  Set all the pointers to NULL        */
195   atomTypeNames=NULL;
196   bondp=NULL;
197   anglep=NULL;
198   improperp=NULL;
199   dihedralp=NULL;
200   crosstermp=NULL;
201   vdwp=NULL;
202   vdw_pairp=NULL;
203   nbthole_pairp=NULL;
204   table_pairp=NULL;
205   bond_array=NULL;
206   angle_array=NULL;
207   dihedral_array=NULL;
208   improper_array=NULL;
209   crossterm_array=NULL;
210   // JLai
211   gromacsPair_array=NULL;
212   // End of JLai
213   vdw_array=NULL;
214   vdw_pair_tree=NULL;
215   nbthole_pair_tree=NULL;
216   tab_pair_tree=NULL;
217   maxDihedralMults=NULL;
218   maxImproperMults=NULL;
219   table_ener = NULL;
220
221   /*  Set all the counts to 0          */
222   NumBondParams=0;
223   NumAngleParams=0;
224   NumDihedralParams=0;
225   NumImproperParams=0;
226   NumCrosstermParams=0;
227   // JLai
228   NumGromacsPairParams=0;
229   // End of JLai
230   NumVdwParams=0;
231   NumVdwPairParams=0;
232   NumNbtholePairParams=0;
233   NumTablePairParams=0;
234   NumCosAngles=0;
235   numenerentries=0;
236 }
237
238 /************************************************************************/
239 /*                  */
240 /*      FUNCTION Parameters        */
241 /*                  */
242 /*  This is the constructor for the class.  It simlpy sets the      */
243 /*  pointers to the list and trees to NULL and the count of all the     */
244 /*  parameters to 0.              */
245 /*  The type (format) of the input parameters (Xplor,Charmm) is set here. */
246 /*                  */
247 /************************************************************************/
248
249 Parameters::Parameters(SimParameters *simParams, StringList *f)
250 {
251   initialize();
252
253   //// get current parameter format
254   if (simParams->paraTypeXplorOn)
255   {
256     paramType = paraXplor;
257   }
258   else if (simParams->paraTypeCharmmOn)
259   {
260     paramType = paraCharmm;
261   }
262   //****** END CHARMM/XPLOR type changes
263   //Test for cos-based angles
264   if (simParams->cosAngles) {
265     cosAngles = true;
266   } else {
267     cosAngles = false;
268   }
269
270   if (simParams->tabulatedEnergies) {
271           CkPrintf("Working on tables\n");
272           read_ener_table(simParams);
273   }
274
275   //****** BEGIN CHARMM/XPLOR type changes
276   /* Set up AllFilesRead flag to FALSE.  Once all of the files    */
277   /* have been read in, then this will be set to true and the     */
278   /* arrays of parameters will be set up        */
279   AllFilesRead = FALSE;
280
281   if (NULL != f) 
282   {
283     do
284     {
285       //****** BEGIN CHARMM/XPLOR type changes
286       if (paramType == paraXplor)
287       {
288         read_parameter_file(f->data);
289       }
290       else if (paramType == paraCharmm)
291       {
292         read_charmm_parameter_file(f->data);
293       }
294       //****** END CHARMM/XPLOR type changes
295       f = f->next;
296     } while ( f != NULL );
297
298     done_reading_files();
299   }
300
301 }
302 /*      END OF FUNCTION Parameters      */
303
304 /************************************************************************/
305 /*                  */
306 /*      FUNCTION ~Parameters        */
307 /*                  */
308 /*  This is the destructor for this class.  It basically just       */
309 /*  frees all of the memory allocated for the parameters.    */
310 /*                  */
311 /************************************************************************/
312
313 Parameters::~Parameters()
314
315 {
316         if (atomTypeNames)
317           delete [] atomTypeNames;
318
319   if (bondp != NULL)
320     free_bond_tree(bondp);
321
322   if (anglep != NULL)
323     free_angle_tree(anglep);
324
325   if (dihedralp != NULL)
326     free_dihedral_list(dihedralp);
327
328   if (improperp != NULL)
329     free_improper_list(improperp);
330
331   if (crosstermp != NULL)
332     free_crossterm_list(crosstermp);
333
334   if (vdwp != NULL)
335     free_vdw_tree(vdwp);
336
337   if (vdw_pairp != NULL)
338     free_vdw_pair_list();
339
340   if (nbthole_pairp != NULL)
341     free_nbthole_pair_list();
342
343   if (bond_array != NULL)
344     delete [] bond_array;
345
346   if (angle_array != NULL)
347     delete [] angle_array;
348
349   if (dihedral_array != NULL)
350     delete [] dihedral_array;
351
352   if (improper_array != NULL)
353     delete [] improper_array;
354
355   if (crossterm_array != NULL)
356     delete [] crossterm_array;
357
358   // JLai
359   if (gromacsPair_array != NULL)
360     delete [] gromacsPair_array;
361   // End of JLai
362
363   if (vdw_array != NULL)
364     delete [] vdw_array;
365   
366   if (tab_pair_tree != NULL)
367     free_table_pair_tree(tab_pair_tree);
368
369   if (vdw_pair_tree != NULL)
370     free_vdw_pair_tree(vdw_pair_tree);
371
372   if (nbthole_pair_tree != NULL)
373     free_nbthole_pair_tree(nbthole_pair_tree);
374
375   if (maxDihedralMults != NULL)
376     delete [] maxDihedralMults;
377
378   if (maxImproperMults != NULL)
379     delete [] maxImproperMults;
380
381   for( int i = 0; i < error_msgs.size(); ++i ) {
382     delete [] error_msgs[i];
383   }
384   error_msgs.resize(0);
385 }
386 /*      END OF FUNCTION ~Parameters      */
387
388 /************************************************************************/
389 /*                  */
390 /*      FUNCTION read_paramter_file      */
391 /*                  */
392 /*   INPUTS:                */
393 /*  fname - name of the parameter file to read      */
394 /*                  */
395 /*  This function reads in a parameter file and adds the parameters */
396 /*   from this file to the current group of parameters.  The basic      */
397 /*   logic of the routine is to read in a line from the file, looks at  */
398 /*   the first word of the line to determine what kind of parameter we  */
399 /*   have, and then call the appropriate routine to add the parameter   */
400 /*   to the parameters that we have.          */
401 /*                  */
402 /************************************************************************/
403
404 void Parameters::read_parameter_file(char *fname)
405
406 {
407   char buffer[512];  //  Buffer to store each line of the file
408   char first_word[512];  //  First word of the current line
409   FILE *pfile;    //  File descriptor for the parameter file
410
411   /*  Check to make sure that we haven't previously been told     */
412   /*  that all the files were read        */
413   if (AllFilesRead)
414   {
415     NAMD_die("Tried to read another parameter file after being told that all files were read!");
416   }
417
418   /*  Try and open the file          */
419   if ( (pfile = Fopen(fname, "r")) == NULL)
420   {
421     char err_msg[256];
422
423     sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
424     NAMD_die(err_msg);
425   }
426
427   /*  Keep reading in lines until we hit the EOF      */
428   while (NAMD_read_line(pfile, buffer) != -1)
429   {
430     /*  Get the first word of the line      */
431     NAMD_find_first_word(buffer, first_word);
432
433     /*  First, screen out things that we ignore, such as    */
434     /*  blank lines, lines that start with '!', lines that  */
435     /*  start with "REMARK", lines that start with set",    */
436     /*  and most of the HBOND parameters which include      */
437     /*  AEXP, REXP, HAEX, AAEX, but not the HBOND statement */
438     /*  which is parsed.                                    */
439     if ((buffer[0] != '!') && 
440         !NAMD_blank_string(buffer) &&
441         (strncasecmp(first_word, "REMARK", 6) != 0) &&
442         (strcasecmp(first_word, "set")!=0) &&
443         (strncasecmp(first_word, "AEXP", 4) != 0) &&
444         (strncasecmp(first_word, "REXP", 4) != 0) &&
445         (strncasecmp(first_word, "HAEX", 4) != 0) &&
446         (strncasecmp(first_word, "AAEX", 4) != 0) &&
447         (strncasecmp(first_word, "NBOND", 5) != 0) &&
448         (strncasecmp(first_word, "CUTNB", 5) != 0) &&
449         (strncasecmp(first_word, "END", 3) != 0) &&
450         (strncasecmp(first_word, "CTONN", 5) != 0) &&
451         (strncasecmp(first_word, "EPS", 3) != 0) &&
452         (strncasecmp(first_word, "VSWI", 4) != 0) &&
453         (strncasecmp(first_word, "NBXM", 4) != 0) &&
454         (strncasecmp(first_word, "INHI", 4) != 0) )
455     {
456       /*  Now, call the appropriate function based    */
457       /*  on the type of parameter we have    */
458       if (strncasecmp(first_word, "bond", 4)==0)
459       {
460         add_bond_param(buffer);
461         NumBondParams++;
462       }
463       else if (strncasecmp(first_word, "angl", 4)==0)
464       {
465         add_angle_param(buffer);
466         NumAngleParams++;
467       }
468       else if (strncasecmp(first_word, "dihe", 4)==0)
469       {
470         add_dihedral_param(buffer, pfile);
471         NumDihedralParams++;
472       }
473       else if (strncasecmp(first_word, "impr", 4)==0)
474       {
475         add_improper_param(buffer, pfile);
476         NumImproperParams++;
477       }
478       else if (strncasecmp(first_word, "nonb", 4)==0)
479       {
480         add_vdw_param(buffer);
481         NumVdwParams++; 
482       }
483       else if (strncasecmp(first_word, "nbfi", 4)==0)
484       {
485         add_vdw_pair_param(buffer);
486         NumVdwPairParams++; 
487       }
488       else if (strncasecmp(first_word, "nbta", 4)==0)
489       {
490
491         if (table_ener == NULL) {
492           continue;
493         }
494
495         add_table_pair_param(buffer);
496         NumTablePairParams++; 
497       }
498       else if (strncasecmp(first_word, "hbon", 4)==0)
499       {
500         add_hb_pair_param(buffer);
501       }
502       else
503       {
504         /*  This is an unknown paramter.        */
505         /*  This is BAD        */
506         char err_msg[512];
507
508         sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
509            fname, buffer);
510         NAMD_die(err_msg);
511       }
512     }
513   }
514
515   /*  Close the file            */
516   Fclose(pfile);
517
518   return;
519 }
520 /*      END OF FUNCTION read_paramter_file    */
521
522 //****** BEGIN CHARMM/XPLOR type changes
523 /************************************************************************/
524 /*                                                                        */
525 /*                        FUNCTION read_charmm_paramter_file                */
526 /*                                                                        */
527 /*   INPUTS:                                                                */
528 /*        fname - name of the parameter file to read                        */
529 /*                                                                        */
530 /*        This function reads in a CAHRMM parameter file and adds the     */ 
531 /*   parameters from this file to the current group of parameters.      */
532 /*   The basic logic of the routine is to first find out what type of   */
533 /*   parameter we have in the file. Then look at each line in turn      */
534 /*   and call the appropriate routine to add the parameters until we hit*/
535 /*   a new type of parameter or EOF.                                    */
536 /*                                                                        */
537 /************************************************************************/
538
539 void Parameters::read_charmm_parameter_file(char *fname)
540
541 {
542   int  par_type=0;         //  What type of parameter are we currently
543                            //  dealing with? (vide infra)
544   int  skipline;           //  skip this line?
545   int  skipall = 0;        //  skip rest of file;
546   char buffer[512];           //  Buffer to store each line of the file
547   char first_word[512];           //  First word of the current line
548   FILE *pfile;                   //  File descriptor for the parameter file
549
550   /*  Check to make sure that we haven't previously been told     */
551   /*  that all the files were read                                */
552   if (AllFilesRead)
553   {
554     NAMD_die("Tried to read another parameter file after being told that all files were read!");
555   }
556
557   /*  Try and open the file                                        */
558   if ( (pfile = fopen(fname, "r")) == NULL)
559   {
560     char err_msg[256];
561
562     sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
563     NAMD_die(err_msg);
564   }
565
566   /*  Keep reading in lines until we hit the EOF                        */
567   while (NAMD_read_line(pfile, buffer) != -1)
568   {
569     /*  Get the first word of the line                        */
570     NAMD_find_first_word(buffer, first_word);
571     skipline=0;
572
573     /*  First, screen out things that we ignore.                   */   
574     /*  blank lines, lines that start with '!' or '*', lines that  */
575     /*  start with "END".                                          */
576     if (!NAMD_blank_string(buffer) &&
577         (strncmp(first_word, "!", 1) != 0) &&
578          (strncmp(first_word, "*", 1) != 0) &&
579         (strncasecmp(first_word, "END", 3) != 0))
580     {
581       if ( skipall ) {
582         iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
583         break;
584       }
585       /*  Now, determine the apropriate parameter type.   */
586       if (strncasecmp(first_word, "bond", 4)==0)
587       {
588         par_type=1; skipline=1;
589       }
590       else if (strncasecmp(first_word, "angl", 4)==0)
591       {
592         par_type=2; skipline=1;
593       }
594       else if (strncasecmp(first_word, "thet", 4)==0)
595       {
596         par_type=2; skipline=1;
597       }
598       else if (strncasecmp(first_word, "dihe", 4)==0)
599       {
600         par_type=3; skipline=1;
601       }
602       else if (strncasecmp(first_word, "phi", 3)==0)
603       {
604         par_type=3; skipline=1;
605       }
606       else if (strncasecmp(first_word, "impr", 4)==0)
607       {
608         par_type=4; skipline=1;
609       }
610       else if (strncasecmp(first_word, "imph", 4)==0)
611       {
612         par_type=4; skipline=1;
613       }
614       else if (strncasecmp(first_word, "nbon", 4)==0)
615       {
616         par_type=5; skipline=1;
617       }
618       else if (strncasecmp(first_word, "nonb", 4)==0)
619       {
620         par_type=5; skipline=1;
621       }
622       else if (strncasecmp(first_word, "nbfi", 4)==0)
623       {
624         par_type=6; skipline=1;
625       }
626       else if (strncasecmp(first_word, "hbon", 4)==0)
627       {
628         par_type=7; skipline=1;
629       }
630       else if (strncasecmp(first_word, "cmap", 4)==0)
631       {
632         par_type=8; skipline=1;
633       }
634       else if (strncasecmp(first_word, "nbta", 4)==0)
635       {
636         par_type=9; skipline=1;
637       }
638       else if (strncasecmp(first_word, "thol", 4)==0)
639       {
640         par_type=10; skipline=1;
641       }
642       else if (strncasecmp(first_word, "atom", 4)==0)
643       {
644         par_type=11; skipline=1;
645       }
646       else if (strncasecmp(first_word, "ioformat", 8)==0)
647       {
648         skipline=1;
649       }
650       else if (strncasecmp(first_word, "read", 4)==0)
651       {
652         skip_stream_read(buffer, pfile);  skipline=1;
653       }
654       else if (strncasecmp(first_word, "return", 4)==0)
655       {
656         skipall=1;  skipline=1;
657       }
658       else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
659                (strncasecmp(first_word, "grou", 4) == 0) ||
660                (strncasecmp(first_word, "cdie", 4) == 0) ||
661                (strncasecmp(first_word, "shif", 4) == 0) ||
662                (strncasecmp(first_word, "vgro", 4) == 0) ||
663                (strncasecmp(first_word, "vdis", 4) == 0) ||
664                (strncasecmp(first_word, "vswi", 4) == 0) ||
665                (strncasecmp(first_word, "cutn", 4) == 0) ||
666                (strncasecmp(first_word, "ctof", 4) == 0) ||
667                (strncasecmp(first_word, "cton", 4) == 0) ||
668                (strncasecmp(first_word, "eps", 3) == 0) ||
669                (strncasecmp(first_word, "e14f", 4) == 0) ||
670                (strncasecmp(first_word, "wmin", 4) == 0) ||
671                (strncasecmp(first_word, "aexp", 4) == 0) ||
672                (strncasecmp(first_word, "rexp", 4) == 0) ||
673                (strncasecmp(first_word, "haex", 4) == 0) ||
674                (strncasecmp(first_word, "aaex", 4) == 0) ||
675                (strncasecmp(first_word, "noac", 4) == 0) ||
676                (strncasecmp(first_word, "hbno", 4) == 0) ||
677                (strncasecmp(first_word, "cuth", 4) == 0) ||
678                (strncasecmp(first_word, "ctof", 4) == 0) ||
679                (strncasecmp(first_word, "cton", 4) == 0) ||
680                (strncasecmp(first_word, "cuth", 4) == 0) ||
681                (strncasecmp(first_word, "ctof", 4) == 0) ||
682                (strncasecmp(first_word, "cton", 4) == 0) ) 
683       {
684         if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
685         {
686           char err_msg[512];
687
688           sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
689           NAMD_die(err_msg);
690         }
691         else 
692         {
693           skipline = 1;
694         }
695       }        
696       else if (par_type == 0)
697       {
698         /*  This is an unknown paramter.        */
699         /*  This is BAD                                */
700         char err_msg[512];
701
702         sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
703         NAMD_die(err_msg);
704       }
705     }
706     else
707     {
708       skipline=1;
709     }
710
711     if ( (par_type != 0) && (!skipline) )
712     {
713       /*  Now, call the appropriate function based    */
714       /*  on the type of parameter we have                */
715       /*  I know, this should really be a switch ...  */
716       if (par_type == 1)
717       {
718         add_bond_param(buffer);
719         NumBondParams++;
720       }
721       else if (par_type == 2)
722       {
723         add_angle_param(buffer);
724         NumAngleParams++;
725       }
726       else if (par_type == 3)
727       {
728         add_dihedral_param(buffer, pfile);
729         NumDihedralParams++;
730       }
731       else if (par_type == 4)
732       {
733         add_improper_param(buffer, pfile);
734         NumImproperParams++;
735       }
736       else if (par_type == 5)
737       {
738         add_vdw_param(buffer);
739         NumVdwParams++;
740       }
741       else if (par_type == 6)
742       {
743         add_vdw_pair_param(buffer);
744         NumVdwPairParams++; 
745       }
746       else if (par_type == 7)
747       {
748         add_hb_pair_param(buffer);                  
749       }
750       else if (par_type == 8)
751       {
752         add_crossterm_param(buffer, pfile);                  
753         NumCrosstermParams++;
754       }
755       else if (par_type == 9)
756       {
757
758         if (table_ener == NULL) {
759           continue;
760         }
761
762         add_table_pair_param(buffer);                  
763         NumTablePairParams++;
764       }
765       else if (par_type == 10)
766       {
767         add_nbthole_pair_param(buffer);
768         NumNbtholePairParams++;
769       }
770       else if (par_type == 11)
771       {
772         if ( strncasecmp(first_word, "mass", 4) != 0 ) {
773           char err_msg[512];
774           sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
775           NAMD_die(err_msg);
776         }
777       }
778       else
779       {
780         /*  This really should not occour!      */
781         /*  This is an internal error.          */
782         /*  This is VERY BAD                        */
783         char err_msg[512];
784
785         sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
786         NAMD_die(err_msg);
787       }
788     }
789   }
790
791   /*  Close the file                                                */
792   fclose(pfile);
793
794   return;
795 }
796 /*                        END OF FUNCTION read_charmm_paramter_file                */
797 //****** END CHARMM/XPLOR type changes
798
799
800 void Parameters::skip_stream_read(char *buf, FILE *fd) {
801
802   char buffer[513];
803   char first_word[513];
804   char s1[128];
805   char s2[128];
806   int rval = sscanf(buf, "%s %s", s1, s2);
807   if (rval != 2) {
808         char err_msg[512];
809         sprintf(err_msg, "BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
810         NAMD_die(err_msg);
811   }
812   if ( ! strncasecmp(s2,"PARA",4) ) return;  // read parameters
813   
814   iout << iINFO << "SKIPPING " << s2 << " SECTION IN STREAM FILE\n" << endi;
815
816   while (NAMD_read_line(fd, buffer) != -1)
817   {
818     // read until we find "END"
819     NAMD_find_first_word(buffer, first_word);
820     if (!NAMD_blank_string(buffer) &&
821         (strncmp(first_word, "!", 1) != 0) &&
822          (strncmp(first_word, "*", 1) != 0) &&
823          (strncasecmp(first_word, "END", 3) == 0) ) return;
824   }
825
826 }
827
828
829 /************************************************************************/
830 /*                  */
831 /*      FUNCTION add_bond_param        */
832 /*                  */
833 /*   INPUTS:                */
834 /*  buf - Line from parameter file containing bond parameters  */
835 /*                  */
836 /*  This function adds a new bond paramter to the binary tree of    */
837 /*   angle paramters that we have.  If a duplicate is found, a warning  */
838 /*   message is printed and the new parameters are used.    */
839 /*                  */
840 /************************************************************************/
841
842 void Parameters::add_bond_param(char *buf)
843
844 {
845   char atom1name[11];    //  Atom type for atom 1
846   char atom2name[11];    //  Atom type for atom 2
847   Real forceconstant;    //  Force constant for bond
848   Real distance;      //  Rest distance for bond
849   int read_count;      //  Count from sscanf
850   struct bond_params *new_node;  //  New node in tree
851
852   //****** BEGIN CHARMM/XPLOR type changes
853   /*  Use sscanf to parse up the input line      */
854   if (paramType == paraXplor)
855   {
856     /* read XPLOR format */
857     read_count=sscanf(buf, "%*s %s %s %f %f\n", atom1name, atom2name, 
858        &forceconstant, &distance);
859   }
860   else if (paramType == paraCharmm)
861   {
862     /* read CHARMM format */
863     read_count=sscanf(buf, "%s %s %f %f\n", atom1name, atom2name, 
864        &forceconstant, &distance);
865   }
866   //****** END CHARMM/XPLOR type changes
867
868   /*  Check to make sure we found everything we expeceted    */
869   if (read_count != 4)
870   {
871     char err_msg[512];
872
873     if (paramType == paraXplor)
874     {
875       sprintf(err_msg, "BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
876     }
877     else if (paramType == paraCharmm)
878     {
879       sprintf(err_msg, "BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
880     }
881     NAMD_die(err_msg);
882   }
883
884   /*  Allocate a new node            */
885   new_node = new bond_params;
886
887   if (new_node == NULL)
888   {
889     NAMD_die("memory allocation failed in Parameters::add_bond_param\n");
890   }
891
892   /*  Order the atoms so that the atom that comes alphabetically  */
893   /*  first is atom 1.  Since the bond is symmetric, it doesn't   */
894   /*  matter physically which atom is first.  And this allows the */
895   /*  search of the binary tree to be done in a logical manner    */
896   if (strcasecmp(atom1name, atom2name) < 0)
897   {
898     strcpy(new_node->atom1name, atom1name);
899     strcpy(new_node->atom2name, atom2name);
900   }
901   else
902   {
903     strcpy(new_node->atom2name, atom1name);
904     strcpy(new_node->atom1name, atom2name);
905   }
906
907   /*  Assign force constant and distance        */
908   new_node->forceconstant = forceconstant;
909   new_node->distance = distance;
910
911   /*  Set pointers to null          */
912   new_node->left = NULL;
913   new_node->right = NULL;
914
915   /*  Make call to recursive call to actually add the node to the */
916   /*  tree              */
917   bondp=add_to_bond_tree(new_node, bondp);
918
919   return;
920 }
921 /*      END OF FUNCTION add_bond_param      */
922
923 /************************************************************************/
924 /*                  */
925 /*      FUNCTION add_to_bond_tree      */
926 /*                  */
927 /*   INPUTS:                */
928 /*  new_node - Node to add to the tree        */
929 /*  tree - tree to add the node to          */
930 /*                  */
931 /*   OUTPUTS:                */
932 /*  ths function returns a pointer to the new tree with the node    */
933 /*   added to it.  Most of the time it will be the same pointer as was  */
934 /*   passed in, but not if the current tree is empty.      */
935 /*                  */
936 /*  this is a receursive function that adds a node to the binary    */
937 /*   tree used to store bond parameters.        */
938 /*                  */
939 /************************************************************************/
940
941 struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
942              struct bond_params *tree)
943
944 {
945   int compare_code;  //  Results from strcasecmp
946
947   /*  If the tree is currently empty, then the new tree consists  */
948   /*  only of the new node          */
949   if (tree == NULL)
950     return(new_node);
951
952   /*  Compare the atom1 name from the new node and the head of    */
953   /*  the tree              */
954   compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
955
956   /*  Check to see if they are the same        */
957   if (compare_code == 0)
958   {
959     /*  The atom 1 names are the same, compare atom 2  */
960     compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
961
962     /*  If atom 1 AND atom 2 are the same, we have a duplicate */
963     if (compare_code == 0)
964     {
965       /*  We have a duplicate.  So print out a warning*/
966       /*  message.  Then assign the new values to the */
967       /*  tree and free the new_node      */
968       //****** BEGIN CHARMM/XPLOR type changes
969       /* we do not care about identical replacement */
970       if ((tree->forceconstant != new_node->forceconstant) || 
971           (tree->distance != new_node->distance))
972       {
973         iout << "\n" << iWARN << "DUPLICATE BOND ENTRY FOR "
974           << new_node->atom1name << "-"
975           << new_node->atom2name
976           << "\nPREVIOUS VALUES  k=" << tree->forceconstant
977           << "  x0=" << tree->distance
978           << "\n   USING VALUES  k=" << new_node->forceconstant
979           << "  x0=" << new_node->distance
980           << "\n" << endi;
981
982         tree->forceconstant=new_node->forceconstant;
983         tree->distance=new_node->distance;
984       }
985       //****** END CHARMM/XPLOR type changes
986
987       delete new_node;
988
989       return(tree);
990     }
991   }
992
993   /*  We don't have a duplicate, so if the new value is less      */
994   /*  than the head of the tree, add it to the left child,   */
995   /*  otherwise add it to the right child        */
996   if (compare_code < 0)
997   {
998     tree->left = add_to_bond_tree(new_node, tree->left);
999   }
1000   else
1001   {
1002     tree->right = add_to_bond_tree(new_node, tree->right);
1003   }
1004
1005   return(tree);
1006 }
1007 /*    END OF FUNCTION add_to_bond_tree      */
1008
1009 /************************************************************************/
1010 /*                  */
1011 /*      FUNCTION add_angle_param      */
1012 /*                  */
1013 /*   INPUTS:                */
1014 /*  buf - line from paramter file with angle parameters    */
1015 /*                  */
1016 /*  this function adds an angle parameter.  It parses up the input  */
1017 /*   line and then adds it to the binary tree used to store the angle   */
1018 /*   parameters.              */
1019 /*                  */
1020 /************************************************************************/
1021
1022 void Parameters::add_angle_param(char *buf)
1023
1024 {
1025   char atom1name[11];    // Type for atom 1
1026   char atom2name[11];    // Type for atom 2
1027   char atom3name[11];    // Type for atom 3
1028   char norm[4]="xxx";
1029   Real forceconstant;    // Force constant
1030   Real angle;      // Theta 0
1031   Real k_ub;      // Urey-Bradley force constant
1032   Real r_ub;      // Urey-Bradley distance
1033   int read_count;      // count from sscanf
1034   struct angle_params *new_node;  // new node in tree
1035
1036   //****** BEGIN CHARMM/XPLOR type changes
1037   /*  parse up the input line with sscanf        */
1038   if (paramType == paraXplor)
1039   {
1040     /* read XPLOR format */
1041     read_count=sscanf(buf, "%*s %s %s %s %f %f UB %f %f\n", 
1042        atom1name, atom2name, atom3name, &forceconstant, &angle,
1043        &k_ub, &r_ub);
1044   }
1045   else if ((paramType == paraCharmm) && cosAngles) {
1046     read_count=sscanf(buf, "%s %s %s %f %f %3s %f %f\n", 
1047        atom1name, atom2name, atom3name, &forceconstant, &angle, norm,
1048        &k_ub, &r_ub);
1049 //    printf("%s\n", buf);
1050 //    printf("Data: %s %s %s %f %f %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle, norm, k_ub, r_ub);
1051   }  
1052   else if (paramType == paraCharmm)
1053   {
1054     /* read CHARMM format */
1055     read_count=sscanf(buf, "%s %s %s %f %f %f %f\n", 
1056        atom1name, atom2name, atom3name, &forceconstant, &angle,
1057        &k_ub, &r_ub);
1058 //    printf("%s\n", buf);
1059 //    printf("Data: %s %s %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle);
1060   }
1061   //****** END CHARMM/XPLOR type changes
1062
1063   /*  Check to make sure we got what we expected      */
1064   if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
1065   {
1066     char err_msg[512];
1067
1068     if (paramType == paraXplor)
1069     {
1070       sprintf(err_msg, "BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1071     }
1072     else if (paramType == paraCharmm)
1073     {
1074       sprintf(err_msg, "BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1075     }
1076     NAMD_die(err_msg);
1077   }
1078
1079   /*  Allocate the new node          */
1080   new_node = new angle_params;
1081
1082   if (new_node == NULL)
1083   {
1084     NAMD_die("memory allocation failed in Parameters::add_angle_param");
1085   }
1086
1087   /*  As with the bond, we want the atom type is comes first  */
1088   /*  alphbetically first between atom 1 and atom 3 to be in      */
1089   /*  atom 1 so that we can search the tree reliably.    */
1090   if (strcasecmp(atom1name, atom3name) < 0)
1091   {
1092     strcpy(new_node->atom1name, atom1name);
1093     strcpy(new_node->atom2name, atom2name);
1094     strcpy(new_node->atom3name, atom3name);
1095   }
1096   else
1097   {
1098     strcpy(new_node->atom3name, atom1name);
1099     strcpy(new_node->atom2name, atom2name);
1100     strcpy(new_node->atom1name, atom3name);
1101   }
1102
1103   /*  Assign the constants and pointer values      */
1104   new_node->forceconstant = forceconstant;
1105   new_node->angle = angle;
1106
1107   if (cosAngles) {
1108     if (strcasecmp("cos",norm)==0) {
1109 //      iout << "Info: Using cos mode for angle " << buf << endl;
1110       NumCosAngles++;
1111       new_node->normal = 0;
1112     } else {
1113 //      iout << "Info: Using x^2 mode for angle " << buf << endl;
1114       new_node->normal = 1;
1115     }
1116   } else {
1117     new_node->normal = 1;
1118   }
1119
1120   if (read_count == 7)
1121   {
1122     //  Urey-Bradley constants
1123     if (new_node->normal == 0) {
1124       NAMD_die("ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
1125     }
1126     new_node->k_ub = k_ub;
1127     new_node->r_ub = r_ub;
1128   }
1129   else
1130   {
1131     new_node->k_ub = 0.0;
1132     new_node->r_ub = 0.0;
1133   }
1134
1135   new_node->left = NULL;
1136   new_node->right = NULL;
1137
1138   /*  Insert it into the tree          */
1139   anglep = add_to_angle_tree(new_node, anglep);
1140
1141   return;
1142 }
1143 /*      END OF FUNCTION add_angle_param      */
1144
1145 /************************************************************************/
1146 /*                  */
1147 /*      FUNCTION add_to_angle_tree      */
1148 /*                  */
1149 /*   INPUTS:                */
1150 /*  new_node - new node to add to the angle tree      */
1151 /*  tree - tree to add the node to          */
1152 /*                  */
1153 /*   OUTPUTS:                */
1154 /*  the function returns a pointer to the new tree with the node    */
1155 /*   added.  Most of the time, this will be the same as the value passed*/
1156 /*   in, but not in the case where the tree is empty.      */
1157 /*                  */
1158 /*  this is a recursive function that adds an angle parameter  */
1159 /*   to the binary tree storing the angle parameters.  If a duplicate   */
1160 /*   is found, a warning message is printed, the current values in the  */
1161 /*   tree are replaced with the new values, and the new node is free'd  */
1162 /*                  */
1163 /************************************************************************/
1164
1165 struct angle_params *Parameters::add_to_angle_tree(struct angle_params *new_node,
1166              struct angle_params *tree)
1167
1168 {
1169   int compare_code;  //  Return code from strcasecmp
1170
1171   /*  If the tree is empty, then the new_node is the tree    */
1172   if (tree == NULL)
1173     return(new_node);
1174
1175   /*  Compare atom 1 from the new node and the head of the tree   */
1176   compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
1177
1178   if (compare_code == 0)
1179   {
1180     /*  Atom 1 is the same, compare atom 2      */
1181     compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
1182
1183     if (compare_code == 0)
1184     {
1185       /*  Atoms 1 & 2 are the same, compare atom 3  */
1186       compare_code = strcasecmp(new_node->atom3name, 
1187             tree->atom3name);
1188
1189       if (compare_code == 0)
1190       {
1191         /*  All three atoms were the same, this */
1192         /*  is a duplicate.  Print a warning    */
1193         /*  message, replace the current values,*/
1194         /*  and free the new node    */
1195         //****** BEGIN CHARMM/XPLOR type changes
1196         /* we do not care about identical replacement */
1197         if ((tree->forceconstant != new_node->forceconstant) ||
1198             (tree->angle != new_node->angle) ||
1199             (tree->k_ub != new_node->k_ub) ||
1200             (tree->r_ub != new_node->r_ub) || (tree->normal != new_node->normal))
1201         {
1202           iout << "\n" << iWARN << "DUPLICATE ANGLE ENTRY FOR "
1203             << new_node->atom1name << "-"
1204             << new_node->atom2name << "-"
1205             << new_node->atom3name
1206             << "\nPREVIOUS VALUES  k="
1207             << tree->forceconstant << "  theta0="
1208             << tree->angle << " k_ub="
1209             << tree->k_ub << " r_ub="
1210             << tree->r_ub
1211             << "\n   USING VALUES  k="
1212             << new_node->forceconstant << "  theta0="
1213             << new_node->angle << " k_ub="
1214             << new_node->k_ub << " r_ub=" << new_node->r_ub 
1215             << "\n" << endi;
1216
1217           tree->forceconstant=new_node->forceconstant;
1218           tree->angle=new_node->angle;
1219           tree->k_ub=new_node->k_ub;
1220           tree->r_ub=new_node->r_ub;
1221           tree->normal=new_node->normal;
1222         }
1223         //****** END CHARMM/XPLOR type changes
1224
1225         delete new_node;
1226
1227         return(tree);
1228       }
1229     }
1230   }
1231
1232   /*  Didn't find a duplicate, so if the new_node is smaller  */
1233   /*  than the current head, add it to the left child.  Otherwise */
1234   /*  add it to the right child.          */
1235   if (compare_code < 0)
1236   {
1237     tree->left = add_to_angle_tree(new_node, tree->left);
1238   }
1239   else
1240   {
1241     tree->right = add_to_angle_tree(new_node, tree->right);
1242   }
1243
1244   return(tree);
1245 }
1246 /*      END OF FUNCTION add_to_angle_tree    */
1247
1248 /************************************************************************/
1249 /*                  */
1250 /*      FUNCTION add_dihedral_param      */
1251 /*                  */
1252 /*   INPUTS:                */
1253 /*  buf - line from paramter file with dihedral parameters    */
1254 /*                  */
1255 /*  this function adds an dihedral parameter.  It parses up the     */
1256 /*   input line and then adds it to the binary tree used to store the   */
1257 /*   dihedral parameters.            */
1258 /*                  */
1259 /************************************************************************/
1260
1261 void Parameters::add_dihedral_param(char *buf, FILE *fd)
1262
1263 {
1264   char atom1name[11];       //  Type of atom 1
1265   char atom2name[11];       //  Type of atom 2
1266   char atom3name[11];       //  Type of atom 3
1267   char atom4name[11];       //  Type of atom 4
1268   Real forceconstant;       //  Force constant
1269   int periodicity;       //  Periodicity
1270   Real phase_shift;       //  Phase shift
1271   int read_count;         //  Count from sscanf
1272   struct dihedral_params *new_node;  //  New node
1273   int multiplicity;       //  Multiplicity for bonds
1274   int i;           //  Loop counter
1275   char buffer[513];       //  Buffer for new line
1276   int ret_code;         //  Return code
1277
1278   //****** BEGIN CHARMM/XPLOR type changes
1279   /*  Parse up the input line using sscanf      */
1280   if (paramType == paraXplor)
1281   {
1282     /* read XPLOR format */
1283     read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n", 
1284        atom1name, atom2name, atom3name, atom4name, &multiplicity,
1285        &forceconstant, &periodicity, &phase_shift);
1286   }
1287   else if (paramType == paraCharmm)
1288   {
1289     /* read CHARMM format */
1290     read_count=sscanf(buf, "%s %s %s %s %f %d %f\n", 
1291        atom1name, atom2name, atom3name, atom4name,
1292        &forceconstant, &periodicity, &phase_shift);
1293     multiplicity=1; 
1294   }
1295
1296   if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
1297   {
1298     char err_msg[512];
1299
1300     sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1301     NAMD_die(err_msg);
1302   }
1303   else if ( (read_count != 7) && (paramType == paraCharmm) )
1304   {
1305     char err_msg[512];
1306
1307     sprintf(err_msg, "BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1308     NAMD_die(err_msg);
1309   }
1310
1311   if ( (read_count == 4) && (paramType == paraXplor) )
1312   //****** END CHARMM/XPLOR type changes
1313   {
1314     read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n", 
1315           &forceconstant, &periodicity, &phase_shift);
1316
1317     /*  Check to make sure we got what we expected    */
1318     if (read_count != 3)
1319     {
1320       char err_msg[512];
1321
1322       sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1323       NAMD_die(err_msg);
1324     }
1325
1326     multiplicity = 1;
1327   }
1328
1329   if (multiplicity > MAX_MULTIPLICITY)
1330   {
1331     char err_msg[181];
1332
1333     sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
1334        multiplicity, MAX_MULTIPLICITY);
1335     NAMD_die(err_msg);
1336   }
1337
1338   /*  Allocate new node            */
1339   new_node = new dihedral_params;
1340
1341   if (new_node == NULL)
1342   {
1343     NAMD_die("memory allocation failed in Parameters::add_dihedral_param\n");
1344   }
1345
1346   /*  Assign all of the values for this node.  Notice that since  */
1347   /*  the dihedrals and impropers are implemented with a linked   */
1348   /*  list rather than a binary tree, we don't really care about  */
1349   /*  the order of the atoms any more        */
1350   strcpy(new_node->atom1name, atom1name);
1351   strcpy(new_node->atom2name, atom2name);
1352   strcpy(new_node->atom3name, atom3name);
1353   strcpy(new_node->atom4name, atom4name);
1354   new_node->atom1wild = ! strcasecmp(atom1name, "X");
1355   new_node->atom2wild = ! strcasecmp(atom2name, "X");
1356   new_node->atom3wild = ! strcasecmp(atom3name, "X");
1357   new_node->atom4wild = ! strcasecmp(atom4name, "X");
1358   new_node->multiplicity = multiplicity;
1359   if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1360   new_node->values[0].k = forceconstant;
1361   new_node->values[0].n = periodicity;
1362   new_node->values[0].delta = phase_shift;
1363
1364   new_node->next = NULL;
1365
1366   //  If the multiplicity is greater than 1, then read in other parameters
1367   if (multiplicity > 1)
1368   {
1369     for (i=1; i<multiplicity; i++)
1370     {
1371       ret_code = NAMD_read_line(fd, buffer);
1372
1373       //  Get rid of comments at the end of a line
1374       if (ret_code == 0)
1375       {
1376         NAMD_remove_comment(buffer);
1377       }
1378
1379       //  Keep reading lines until we get one that isn't blank
1380       while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
1381       {
1382         ret_code = NAMD_read_line(fd, buffer);
1383       }
1384
1385       if (ret_code != 0)
1386       {
1387         NAMD_die("EOF encoutner in middle of multiple dihedral");
1388       }
1389
1390       read_count=sscanf(buffer, "%f %d %f\n", 
1391             &forceconstant, &periodicity, &phase_shift);
1392
1393       if (read_count != 3)
1394       {
1395         char err_msg[512];
1396
1397         sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
1398         NAMD_die(err_msg);
1399       }
1400
1401       if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1402       new_node->values[i].k = forceconstant;
1403       new_node->values[i].n = periodicity;
1404       new_node->values[i].delta = phase_shift;
1405     }
1406   }
1407
1408   //****** BEGIN CHARMM/XPLOR type changes
1409   /*  Add this node to the list          */
1410   if (paramType == paraXplor)
1411   {
1412     add_to_dihedral_list(new_node); // XPLOR
1413   }
1414   else if (paramType == paraCharmm)
1415   {
1416     add_to_charmm_dihedral_list(new_node); // CHARMM
1417   }
1418  //****** END CHARMM/XPLOR type changes
1419
1420   return;
1421 }
1422 /*      END OF FUNCTION add_dihedral_param    */
1423
1424 /************************************************************************/
1425 /*                  */
1426 /*      FUNCTION add_to_dihedral_list      */
1427 /*                  */
1428 /*   INPUTS:                */
1429 /*  new_node - node that is to be added to dihedral_list    */
1430 /*                  */
1431 /*  this function adds a new dihedral parameter to the linked list  */
1432 /*   of dihedral parameters.  First, it checks for duplicates.  If a    */
1433 /*   duplicate is found, a warning message is printed, the old values   */
1434 /*   are replaced with the new values, and the new node is freed.  If   */
1435 /*   Otherwise, the node is added to the list.  This list is arranged   */
1436 /*   so that bods with wildcards are placed at the tail of the list.    */
1437 /*   This will guarantee that if we just do a linear search, we will    */
1438 /*   always find an exact match before a wildcard match.    */
1439 /*                  */
1440 /************************************************************************/
1441
1442 void Parameters::add_to_dihedral_list(
1443         struct dihedral_params *new_node)
1444
1445 {
1446   static struct dihedral_params *ptr;   //  position within list
1447   static struct dihedral_params *tail;  //  Pointer to the end of 
1448                 //  the list so we can add
1449                 //  entries to the end of the
1450                 //  list in constant time
1451   int i;              //  Loop counter
1452
1453   /*  If the list is currently empty, then the new node is the list*/
1454   if (dihedralp == NULL)
1455   {
1456     dihedralp=new_node;
1457     tail=new_node;
1458
1459     return;
1460   }
1461
1462   /*  The list isn't empty, so check for a duplicate    */
1463   ptr=dihedralp;
1464
1465   while (ptr != NULL)
1466   {
1467     if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1468            (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1469            (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1470            (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1471          ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1472            (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1473            (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1474            (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
1475     {
1476       /*  Found a duplicate        */
1477       //****** BEGIN CHARMM/XPLOR type changes
1478       /* we do not care about identical replacement */
1479       int echoWarn=0;  // echo warning messages ?
1480
1481       if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
1482       
1483       if (!echoWarn)
1484       {
1485         for (i=0; i<ptr->multiplicity; i++)
1486         {
1487           if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
1488           if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
1489           if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
1490         }
1491       }
1492
1493       if (echoWarn)
1494       {
1495         iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
1496           << ptr->atom1name << "-"
1497           << ptr->atom2name << "-"
1498           << ptr->atom3name << "-"
1499           << ptr->atom4name
1500           << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
1501         
1502         for (i=0; i<ptr->multiplicity; i++)
1503         {
1504           iout     << "  k=" << ptr->values[i].k
1505                    << "  n=" << ptr->values[i].n
1506                    << "  delta=" << ptr->values[i].delta;
1507         }
1508
1509         iout << "\nUSING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
1510
1511         for (i=0; i<new_node->multiplicity; i++)
1512         {
1513           iout <<     "  k=" << new_node->values[i].k
1514                    << "  n=" << new_node->values[i].n
1515                    << "  delta=" << new_node->values[i].delta;
1516         }
1517
1518         iout << endi;
1519
1520         ptr->multiplicity = new_node->multiplicity;
1521
1522         for (i=0; i<new_node->multiplicity; i++)
1523         {
1524           ptr->values[i].k = new_node->values[i].k;
1525           ptr->values[i].n = new_node->values[i].n;
1526           ptr->values[i].delta = new_node->values[i].delta;
1527         }
1528
1529       }
1530       //****** END CHARMM/XPLOR type changes
1531
1532       delete new_node;
1533
1534       return;
1535     }
1536
1537     ptr=ptr->next;
1538   }
1539
1540   /*  Check to see if we have any wildcards.  Since specific  */
1541   /*  entries are to take precedence, we'll put anything without  */
1542   /*  wildcards at the begining of the list and anything with     */
1543   /*  wildcards at the end of the list.  Then, we can just do a   */
1544   /*  linear search for a bond and be guaranteed to have specific */
1545   /*  entries take precendence over over wildcards          */
1546   if ( new_node->atom1wild ||
1547        new_node->atom2wild ||
1548        new_node->atom3wild ||
1549        new_node->atom4wild )
1550   {
1551     /*  add to the end of the list        */
1552     tail->next=new_node;
1553     tail=new_node;
1554
1555     return;
1556   }
1557   else
1558   {
1559     /*  add to the head of the list        */
1560     new_node->next=dihedralp;
1561     dihedralp=new_node;
1562
1563     return;
1564   }
1565
1566 }
1567 /*    END OF FUNCTION add_to_dihedral_list      */
1568
1569 //****** BEGIN CHARMM/XPLOR type changes
1570 /************************************************************************/
1571 /*                                                                        */
1572 /*                        FUNCTION add_to_charmm_dihedral_list                */
1573 /*                                                                        */
1574 /*   INPUTS:                                                                */
1575 /*        new_node - node that is to be added to dihedral_list                */
1576 /*                                                                        */
1577 /*        this function adds a new dihedral parameter to the linked list  */
1578 /*   of dihedral parameters in CHARMM format.                           */
1579 /*   First, it checks for duplicates.  If a duplicate is found, a       */
1580 /*   warning message is printed. If the periodicity is the same as of   */
1581 /*   a previous dihedral the old values are replaced with the new       */
1582 /*   values, otherwise, the dihedral is added and the multiplicity is   */
1583 /*   increased.                                                         */
1584 /*   Otherwise, the node is added to the list.  This list is arranged   */
1585 /*   so that bonds with wildcards are placed at the tail of the list.   */
1586 /*   This will guarantee that if we just do a linear search, we will    */
1587 /*   always find an exact match before a wildcard match.                */
1588 /*                                                                        */
1589 /************************************************************************/
1590
1591 void Parameters::add_to_charmm_dihedral_list(
1592                                 struct dihedral_params *new_node)
1593
1594 {
1595         static struct dihedral_params *ptr;   //  position within list
1596         static struct dihedral_params *tail;  //  Pointer to the end of 
1597                                               //  the list so we can add
1598                                               //  entries to the end of the
1599                                               //  list in constant time
1600         int i;                                      //  Loop counter
1601         int replace;                          //  replace values?
1602
1603         // keep track of the last dihedral param read to avoid spurious
1604         // error messages.
1605         static struct dihedral_params last_dihedral; 
1606
1607         /*  If the list is currently empty, then the new node is the list*/
1608         if (dihedralp == NULL)
1609         {
1610                 dihedralp=new_node;
1611                 tail=new_node;
1612                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1613
1614                 return;
1615         }
1616
1617         /*  The list isn't empty, so check for a duplicate                */
1618         ptr=dihedralp;
1619
1620         while (ptr != NULL)
1621         {
1622                 int same_as_last = 0;
1623                 if (  ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1624                        (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1625                        (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1626                        (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1627                      ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1628                        (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1629                        (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1630                        (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) )
1631                        )
1632                 {
1633                         /*  Found a duplicate                                */
1634                         
1635                         // check for same_as_last.  Note: don't believe the echoWarn crap; it controls
1636                         // not just whether we print warning messages, but whether we actually change
1637                         // values or not!  
1638
1639                         if ( ( !strcmp(ptr->atom1name, last_dihedral.atom1name) && 
1640                                !strcmp(ptr->atom2name, last_dihedral.atom2name) &&
1641                                !strcmp(ptr->atom3name, last_dihedral.atom3name) &&
1642                                !strcmp(ptr->atom4name, last_dihedral.atom4name)))
1643                           same_as_last = 1;
1644
1645                         //****** BEGIN CHARMM/XPLOR type changes
1646                         /* we do not care about identical replacement */
1647                         int echoWarn=1;  // echo warning messages ?
1648
1649                         // ptr->multiplicity will always be >= new_node->multiplicity
1650                         for (i=0; i<ptr->multiplicity; i++)
1651                         {
1652                           if ((ptr->values[i].k == new_node->values[0].k) && 
1653                               (ptr->values[i].n == new_node->values[0].n) &&
1654                               (ptr->values[i].delta == new_node->values[0].delta)) 
1655                           {
1656                             // found an identical replacement
1657                             echoWarn=0; 
1658                             break;
1659                           }
1660
1661                         }
1662                   
1663                         if (echoWarn)
1664                         {
1665                           if (!same_as_last) {
1666                             iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
1667                                  << ptr->atom1name << "-"
1668                                  << ptr->atom2name << "-"
1669                                  << ptr->atom3name << "-"
1670                                  << ptr->atom4name
1671                                  << "\nPREVIOUS VALUES MULTIPLICITY: " << ptr->multiplicity << "\n";
1672                           }
1673                           replace=0;
1674                           
1675                           for (i=0; i<ptr->multiplicity; i++)
1676                           {
1677                             if (!same_as_last) {
1678                               iout << "  k=" << ptr->values[i].k
1679                                    << "  n=" << ptr->values[i].n
1680                                    << "  delta=" << ptr->values[i].delta << "\n";
1681                             }
1682                             if (ptr->values[i].n == new_node->values[0].n)
1683                             {
1684                               iout << iWARN << "IDENTICAL PERIODICITY! REPLACING OLD VALUES BY: \n";
1685                               ptr->values[i].k = new_node->values[0].k;
1686                               ptr->values[i].delta = new_node->values[0].delta;
1687                               iout << "  k=" << ptr->values[i].k
1688                                    << "  n=" << ptr->values[i].n
1689                                    << "  delta=" << ptr->values[i].delta<< "\n";
1690                               replace=1;
1691                               break;
1692                             }
1693                           }
1694
1695                           if (!replace)
1696                           {
1697                             ptr->multiplicity += 1;
1698
1699                             if (ptr->multiplicity > MAX_MULTIPLICITY)
1700                             {
1701                               char err_msg[181];
1702
1703                               sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
1704                                       ptr->multiplicity, MAX_MULTIPLICITY);
1705                               NAMD_die(err_msg);
1706                             }
1707                             if (!same_as_last) 
1708                               iout << "INCREASING MULTIPLICITY TO: " << ptr->multiplicity << "\n";
1709
1710                             i= ptr->multiplicity - 1; 
1711                             ptr->values[i].k = new_node->values[0].k;
1712                             ptr->values[i].n = new_node->values[0].n;
1713                             ptr->values[i].delta = new_node->values[0].delta;
1714
1715                             if (!same_as_last) 
1716                               iout << "  k=" << ptr->values[i].k
1717                                    << "  n=" << ptr->values[i].n
1718                                    << "  delta=" << ptr->values[i].delta<< "\n";
1719                           }
1720                         
1721                           iout << endi;
1722                         } 
1723                         //****** END CHARMM/XPLOR type changes
1724
1725                         memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1726                         delete new_node;
1727
1728                         return;
1729                 }
1730
1731                 ptr=ptr->next;
1732         }
1733
1734         /*  CHARMM and XPLOR wildcards for dihedrals are luckily the same */
1735         /*  Check to see if we have any wildcards.  Since specific        */
1736         /*  entries are to take precedence, we'll put anything without  */
1737         /*  wildcards at the begining of the list and anything with     */
1738         /*  wildcards at the end of the list.  Then, we can just do a   */
1739         /*  linear search for a bond and be guaranteed to have specific */
1740         /*  entries take precendence over over wildcards                */
1741         if ( new_node->atom1wild ||
1742              new_node->atom2wild ||
1743              new_node->atom3wild ||
1744              new_node->atom4wild )
1745         {
1746                 /*  add to the end of the list                                */
1747                 tail->next=new_node;
1748                 tail=new_node;
1749
1750                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1751                 return;
1752         }
1753         else
1754         {
1755                 /*  add to the head of the list                                */
1756                 new_node->next=dihedralp;
1757                 dihedralp=new_node;
1758
1759                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1760                 return;
1761         }
1762
1763 }
1764 /*                END OF FUNCTION add_to_charmm_dihedral_list                */
1765 //****** END CHARMM/XPLOR type changes
1766
1767 /************************************************************************/
1768 /*                  */
1769 /*      FUNCTION add_improper_param      */
1770 /*                  */
1771 /*   INPUTS:                */
1772 /*  buf - line from paramter file with improper parameters    */
1773 /*                  */
1774 /*  this function adds an improper parameter.  It parses up the     */
1775 /*   input line and then adds it to the binary tree used to store the   */
1776 /*   improper parameters.            */
1777 /*                  */
1778 /************************************************************************/
1779
1780 void Parameters::add_improper_param(char *buf, FILE *fd)
1781
1782 {
1783   char atom1name[11];       //  Atom 1 type
1784   char atom2name[11];       //  Atom 2 type
1785   char atom3name[11];       //  Atom 3 type
1786   char atom4name[11];       //  Atom 4 type
1787   Real forceconstant;       //  Force constant 
1788   int periodicity;       //  Periodicity
1789   Real phase_shift;       //  Phase shift
1790   int read_count;         //  Count from sscanf
1791   struct improper_params *new_node;  //  New node
1792   int multiplicity;       //  Multiplicity for bonds
1793   int i;           //  Loop counter
1794   char buffer[513];       //  Buffer for new line
1795   int ret_code;         //  Return code
1796
1797   //****** BEGIN CHARMM/XPLOR type changes
1798   /*  Parse up the line with sscanf                                */
1799   if (paramType == paraXplor)
1800   {
1801     /* read XPLOR format */
1802     read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n", 
1803        atom1name, atom2name, atom3name, atom4name, &multiplicity, 
1804        &forceconstant, &periodicity, &phase_shift);
1805   }
1806   else if (paramType == paraCharmm)
1807   {
1808     /* read CHARMM format */
1809     read_count=sscanf(buf, "%s %s %s %s %f %d %f\n", 
1810        atom1name, atom2name, atom3name, atom4name,  
1811        &forceconstant, &periodicity, &phase_shift); 
1812     multiplicity=1;      
1813   }
1814
1815   if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
1816   {
1817     char err_msg[512];
1818
1819     sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
1820     NAMD_die(err_msg);
1821   }
1822   else if ( (read_count != 7) && (paramType == paraCharmm) )
1823   {
1824     char err_msg[512];
1825
1826     sprintf(err_msg, "BAD IMPROPER FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
1827     NAMD_die(err_msg);
1828   }
1829
1830   if ( (read_count == 4) && (paramType == paraXplor) )
1831   //****** END CHARMM/XPLOR type changes
1832   {
1833     read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n", 
1834           &forceconstant, &periodicity, &phase_shift);
1835
1836     /*  Check to make sure we got what we expected    */
1837     if (read_count != 3)
1838     {
1839       char err_msg[512];
1840
1841       sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1842       NAMD_die(err_msg);
1843     }
1844
1845     multiplicity = 1;
1846   }
1847
1848   if (multiplicity > MAX_MULTIPLICITY)
1849   {
1850     char err_msg[181];
1851
1852     sprintf(err_msg, "Multiple improper with multiplicity of %d greater than max of %d",
1853        multiplicity, MAX_MULTIPLICITY);
1854     NAMD_die(err_msg);
1855   }
1856
1857   /*  Allocate a new node            */
1858   new_node = new improper_params;
1859
1860   if (new_node == NULL)
1861   {
1862     NAMD_die("memory allocation failed in Parameters::add_improper_param");
1863   }
1864
1865   /*  Assign the values for this bond.  As with the dihedrals,    */
1866   /*  the atom order doesn't matter        */
1867   strcpy(new_node->atom1name, atom1name);
1868   strcpy(new_node->atom2name, atom2name);
1869   strcpy(new_node->atom3name, atom3name);
1870   strcpy(new_node->atom4name, atom4name);
1871   new_node->multiplicity = multiplicity;
1872   if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1873   new_node->values[0].k = forceconstant;
1874   new_node->values[0].n = periodicity;
1875   new_node->values[0].delta = phase_shift;
1876
1877   new_node->next = NULL;
1878
1879   //  Check to see if this improper has multiple values
1880   if (multiplicity > 1)
1881   {
1882     //  Loop through and read the other values
1883     for (i=1; i<multiplicity; i++)
1884     {
1885       ret_code = NAMD_read_line(fd, buffer);
1886
1887       //  Strip off comments at the end of the line
1888       if (ret_code == 0)
1889       {
1890         NAMD_remove_comment(buffer);
1891       }
1892
1893       //  Skip blank lines
1894       while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
1895       {
1896         ret_code = NAMD_read_line(fd, buffer);
1897       }
1898
1899       if (ret_code != 0)
1900       {
1901         NAMD_die("EOF encoutner in middle of multiple improper");
1902       }
1903
1904       //  Get the values from the line
1905       read_count=sscanf(buffer, "%f %d %f\n", 
1906             &forceconstant, &periodicity, &phase_shift);
1907
1908       if (read_count != 3)
1909       {
1910         char err_msg[512];
1911
1912         sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
1913         NAMD_die(err_msg);
1914       }
1915
1916       if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1917       new_node->values[i].k = forceconstant;
1918       new_node->values[i].n = periodicity;
1919       new_node->values[i].delta = phase_shift;
1920     }
1921   }
1922
1923   /*  Add the paramter to the list        */
1924   add_to_improper_list(new_node);  // works for both XPLOR & CHARMM
1925
1926   return;
1927 }
1928 /*      END OF FUNCTION add_improper_param    */
1929
1930 /************************************************************************/
1931 /*                  */
1932 /*      FUNCTION add_to_improper_list      */
1933 /*                  */
1934 /*   INPUTS:                */
1935 /*  new_node - node that is to be added to imporper_list    */
1936 /*                  */
1937 /*  this function adds a new dihedral parameter to the linked list  */
1938 /*   of improper parameters.  First, it checks for duplicates.  If a    */
1939 /*   duplicate is found, a warning message is printed, the old values   */
1940 /*   are replaced with the new values, and the new node is freed.  If   */
1941 /*   Otherwise, the node is added to the list.  This list is arranged   */
1942 /*   so that bods with wildcards are placed at the tail of the list.    */
1943 /*   This will guarantee that if we just do a linear search, we will    */
1944 /*   always find an exact match before a wildcard match.    */
1945 /*                  */
1946 /************************************************************************/
1947
1948 void Parameters::add_to_improper_list(struct improper_params *new_node)
1949
1950 {
1951   int i;              //  Loop counter
1952   static struct improper_params *ptr;   //  position within list
1953   static struct improper_params *tail;  //  Pointer to the end of 
1954                 //  the list so we can add
1955                 //  entries to the end of the
1956                 //  list in constant time
1957
1958   /*  If the list is currently empty, then the new node is the list*/
1959   if (improperp == NULL)
1960   {
1961     improperp=new_node;
1962     tail=new_node;
1963
1964     return;
1965   }
1966
1967   /*  The list isn't empty, so check for a duplicate    */
1968   ptr=improperp;
1969
1970   while (ptr != NULL)
1971   {
1972     if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1973            (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1974            (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1975            (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1976          ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1977            (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1978            (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1979            (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
1980     {
1981       /*  Found a duplicate        */
1982       //****** BEGIN CHARMM/XPLOR type changes
1983       /* we do not care about identical replacement */
1984       int echoWarn=0;  // echo warning messages ?
1985
1986       if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
1987       
1988       if (!echoWarn)
1989       {
1990         for (i=0; i<ptr->multiplicity; i++)
1991         {
1992           if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
1993           if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
1994           if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
1995         }
1996       }
1997
1998       if (echoWarn)
1999       {
2000         iout << "\n" << iWARN << "DUPLICATE IMPROPER DIHEDRAL ENTRY FOR "
2001           << ptr->atom1name << "-"
2002           << ptr->atom2name << "-"
2003           << ptr->atom3name << "-"
2004           << ptr->atom4name
2005           << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
2006         
2007         for (i=0; i<ptr->multiplicity; i++)
2008         {
2009           iout <<     "  k=" << ptr->values[i].k
2010                    << "  n=" << ptr->values[i].n
2011                    << "  delta=" << ptr->values[i].delta;
2012         }
2013
2014         iout << "\n" << "USING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
2015
2016         for (i=0; i<new_node->multiplicity; i++)
2017         {
2018           iout <<     "  k=" << new_node->values[i].k
2019                    << "  n=" << new_node->values[i].n
2020                    << "  delta=" << new_node->values[i].delta;
2021         }
2022
2023         iout << endi;
2024
2025         ptr->multiplicity = new_node->multiplicity;
2026
2027         for (i=0; i<new_node->multiplicity; i++)
2028         {
2029           ptr->values[i].k = new_node->values[i].k;
2030           ptr->values[i].n = new_node->values[i].n;
2031           ptr->values[i].delta = new_node->values[i].delta;
2032         }
2033       }
2034       //****** END CHARMM/XPLOR type changes
2035
2036       delete new_node;
2037
2038       return;
2039     }
2040
2041     ptr=ptr->next;
2042   }
2043
2044   /*  Check to see if we have any wildcards.  Since specific  */
2045   /*  entries are to take precedence, we'll put anything without  */
2046   /*  wildcards at the begining of the list and anything with     */
2047   /*  wildcards at the end of the list.  Then, we can just do a   */
2048   /*  linear search for a bond and be guaranteed to have specific */
2049   /*  entries take precendence over over wildcards          */
2050   if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
2051        (strcasecmp(new_node->atom2name, "X") == 0) ||
2052        (strcasecmp(new_node->atom3name, "X") == 0) ||
2053        (strcasecmp(new_node->atom4name, "X") == 0) )
2054   {
2055     /*  add to the end of the list        */
2056     tail->next=new_node;
2057     tail=new_node;
2058
2059     return;
2060   }
2061   else
2062   {
2063     /*  add to the head of the list        */
2064     new_node->next=improperp;
2065     improperp=new_node;
2066
2067     return;
2068   }
2069 }
2070 /*    END OF FUNCTION add_to_improper_list      */
2071
2072 /************************************************************************/
2073 /*                  */
2074 /*      FUNCTION add_crossterm_param      */
2075 /*                  */
2076 /*   INPUTS:                */
2077 /*  buf - line from paramter file with crossterm parameters    */
2078 /*                  */
2079 /*  this function adds an crossterm parameter.  It parses up the     */
2080 /*   input line and then adds it to the binary tree used to store the   */
2081 /*   crossterm parameters.            */
2082 /*                  */
2083 /************************************************************************/
2084
2085 void Parameters::add_crossterm_param(char *buf, FILE *fd)
2086
2087 {
2088   char atom1name[11];       //  Atom 1 type
2089   char atom2name[11];       //  Atom 2 type
2090   char atom3name[11];       //  Atom 3 type
2091   char atom4name[11];       //  Atom 4 type
2092   char atom5name[11];       //  Atom 1 type
2093   char atom6name[11];       //  Atom 2 type
2094   char atom7name[11];       //  Atom 3 type
2095   char atom8name[11];       //  Atom 4 type
2096   int dimension;
2097   int read_count;         //  Count from sscanf
2098   struct crossterm_params *new_node;  //  New node
2099   char buffer[513];       //  Buffer for new line
2100   int ret_code;         //  Return code
2101
2102   /* read CHARMM format */
2103   read_count=sscanf(buf, "%s %s %s %s %s %s %s %s %d\n", 
2104      atom1name, atom2name, atom3name, atom4name,  
2105      atom5name, atom6name, atom7name, atom8name,  
2106      &dimension); 
2107
2108   if ( (read_count != 9) || dimension < 1 || dimension > 1000 )
2109   {
2110     char err_msg[512];
2111
2112     sprintf(err_msg, "BAD CMAP FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2113     NAMD_die(err_msg);
2114   }
2115
2116   /*  Allocate a new node            */
2117   new_node = new crossterm_params(dimension);
2118
2119   /*  Assign the values for this bond.  As with the dihedrals,    */
2120   /*  the atom order doesn't matter        */
2121   strcpy(new_node->atom1name, atom1name);
2122   strcpy(new_node->atom2name, atom2name);
2123   strcpy(new_node->atom3name, atom3name);
2124   strcpy(new_node->atom4name, atom4name);
2125   strcpy(new_node->atom5name, atom5name);
2126   strcpy(new_node->atom6name, atom6name);
2127   strcpy(new_node->atom7name, atom7name);
2128   strcpy(new_node->atom8name, atom8name);
2129
2130   new_node->next = NULL;
2131
2132   int nterms = dimension * dimension;
2133   int nread = 0;
2134
2135   //  Loop through and read the other values
2136   while ( nread < nterms ) {
2137     ret_code = NAMD_read_line(fd, buffer);
2138
2139     //  Strip off comments at the end of the line
2140     if (ret_code == 0) {
2141       NAMD_remove_comment(buffer);
2142     }
2143
2144     //  Skip blank lines
2145     while ( (ret_code == 0) && (NAMD_blank_string(buffer)) ) {
2146       ret_code = NAMD_read_line(fd, buffer);
2147       if (ret_code == 0) {
2148         NAMD_remove_comment(buffer);
2149       }
2150     }
2151
2152     if (ret_code != 0) {
2153       NAMD_die("EOF encoutner in middle of CMAP");
2154     }
2155
2156     //  Get the values from the line
2157     read_count=sscanf(buffer, "%lf %lf %lf %lf %lf %lf %lf %lf\n",
2158         new_node->values + nread,
2159         new_node->values + nread+1,
2160         new_node->values + nread+2,
2161         new_node->values + nread+3,
2162         new_node->values + nread+4,
2163         new_node->values + nread+5,
2164         new_node->values + nread+6,
2165         new_node->values + nread+7);
2166
2167     nread += read_count;
2168
2169     if (read_count == 0 || nread > nterms) {
2170       char err_msg[512];
2171
2172       sprintf(err_msg, "BAD CMAP FORMAT IN PARAMETER FILE\nLINE=*%s*\n", buffer);
2173       NAMD_die(err_msg);
2174     }
2175   }
2176
2177   /*  Add the paramter to the list        */
2178   add_to_crossterm_list(new_node);
2179
2180   return;
2181 }
2182 /*      END OF FUNCTION add_crossterm_param    */
2183
2184 /************************************************************************/
2185 /*                  */
2186 /*      FUNCTION add_to_crossterm_list      */
2187 /*                  */
2188 /*   INPUTS:                */
2189 /*  new_node - node that is to be added to crossterm_list    */
2190 /*                  */
2191 /*  this function adds a new crossterm parameter to the linked list  */
2192 /*   of crossterm parameters.  First, it checks for duplicates.  If a    */
2193 /*   duplicate is found, a warning message is printed, the old values   */
2194 /*   are replaced with the new values, and the new node is freed.  If   */
2195 /*   Otherwise, the node is added to the list.  This list is arranged   */
2196 /*   so that bods with wildcards are placed at the tail of the list.    */
2197 /*   This will guarantee that if we just do a linear search, we will    */
2198 /*   always find an exact match before a wildcard match.    */
2199 /*                  */
2200 /************************************************************************/
2201
2202 void Parameters::add_to_crossterm_list(struct crossterm_params *new_node)
2203
2204 {
2205   int i;              //  Loop counter
2206   static struct crossterm_params *ptr;   //  position within list
2207   static struct crossterm_params *tail;  //  Pointer to the end of 
2208                 //  the list so we can add
2209                 //  entries to the end of the
2210                 //  list in constant time
2211
2212   /*  If the list is currently empty, then the new node is the list*/
2213   if (crosstermp == NULL)
2214   {
2215     crosstermp=new_node;
2216     tail=new_node;
2217
2218     return;
2219   }
2220
2221   /*  The list isn't empty, so check for a duplicate    */
2222   ptr=crosstermp;
2223
2224   while (ptr != NULL)
2225   {
2226     if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
2227            (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
2228            (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
2229            (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) &&
2230            (strcasecmp(new_node->atom5name, ptr->atom5name) == 0) &&
2231            (strcasecmp(new_node->atom6name, ptr->atom6name) == 0) &&
2232            (strcasecmp(new_node->atom7name, ptr->atom7name) == 0) &&
2233            (strcasecmp(new_node->atom8name, ptr->atom8name) == 0) ) )
2234     {
2235       /*  Found a duplicate        */
2236       /* we do not care about identical replacement */
2237       int echoWarn=0;  // echo warning messages ?
2238
2239       if (ptr->dimension != new_node->dimension) {echoWarn=1;}
2240       
2241       if (!echoWarn)
2242       {
2243         int nvals = ptr->dimension * ptr->dimension;
2244         for (i=0; i<nvals; i++)
2245         {
2246           if (ptr->values[i] != new_node->values[i]) {echoWarn=1; break;}
2247         }
2248       }
2249
2250       if (echoWarn)
2251       {
2252         iout << "\n" << iWARN << "DUPLICATE CMAP ENTRY FOR "
2253           << ptr->atom1name << "-"
2254           << ptr->atom2name << "-"
2255           << ptr->atom3name << "-"
2256           << ptr->atom4name << " "
2257           << ptr->atom5name << "-"
2258           << ptr->atom6name << "-"
2259           << ptr->atom7name << "-"
2260           << ptr->atom8name << ", USING NEW VALUES\n";
2261         
2262         iout << endi;
2263
2264         ptr->dimension = new_node->dimension;
2265
2266         BigReal *tmpvalues = ptr->values;
2267         ptr->values = new_node->values;
2268         new_node->values = tmpvalues;
2269       }
2270
2271       delete new_node;
2272
2273       return;
2274     }
2275
2276     ptr=ptr->next;
2277   }
2278
2279   /*  Check to see if we have any wildcards.  Since specific  */
2280   /*  entries are to take precedence, we'll put anything without  */
2281   /*  wildcards at the begining of the list and anything with     */
2282   /*  wildcards at the end of the list.  Then, we can just do a   */
2283   /*  linear search for a bond and be guaranteed to have specific */
2284   /*  entries take precendence over over wildcards          */
2285   if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
2286        (strcasecmp(new_node->atom2name, "X") == 0) ||
2287        (strcasecmp(new_node->atom3name, "X") == 0) ||
2288        (strcasecmp(new_node->atom4name, "X") == 0) ||
2289        (strcasecmp(new_node->atom5name, "X") == 0) ||
2290        (strcasecmp(new_node->atom6name, "X") == 0) ||
2291        (strcasecmp(new_node->atom7name, "X") == 0) ||
2292        (strcasecmp(new_node->atom8name, "X") == 0) )
2293   {
2294     /*  add to the end of the list        */
2295     tail->next=new_node;
2296     tail=new_node;
2297
2298     return;
2299   }
2300   else
2301   {
2302     /*  add to the head of the list        */
2303     new_node->next=crosstermp;
2304     crosstermp=new_node;
2305
2306     return;
2307   }
2308 }
2309 /*    END OF FUNCTION add_to_crossterm_list      */
2310
2311 /************************************************************************/
2312 /*                  */
2313 /*      FUNCTION add_vdw_param        */
2314 /*                  */
2315 /*  INPUTS:                */
2316 /*  buf - line containing the vdw information      */
2317 /*                  */
2318 /*  add_vdw_param adds a vdw parameter for an atom to the current   */
2319 /*  binary tree of values.            */
2320 /*                  */
2321 /************************************************************************/
2322
2323 void Parameters::add_vdw_param(char *buf)
2324
2325 {
2326   char atomname[11];    //  atom type of paramter
2327   Real sigma;      //  sigma value for this atom
2328   Real epsilon;      //  epsilon value for this atom
2329   Real sigma14;      //  sigma value for 1-4 interactions
2330   Real epsilon14;      //  epsilon value for 1-4 interactions
2331   Real sqrt26;         //  2^(1/6)
2332   int read_count;      //  count returned by sscanf
2333   struct vdw_params *new_node;  //  new node for tree
2334
2335   //****** BEGIN CHARMM/XPLOR type changes
2336   /*  Parse up the line with sscanf        */
2337   if (paramType == paraXplor)
2338   {
2339     /* read XPLOR format */
2340     read_count=sscanf(buf, "%*s %s %f %f %f %f\n", atomname, 
2341        &epsilon, &sigma, &epsilon14, &sigma14);
2342   }
2343   else if (paramType == paraCharmm)
2344   {
2345     /* read CHARMM format */
2346     read_count=sscanf(buf, "%s %*f %f %f %*f %f %f\n", atomname, 
2347        &epsilon, &sigma, &epsilon14, &sigma14);
2348   }
2349
2350   /*  Check to make sure we got what we expected      */
2351   if ((read_count != 5) && (paramType == paraXplor))
2352   {
2353     char err_msg[512];
2354
2355     sprintf(err_msg, "BAD vdW FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
2356     NAMD_die(err_msg);
2357   }
2358   else if ( ((read_count != 5) && (read_count != 3)) && (paramType == paraCharmm))
2359   {
2360     char err_msg[512];
2361
2362     sprintf(err_msg, "BAD vdW FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
2363     NAMD_die(err_msg);
2364   }
2365
2366   if (paramType == paraCharmm)
2367   {
2368     // convert CHARMM to XPLOR format
2369     epsilon*=-1.;
2370     sqrt26=pow(2.,(1./6.));
2371     sigma=2.*sigma/sqrt26; 
2372
2373     if (read_count == 3)
2374     {
2375       epsilon14=epsilon;
2376       sigma14=sigma;
2377     }
2378     else
2379     {
2380       epsilon14*=-1.;
2381       sigma14=2.*sigma14/sqrt26; 
2382     }
2383   }
2384   //****** END CHARMM/XPLOR type changes
2385
2386   if ( epsilon < 0. || epsilon14 < 0. ) {
2387     iout << iWARN << "Ignoring VDW parameter with negative epsilon:\n"
2388         << buf << "\n" << endi;
2389     return;
2390   }
2391
2392   /*  Allocate a new node            */
2393   new_node = new vdw_params;
2394
2395   if (new_node == NULL)
2396   {
2397     NAMD_die("memory allocation failed in Parameters::add_vdw_param");
2398   }
2399
2400   /*  Assign the values to the new node        */
2401   strcpy(new_node->atomname, atomname);
2402   new_node->sigma = sigma;
2403   new_node->sigma14 = sigma14;
2404   new_node->epsilon = epsilon;
2405   new_node->epsilon14 = epsilon14;
2406
2407   new_node->left = NULL;
2408   new_node->right = NULL;
2409
2410   /*  Add the new node into the tree        */
2411   vdwp=add_to_vdw_tree(new_node, vdwp);
2412
2413   return;
2414 }
2415 /*      END OF FUNCTION add_vdw_param      */
2416
2417 /************************************************************************/
2418 /*                  */
2419 /*      FUNCTION add_to_vdw_tree      */
2420 /*                  */
2421 /*   INPUTS:                */
2422 /*  new_node - node to add to tree          */
2423 /*  tree - tree to add the node to          */
2424 /*                  */
2425 /*   OUTPUTS:                */
2426 /*  the function returns a pointer to the tree with the node added  */
2427 /*                  */
2428 /*  this function adds a vdw to the binary tree containing the      */
2429 /*   parameters.              */
2430 /*                  */
2431 /************************************************************************/
2432
2433 struct vdw_params *Parameters::add_to_vdw_tree(struct vdw_params *new_node,
2434              struct vdw_params *tree)
2435
2436 {
2437   int compare_code;  //  Return code from strcasecmp
2438
2439   /*  If the tree is currently empty, the new node is the tree    */
2440   if (tree == NULL)
2441     return(new_node);
2442
2443   compare_code = strcasecmp(new_node->atomname, tree->atomname);
2444
2445   /*  Check to see if we have a duplicate        */
2446   if (compare_code==0)
2447   {
2448     /*  We have a duplicate.  So print out a warning   */
2449     /*  message, copy the new values into the current node  */
2450     /*  of the tree, and then free the new_node    */
2451     if ((tree->sigma != new_node->sigma) || 
2452         (tree->epsilon != new_node->epsilon) ||
2453         (tree->sigma14 != new_node->sigma14) ||
2454         (tree->epsilon14 != new_node->epsilon14))
2455     {
2456       iout << iWARN << "DUPLICATE vdW ENTRY FOR " << tree->atomname
2457         << "\nPREVIOUS VALUES  sigma=" << tree->sigma
2458         << " epsilon=" << tree->epsilon
2459         << " sigma14=" << tree->sigma14
2460         << " epsilon14=" << tree->epsilon14
2461         << "\n   USING VALUES  sigma=" << new_node->sigma
2462         << " epsilon=" << new_node->epsilon
2463         << " sigma14=" << new_node->sigma14
2464         << " epsilon14=" << new_node->epsilon14
2465         << "\n" << endi;
2466
2467       tree->sigma=new_node->sigma;
2468       tree->epsilon=new_node->epsilon;
2469       tree->sigma14=new_node->sigma14;
2470       tree->epsilon14=new_node->epsilon14;
2471     }
2472
2473     delete new_node;
2474
2475     return(tree);
2476   }
2477
2478   /*  Otherwise, if the new node is less than the head of    */
2479   /*  the tree, add it to the left child, and if it is greater  */
2480   /*  add it to the right child          */
2481   if (compare_code < 0)
2482   {
2483     tree->left = add_to_vdw_tree(new_node, tree->left);
2484   }
2485   else
2486   {
2487     tree->right = add_to_vdw_tree(new_node, tree->right);
2488   }
2489
2490   return(tree);
2491 }
2492 /*      END OF FUNCTION add_to_vdw_tree      */
2493
2494 /************************************************************************/
2495 /*                  */
2496 /*      FUNCTION add_table_pair_param      */
2497 /*                  */
2498 /*   INPUTS:                */
2499 /*  buf - line containing the table_pair information      */
2500 /*                  */
2501 /*  this function adds a table_pair parameter to the current          */
2502 /*   parameters.              */
2503 /*                  */
2504 /************************************************************************/
2505
2506 void Parameters::add_table_pair_param(char *buf)
2507
2508 {
2509   char atom1name[11];      //  Atom 1 name
2510   char atom2name[11];      //  Atom 2 name
2511   char tabletype[11];      // Name of pair interaction
2512   int K;           // Table entry type for pair
2513   int read_count;        //  count from sscanf
2514   struct table_pair_params *new_node;  //  new node
2515
2516   /*  Parse up the input line using sscanf      */
2517   read_count=sscanf(buf, "%s %s %s\n", atom1name, 
2518   atom2name, tabletype);
2519
2520   /*  Check to make sure we got what we expected      */
2521   if ((read_count != 3))
2522   {
2523     char err_msg[512];
2524
2525     sprintf(err_msg, "BAD TABLE PAIR FORMAT IN PARAMETER FILE\nLINE=*%s*", buf);
2526     NAMD_die(err_msg);
2527   }
2528
2529   /*  Allocate a new node            */
2530   new_node = new table_pair_params;
2531
2532   if (new_node == NULL)
2533   {
2534     NAMD_die("memory allocation failed in Parameters::add_table_pair_param\n");
2535   }
2536
2537   strcpy(new_node->atom1name, atom1name);
2538   strcpy(new_node->atom2name, atom2name);
2539
2540   /* Find the proper table type for this pairing */
2541   K = get_int_table_type(tabletype);
2542 //  printf("Looking for type %s; got %i\n", tabletype, K);
2543   if (K < 0) {
2544     char err_msg[512];
2545     sprintf(err_msg, "Couldn't find table parameters for table interaction %s!\n", tabletype);
2546     NAMD_die(err_msg);
2547   }
2548
2549   /*  Assign values to this node          */
2550   new_node->K = K;
2551
2552   new_node->next = NULL;
2553
2554   /*  Add this node to the tree          */
2555 //  printf("Adding pair parameter with index %i\n", K);
2556   add_to_table_pair_list(new_node);
2557
2558   return;
2559 }
2560 /*      END OF FUNCTION add_vdw_par_param    */
2561
2562 /************************************************************************/
2563 /*                  */
2564 /*      FUNCTION add_vdw_pair_param      */
2565 /*                  */
2566 /*   INPUTS:                */
2567 /*  buf - line containing the vdw_pair information      */
2568 /*                  */
2569 /*  this function adds a vdw_pair parameter to the current          */
2570 /*   parameters.              */
2571 /*                  */
2572 /************************************************************************/
2573
2574 void Parameters::add_vdw_pair_param(char *buf)
2575
2576 {
2577   char atom1name[11];      //  Atom 1 name
2578   char atom2name[11];      //  Atom 2 name
2579   Real A;          //  A value for pair
2580   Real B;          //  B value for pair
2581   Real A14;        //  A value for 1-4 ints
2582   Real B14;        //  B value for 1-4 ints
2583   int read_count;        //  count from sscanf
2584   struct vdw_pair_params *new_node;  //  new node
2585
2586   /*  Parse up the input line using sscanf      */
2587   if (paramType == paraXplor)
2588   {
2589     /* read XPLOR format */
2590     read_count=sscanf(buf, "%*s %s %s %f %f %f %f\n", atom1name, 
2591        atom2name, &A, &B, &A14, &B14);
2592   }
2593   else if (paramType == paraCharmm)
2594   {
2595     Real well, rmin, well14, rmin14;
2596     /* read CHARMM format */
2597     read_count=sscanf(buf, "%s %s %f %f %f %f\n", atom1name, 
2598        atom2name, &well, &rmin, &well14, &rmin14);
2599     if ( read_count == 4 ) { well14 = well; rmin14 = rmin; }
2600     A = -1. * well * pow(rmin, 12.);
2601     B = -2. * well * pow(rmin, 6.);
2602     A14 = -1. * well14 * pow(rmin14, 12.);
2603     B14 = -2. * well14 * pow(rmin14, 6.);
2604   }
2605
2606   /*  Check to make sure we got what we expected      */
2607   if ((read_count != 6) && (paramType == paraXplor))
2608   {
2609     char err_msg[512];
2610
2611     sprintf(err_msg, "BAD vdW PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
2612     NAMD_die(err_msg);
2613   }
2614   if ((read_count != 4) && (read_count != 6) && (paramType == paraCharmm))
2615   {
2616     char err_msg[512];
2617
2618     sprintf(err_msg, "BAD vdW PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2619     NAMD_die(err_msg);
2620   }
2621
2622
2623   /*  Allocate a new node            */
2624   new_node = new vdw_pair_params;
2625
2626   if (new_node == NULL)
2627   {
2628     NAMD_die("memory allocation failed in Parameters::add_vdw_pair_param\n");
2629   }
2630
2631   strcpy(new_node->atom1name, atom1name);
2632   strcpy(new_node->atom2name, atom2name);
2633
2634   /*  Assign values to this node          */
2635   new_node->A = A;
2636   new_node->A14 = A14;
2637   new_node->B = B;
2638   new_node->B14 = B14;
2639
2640   new_node->next = NULL;
2641
2642   /*  Add this node to the tree          */
2643   add_to_vdw_pair_list(new_node);
2644
2645   return;
2646 }
2647 /*      END OF FUNCTION add_vdw_par_param    */
2648
2649 /************************************************************************/
2650 /*                  */
2651 /*      FUNCTION add_nbthole_pair_param      */
2652 /*                  */
2653 /*   INPUTS:                */
2654 /*  buf - line containing the nbthole_pair information      */
2655 /*                  */
2656 /*  this function adds a nbthole_pair parameter to the current          */
2657 /*   parameters.              */
2658 /*                  */
2659 /************************************************************************/
2660
2661 void Parameters::add_nbthole_pair_param(char *buf)
2662
2663 {
2664   char atom1name[11];      //  Atom 1 name
2665   char atom2name[11];      //  Atom 2 name
2666   Real tholeij;            //  nonbonded thole pair thole
2667   int read_count;        //  count from sscanf
2668   struct nbthole_pair_params *new_node;  //  new node
2669
2670   /*  Parse up the input line using sscanf      */
2671   if (paramType == paraCharmm)
2672   {
2673     /* read CHARMM format */
2674     read_count=sscanf(buf, "%s %s %f\n", atom1name,
2675        atom2name, &tholeij);
2676   }
2677
2678   /*  Check to make sure we got what we expected      */
2679   if ((read_count != 3) && (paramType == paraCharmm))
2680   {
2681     char err_msg[512];
2682
2683     sprintf(err_msg, "BAD NBTHOLE PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2684     NAMD_die(err_msg);
2685   }
2686
2687
2688   /*  Allocate a new node            */
2689   new_node = new nbthole_pair_params;
2690
2691   if (new_node == NULL)
2692   {
2693     NAMD_die("memory allocation failed in Parameters::nbthole_pair_param\n");
2694   }
2695
2696   strcpy(new_node->atom1name, atom1name);
2697   strcpy(new_node->atom2name, atom2name);
2698
2699   /*  Assign values to this node          */
2700   new_node->tholeij = tholeij;
2701
2702   new_node->next = NULL;
2703
2704   /*  Add this node to the tree          */
2705   add_to_nbthole_pair_list(new_node);
2706
2707   return;
2708 }
2709 /*      END OF FUNCTION add_nbthole_par_param    */
2710
2711 /************************************************************************/
2712 /*                  */
2713 /*      FUNCTION add_hb_pair_param      */
2714 /*                  */
2715 /*   INPUTS:                */
2716 /*  buf - line containing the hydrogen bond information    */
2717 /*                  */
2718 /*  this function adds data for a hydrogen bond interaction pair    */
2719 /*   to the hbondParams object.                                         */
2720 /*                  */
2721 /************************************************************************/
2722
2723 void Parameters::add_hb_pair_param(char *buf)
2724
2725 {
2726 #if 0
2727   char a1n[11];      //  Atom 1 name
2728   char a2n[11];      //  Atom 2 name
2729   Real A, B;      //  A, B value for pair
2730
2731   //****** BEGIN CHARMM/XPLOR type changes
2732   //// luckily the format and units are the same CHARMM is just missing the HBON marker
2733   /*  Parse up the input line using sscanf      */
2734   if (paramType == paraXplor) {
2735     if (sscanf(buf, "%*s %s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
2736       char err_msg[512];
2737       sprintf(err_msg, "BAD HBOND PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
2738       NAMD_die(err_msg);
2739     }
2740   }
2741   else if (paramType == paraCharmm) {
2742     if (sscanf(buf, "%s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
2743       char err_msg[512];
2744       sprintf(err_msg, "BAD HBOND PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2745       NAMD_die(err_msg);
2746     }
2747   }
2748   //****** END CHARMM/XPLOR type changes
2749
2750   /*  add data */
2751   if (hbondParams.add_hbond_pair(a1n, a2n, A, B) == FALSE) {
2752     iout << "\n" << iWARN << "Duplicate HBOND parameters for types " << a1n
2753     << " and " << a2n << " found; using latest values." << "\n" << endi;
2754   }
2755 #endif
2756 }
2757 /*      END OF FUNCTION add_hb_par_param    */
2758
2759 /************************************************************************/
2760 /*                  */
2761 /*      FUNCTION add_to_table_pair_list      */
2762 /*                  */
2763 /*   INPUTS:                */
2764 /*  new_node - node to be added to list        */
2765 /*                  */
2766 /*  This function adds a link to the end of the table_pair_list list  */
2767 /*                  */
2768 /************************************************************************/
2769
2770 void Parameters::add_to_table_pair_list(struct table_pair_params *new_node)
2771
2772 {
2773      static struct table_pair_params *tail=NULL;
2774   struct table_pair_params *ptr;
2775   int compare_code;
2776   
2777
2778   //  If the list was empty, then just make the new node the list
2779   if (table_pairp == NULL)
2780   {
2781      table_pairp = new_node;
2782      tail = new_node;
2783      return;
2784   }
2785   
2786   ptr = table_pairp;
2787
2788   //  Now check the list to see if we have a duplicate entry
2789   while (ptr!=NULL)
2790   {
2791       /*  Compare atom 1            */
2792       compare_code = strncasecmp(new_node->atom1name, ptr->atom1name, 5);
2793       
2794       if (compare_code == 0)
2795       {
2796     /*  Atom 1 is the same, compare atom 2      */
2797     compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
2798
2799     if (compare_code==0)
2800     {
2801       /*  Found a duplicate.  Print out a warning   */
2802       /*  message, assign the values to the current   */
2803       /*  node in the tree, and then free the new_node*/
2804       iout << iWARN << "DUPLICATE TABLE PAIR ENTRY FOR "
2805         << new_node->atom1name << "-"
2806         << new_node->atom2name
2807         << "\n" << endi;
2808
2809         ptr->K=new_node->K;
2810
2811       delete new_node;
2812
2813       return;
2814     }
2815       }
2816       
2817       ptr = ptr->next;
2818   }
2819
2820   //  We didn't find a duplicate, so add this node to the end
2821   //  of the list
2822   tail->next = new_node;
2823   tail = new_node;
2824 }
2825 /*      END OF FUNCTION add_to_vdw_pair_list    */
2826
2827 /************************************************************************/
2828 /*                  */
2829 /*      FUNCTION add_to_vdw_pair_list      */
2830 /*                  */
2831 /*   INPUTS:                */
2832 /*  new_node - node to be added to list        */
2833 /*                  */
2834 /*  This function adds a link to the end of the vdw_pair_list list  */
2835 /*                  */
2836 /************************************************************************/
2837
2838 void Parameters::add_to_vdw_pair_list(struct vdw_pair_params *new_node)
2839
2840 {
2841      static struct vdw_pair_params *tail=NULL;
2842   struct vdw_pair_params *ptr;
2843   int compare_code;
2844   
2845
2846   //  If the list was empty, then just make the new node the list
2847   if (vdw_pairp == NULL)
2848   {
2849      vdw_pairp = new_node;
2850      tail = new_node;
2851      return;
2852   }
2853   
2854   ptr = vdw_pairp;
2855
2856   //  Now check the list to see if we have a duplicate entry
2857   while (ptr!=NULL)
2858   {
2859       /*  Compare atom 1            */
2860       compare_code = strcasecmp(new_node->atom1name, ptr->atom1name);
2861       
2862       if (compare_code == 0)
2863       {
2864     /*  Atom 1 is the same, compare atom 2      */
2865     compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
2866
2867     if (compare_code==0)
2868     {
2869       /*  Found a duplicate.  Print out a warning   */
2870       /*  message, assign the values to the current   */
2871       /*  node in the tree, and then free the new_node*/
2872       if ((ptr->A != new_node->A) ||
2873           (ptr->B != new_node->B) ||
2874           (ptr->A14 != new_node->A14) ||
2875           (ptr->B14 != new_node->B14))
2876       {
2877         iout << iWARN << "DUPLICATE vdW PAIR ENTRY FOR "
2878           << new_node->atom1name << "-"
2879           << new_node->atom2name
2880           << "\nPREVIOUS VALUES  A=" << ptr->A
2881           << " B=" << ptr->B
2882           << " A14=" << ptr->A14
2883           << " B14" << ptr->B14
2884           << "\n   USING VALUES  A=" << new_node->A
2885           << " B=" << new_node->B
2886           << " A14=" << new_node->A14
2887           << " B14" << new_node->B14
2888           << "\n" << endi;
2889
2890         ptr->A=new_node->A;
2891         ptr->B=new_node->B;
2892         ptr->A14=new_node->A14;
2893         ptr->B14=new_node->B14;
2894       }
2895
2896       delete new_node;
2897
2898       return;
2899     }
2900       }
2901       
2902       ptr = ptr->next;
2903   }
2904
2905   //  We didn't find a duplicate, so add this node to the end
2906   //  of the list
2907   tail->next = new_node;
2908   tail = new_node;
2909 }
2910 /*      END OF FUNCTION add_to_vdw_pair_list    */
2911
2912 /************************************************************************/
2913 /*                  */
2914 /*      FUNCTION add_to_nbthole_pair_list      */
2915 /*                  */
2916 /*   INPUTS:                */
2917 /*  new_node - node to be added to list        */
2918 /*                  */
2919 /*  This function adds a link to the end of the nbthole_pair_list list  */
2920 /*                  */
2921 /************************************************************************/
2922
2923 void Parameters::add_to_nbthole_pair_list(struct nbthole_pair_params *new_node)
2924
2925 {
2926      static struct nbthole_pair_params *tail=NULL;
2927   struct nbthole_pair_params *ptr;
2928   int compare_code;
2929
2930
2931   //  If the list was empty, then just make the new node the list
2932   if (nbthole_pairp == NULL)
2933   {
2934      nbthole_pairp = new_node;
2935      tail = new_node;
2936      return;
2937   }
2938
2939   ptr = nbthole_pairp;
2940
2941   tail->next = new_node;
2942   tail = new_node;
2943 }
2944 /*      END OF FUNCTION add_to_nbthole_pair_list    */
2945
2946 /************************************************************************/
2947 /*                  */
2948 /*      FUNCTION done_reading_files      */
2949 /*                  */
2950 /*  This function is used to signal the Parameters object that all  */
2951 /*  of the parameter files have been read.  Once the object knows this, */
2952 /*  it can set un indexes for all the parameters and transfer the values*/
2953 /*  to linear arrays.  This will allow constant time access from this   */
2954 /*  point on.                */
2955 /*                  */
2956 /************************************************************************/
2957
2958 void Parameters::done_reading_files()
2959
2960 {
2961   AllFilesRead = TRUE;
2962
2963   //  Allocate space for all of the arrays
2964   if (NumBondParams)
2965   {
2966     bond_array = new BondValue[NumBondParams];
2967
2968     if (bond_array == NULL)
2969     {
2970       NAMD_die("memory allocation of bond_array failed!");
2971     }
2972     memset(bond_array, 0, NumBondParams*sizeof(BondValue));
2973   }
2974
2975   if (NumAngleParams)
2976   {
2977     angle_array = new AngleValue[NumAngleParams];
2978
2979     if (angle_array == NULL)
2980     {
2981       NAMD_die("memory allocation of angle_array failed!");
2982     }
2983     memset(angle_array, 0, NumAngleParams*sizeof(AngleValue));
2984     for ( Index i=0; i<NumAngleParams; ++i ) {
2985       angle_array[i].normal = 1;
2986     }
2987   }
2988
2989   if (NumDihedralParams)
2990   {
2991     dihedral_array = new DihedralValue[NumDihedralParams];
2992
2993     if (dihedral_array == NULL)
2994     {
2995       NAMD_die("memory allocation of dihedral_array failed!");
2996     }
2997     memset(dihedral_array, 0, NumDihedralParams*sizeof(DihedralValue));
2998   }
2999
3000   if (NumImproperParams)
3001   {
3002     improper_array = new ImproperValue[NumImproperParams];
3003
3004     if (improper_array == NULL)
3005     {
3006       NAMD_die("memory allocation of improper_array failed!");
3007     }
3008     memset(improper_array, 0, NumImproperParams*sizeof(ImproperValue));
3009   }
3010
3011   if (NumCrosstermParams)
3012   {
3013     crossterm_array = new CrosstermValue[NumCrosstermParams];
3014     memset(crossterm_array, 0, NumCrosstermParams*sizeof(CrosstermValue));
3015   }
3016
3017   // JLai
3018   if (NumGromacsPairParams)
3019   {
3020     gromacsPair_array = new GromacsPairValue[NumGromacsPairParams];
3021     memset(gromacsPair_array, 0, NumGromacsPairParams*sizeof(GromacsPairValue)); 
3022   }
3023   // End of JLai
3024
3025   if (NumVdwParams)
3026   {
3027           atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
3028     vdw_array = new VdwValue[NumVdwParams];
3029     
3030     if (vdw_array == NULL)
3031     {
3032       NAMD_die("memory allocation of vdw_array failed!");
3033     }
3034   }
3035   if (NumNbtholePairParams)
3036   {
3037     nbthole_array = new NbtholePairValue[NumNbtholePairParams];
3038
3039     if(nbthole_array == NULL)
3040     {
3041       NAMD_die("memory allocation of nbthole_array failed!");
3042     }
3043   }
3044   //  Assign indexes to each of the parameters and populate the
3045   //  arrays using the binary trees and linked lists that we have
3046   //  already read in
3047
3048   // Note that if parameters have been overwritten (matching
3049   // atom patterns but different parameter values) the tree
3050   // contains fewer elements than Num...Params would suggest.
3051   // The arrays are initialized above because the end values
3052   // may not be occupied.  Modifying the Num...Params values
3053   // would break backwards compatibility of memopt extraBonds.
3054
3055   index_bonds(bondp, 0);
3056   index_angles(anglep, 0);
3057   NumVdwParamsAssigned = index_vdw(vdwp, 0);
3058   index_dihedrals();
3059   index_impropers();
3060   index_crossterms();
3061   convert_nbthole_pairs();
3062   //  Convert the vdw pairs
3063   convert_vdw_pairs();
3064   convert_table_pairs();
3065 }
3066 /*      END OF FUNCTION done_reading_files    */
3067
3068 /************************************************************************/
3069 /*                  */
3070 /*      FUNCTION index_bonds        */
3071 /*                  */
3072 /*   INPUTS:                */
3073 /*  tree - The tree that is to be indexed        */
3074 /*  index - index to start with          */
3075 /*                  */
3076 /*  This is a recursive routine that will traverse the binary tree  */
3077 /*   of bond parameters, assigning an index to each one, and copying    */
3078 /*   the data from the binary tree to the array that will be used from  */
3079 /*   here on.                */
3080 /*                  */
3081 /************************************************************************/
3082
3083 Index Parameters::index_bonds(struct bond_params *tree, Index index)
3084
3085 {
3086   //  Tree is empty, do nothing
3087   if (tree==NULL)
3088     return(index);
3089
3090   //  If I have a left subtree, index it first
3091   if (tree->left != NULL)
3092   {
3093     index=index_bonds(tree->left, index);
3094   }
3095
3096   //  Now assign an index to top node and populate array
3097   tree->index = index;
3098   bond_array[index].k = tree->forceconstant;
3099   bond_array[index].x0 = tree->distance;
3100   index++;
3101
3102   //  If I have a right subtree, index it
3103   if (tree->right != NULL)
3104   {
3105     index=index_bonds(tree->right, index);
3106   }
3107
3108   return(index);
3109 }
3110 /*      END OF FUNCTION index_bonds      */
3111
3112 /************************************************************************/
3113 /*                  */
3114 /*      FUNCTION index_angles        */
3115 /*                  */
3116 /*   INPUTS:                */
3117 /*  tree - The tree that is to be indexed        */
3118 /*  index - index to start with          */
3119 /*                  */
3120 /*  This is a recursive routine that will traverse the binary tree  */
3121 /*   of angle parameters, assigning an index to each one, and copying   */
3122 /*   the data from the binary tree to the array that will be used from  */
3123 /*   here on.                */
3124 /*                  */
3125 /************************************************************************/
3126
3127 Index Parameters::index_angles(struct angle_params *tree, Index index)
3128
3129 {
3130   //  Tree is empty, do nothing
3131   if (tree==NULL)
3132     return(index);
3133
3134   //  If I have a left subtree, index it first
3135   if (tree->left != NULL)
3136   {
3137     index=index_angles(tree->left, index);
3138   }
3139
3140   //  Now assign an index to top node and populate array
3141   tree->index = index;
3142
3143   angle_array[index].k = tree->forceconstant;
3144   angle_array[index].k_ub = tree->k_ub;
3145   angle_array[index].r_ub = tree->r_ub;
3146   angle_array[index].normal = tree->normal;
3147
3148   //  Convert the angle to radians before storing it
3149   angle_array[index].theta0 = (tree->angle*PI)/180.0;
3150   index++;
3151
3152   //  If I have a right subtree, index it
3153   if (tree->right != NULL)
3154   {
3155     index=index_angles(tree->right, index);
3156   }
3157
3158   return(index);
3159 }
3160 /*      END OF FUNCTION index_angles      */
3161
3162 /************************************************************************/
3163 /*                  */
3164 /*      FUNCTION index_dihedrals      */
3165 /*                  */
3166 /*  This function walks down the linked list of dihedral parameters */
3167 /*  and assigns an index to each one.  It also copies the data from this*/
3168 /*  linked list to the arrays that will be used from here on out  */
3169 /*                  */
3170 /************************************************************************/
3171
3172 void Parameters::index_dihedrals()
3173
3174 {
3175   struct dihedral_params *ptr;  //  Current location in list
3176   Index index=0;      //  Current index value
3177   int i;        //  Loop counter
3178
3179   //  Allocate an array to hold the multiplicity present in the
3180   //  parameter file for each bond.  This will be used to check
3181   //  the multiplicities that are detected in the psf file
3182
3183   //  This is kind of ugly, but necessary because of the way that
3184   //  X-PLOR psf files deal with Charmm22 parameters.  The way
3185   //  that multiple periodicities are specified is by having
3186   //  the bonds appear multiple times in the psf file.  This even
3187   //  if a bond type has multiple parameters defined, they
3188   //  will be used if the bond appears multiple times in the
3189   //  psf file.  So we need to store the number of parameters
3190   //  we have to make sure the psf file doesn't ask for more
3191   //  parameters than we really have, and we also need to track
3192   //  how many times the bond appears in the psf file so that
3193   //  we can decide how many parameters to actually use.
3194   //  This is different for CHARMM parameter files as stated below!
3195   maxDihedralMults = new int[NumDihedralParams];
3196
3197   if (maxDihedralMults == NULL)
3198   {
3199     NAMD_die("memory allocation failed in Parameters::index_dihedrals()");
3200   }
3201   
3202   //  Start at the head
3203   ptr = dihedralp;
3204
3205   while (ptr != NULL)
3206   {
3207     //  Copy data to array and assign index
3208
3209     //  Save the multiplicity in another array
3210     maxDihedralMults[index] = ptr->multiplicity;
3211
3212
3213     //****** BEGIN CHARMM/XPLOR type changes
3214     if (paramType == paraXplor)
3215     {
3216       //  Assign the multiplicity in the actual structure a bogus value
3217       //  that we will update in assign_dihedral_index
3218       dihedral_array[index].multiplicity = -1;
3219     }
3220     else if (paramType == paraCharmm)
3221     {
3222       // In a CHARMM psf file each dihedral will be only listed once
3223       // even if it has multiple terms. There is no point in comparing
3224       // to the psf information
3225       dihedral_array[index].multiplicity = ptr->multiplicity;
3226     } 
3227     //****** END CHARMM/XPLOR type changes
3228
3229     for (i=0; i<ptr->multiplicity; i++)
3230     {
3231       dihedral_array[index].values[i].k = ptr->values[i].k;
3232       dihedral_array[index].values[i].n = ptr->values[i].n;
3233
3234       //  Convert the angle to radians before storing it
3235       dihedral_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
3236     }
3237
3238     ptr->index = index;
3239
3240     index++;
3241     ptr=ptr->next;
3242   }
3243 }
3244 /*      END OF FUNCTION index_dihedrals      */
3245
3246 /************************************************************************/
3247 /*                  */
3248 /*      FUNCTION index_impropers      */
3249 /*                  */
3250 /*  This function walks down the linked list of improper parameters */
3251 /*  and assigns an index to each one.  It also copies the data from this*/
3252 /*  linked list to the arrays that will be used from here on out  */
3253 /*                  */
3254 /************************************************************************/
3255
3256 void Parameters::index_impropers()
3257
3258 {
3259   struct improper_params *ptr;  //  Current place in list
3260   Index index=0;      //  Current index value
3261   int i;        //  Loop counter
3262
3263   //  Allocate an array to hold the multiplicity present in the
3264   //  parameter file for each bond.  This will be used to check
3265   //  the multiplicities that are detected in the psf file
3266
3267   //  This is kind of ugly, but necessary because of the way that
3268   //  X-PLOR psf files deal with Charmm22 parameters.  The way
3269   //  that multiple periodicities are specified is by having
3270   //  the bonds appear multiple times in the psf file.  This even
3271   //  if a bond type has multiple parameters defined, they
3272   //  will be used if the bond appears multiple times in the
3273   //  psf file.  So we need to store the number of parameters
3274   //  we have to make sure the psf file doesn't ask for more
3275   //  parameters than we really have, and we also need to track
3276   //  how many times the bond appears in the psf file so that
3277   //  we can decide how many parameters to actually use.
3278   maxImproperMults = new int[NumImproperParams];
3279
3280   if (maxImproperMults == NULL)
3281   {
3282     NAMD_die("memory allocation failed in Parameters::index_impropers()");
3283   }
3284   
3285   //  Start at the head
3286   ptr = improperp;
3287
3288   while (ptr != NULL)
3289   {
3290     //  Copy data to array and assign index
3291
3292     //  Save the multiplicity in another array
3293     maxImproperMults[index] = ptr->multiplicity;
3294
3295
3296     //****** BEGIN CHARMM/XPLOR type changes
3297     if (paramType == paraXplor)
3298     {
3299       //  Assign the multiplicity in the actual structure a bogus value
3300       //  that we will update in assign_improper_index
3301       improper_array[index].multiplicity = -1;
3302     }
3303     else if (paramType == paraCharmm)
3304     {
3305       // In a CHARMM psf file each improper will be only listed once
3306       // even if it has multiple terms. There is no point in comparing
3307       // to the psf information
3308       improper_array[index].multiplicity = ptr->multiplicity;
3309     } 
3310     //****** END CHARMM/XPLOR type changes
3311
3312     for (i=0; i<ptr->multiplicity; i++)
3313     {
3314       improper_array[index].values[i].k = ptr->values[i].k;
3315       improper_array[index].values[i].n = ptr->values[i].n;
3316
3317       //  Convert the angle to radians before storing it
3318       improper_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
3319     }
3320
3321     ptr->index=index;
3322
3323     index++;
3324     ptr=ptr->next;
3325   }
3326 }
3327 /*      END OF FUNCTION index_impropers      */
3328
3329
3330 /************************************************************************/
3331 /*                  */
3332 /*      FUNCTION index_crossterms      */
3333 /*                  */
3334 /*  This function walks down the linked list of crossterm parameters */
3335 /*  and assigns an index to each one.  It also copies the data from this*/
3336 /*  linked list to the arrays that will be used from here on out  */
3337 /*                  */
3338 /************************************************************************/
3339
3340 void crossterm_setup(CrosstermData *);
3341
3342 void Parameters::index_crossterms()
3343
3344 {
3345   struct crossterm_params *ptr;  //  Current place in list
3346   Index index=0;      //  Current index value
3347   int i,j,k;        //  Loop counter
3348
3349   //  Start at the head
3350   ptr = crosstermp;
3351
3352   while (ptr != NULL)
3353   {
3354     //  Copy data to array and assign index
3355
3356     int N = CrosstermValue::dim - 1;
3357
3358     if ( ptr->dimension != N ) {
3359       NAMD_die("Sorry, only CMAP dimension of 24 is supported");
3360     }
3361
3362     k = 0;
3363     for (i=0; i<N; i++) {
3364       for (j=0; j<N; j++) {
3365         crossterm_array[index].c[i][j].d00 = ptr->values[k];
3366         ++k;
3367       }
3368     }
3369     for (i=0; i<N; i++) {
3370         crossterm_array[index].c[i][N].d00 = 
3371                                 crossterm_array[index].c[i][0].d00;
3372         crossterm_array[index].c[N][i].d00 = 
3373                                 crossterm_array[index].c[0][i].d00;
3374     }
3375     crossterm_array[index].c[N][N].d00 = 
3376                                 crossterm_array[index].c[0][0].d00;
3377
3378     crossterm_setup(&crossterm_array[index].c[0][0]);
3379
3380     ptr->index=index;
3381
3382     index++;
3383     ptr=ptr->next;
3384   }
3385 }
3386 /*      END OF FUNCTION index_crossterms      */
3387
3388 /************************************************************************/
3389 /*                  */
3390 /*      FUNCTION index_vdw        */
3391 /*                  */
3392 /*   INPUTS:                */
3393 /*  tree - The tree that is to be indexed        */
3394 /*  index - index to start with          */
3395 /*                  */
3396 /*  This is a recursive routine that will traverse the binary tree  */
3397 /*   of vdw parameters, assigning an index to each one, and copying     */
3398 /*   the data from the binary tree to the array that will be used from  */
3399 /*   here on.                */
3400 /*                  */
3401 /************************************************************************/
3402
3403 Index Parameters::index_vdw(struct vdw_params *tree, Index index)
3404
3405 {
3406   //  If the tree is empty, do nothing
3407   if (tree==NULL)
3408     return(index);
3409
3410   //  If I have a left subtree, populate it first
3411   if (tree->left != NULL)
3412   {
3413     index=index_vdw(tree->left, index);
3414   }
3415
3416   //  Assign the index and copy the data to the array
3417   tree->index = index;
3418
3419   vdw_array[index].sigma = tree->sigma;
3420   vdw_array[index].epsilon = tree->epsilon;
3421   vdw_array[index].sigma14 = tree->sigma14;
3422   vdw_array[index].epsilon14 = tree->epsilon14;
3423
3424   char *nameloc = atom_type_name(index);
3425   strncpy(nameloc, tree->atomname, MAX_ATOMTYPE_CHARS);
3426   nameloc[MAX_ATOMTYPE_CHARS] = '\0';
3427
3428 //  iout << iWARN << "Parameters: Stored name for type " << index << ": '";
3429 //      iout << iWARN << nameloc << "'" << "\n" << endi;
3430
3431   index++;
3432
3433   //  If I have a right subtree, index it
3434   if (tree->right != NULL)
3435   {
3436     index=index_vdw(tree->right, index);
3437   }
3438
3439   return(index);
3440 }
3441 /*      END OF FUNCTION index_vdw      */
3442
3443 /************************************************************************/
3444 /*                  */
3445 /*      FUNCTION assign_vdw_index      */
3446 /*                  */
3447 /*   INPUTS:                */
3448 /*  atomtype - atom type to find          */
3449 /*  atom_ptr - pointer to the atom structure to find vdw paramters  */
3450 /*       for              */
3451 /*                  */
3452 /*   OUTPUTS:                */
3453 /*  the vdw_index field of the atom structure is populated    */
3454 /*                  */
3455 /*  This function searches the binary tree of vdw parameters so     */
3456 /*   that an index can be assigned to this atom.  If the parameter is   */
3457 /*   is found, then the index is assigned.  If the parameter is not     */
3458 /*   found, then NAMD terminates.          */
3459 /*                  */
3460 /************************************************************************/
3461 #ifdef MEM_OPT_VERSION
3462 void Parameters::assign_vdw_index(char *atomtype, AtomCstInfo *atom_ptr)
3463 #else    
3464 void Parameters::assign_vdw_index(char *atomtype, Atom *atom_ptr)
3465 #endif
3466 {
3467   struct vdw_params *ptr;    //  Current position in trees
3468   int found=0;      //  Flag 1->found match
3469   int comp_code;      //  return code from strcasecmp
3470
3471   /*  Check to make sure the files have all been read    */
3472   if (!AllFilesRead)
3473   {
3474     NAMD_die("Tried to assign vdw index before all parameter files were read");
3475   }
3476
3477   /*  Start at the top            */
3478   ptr=vdwp;
3479
3480   /*  While we haven't found a match, and we haven't reached      */
3481   /*  the bottom of the tree, compare the atom passed in with     */
3482   /*  the current value and decide if we have a match, or if not, */
3483   /*  which way to go            */
3484   while (!found && (ptr!=NULL))
3485   {
3486     comp_code = strcasecmp(atomtype, ptr->atomname);
3487
3488     if (comp_code == 0)
3489     {
3490       /*  Found a match!        */
3491       atom_ptr->vdw_type=ptr->index;
3492       found=1;
3493     }
3494     else if (comp_code < 0)
3495     {
3496       /*  Go to the left        */
3497       ptr=ptr->left;
3498     }
3499     else
3500     {
3501       /*  Go to the right        */
3502       ptr=ptr->right;
3503     }
3504   }
3505
3506   //****** BEGIN CHARMM/XPLOR type changes
3507   if (!found)
3508   {
3509     // since CHARMM allows wildcards "*" in vdw typenames
3510     // we have to look again if necessary, this way, if
3511     // we already had an exact match, this is never executed
3512     size_t windx;                      //  wildcard index
3513
3514     /*  Start again at the top                                */
3515     ptr=vdwp;
3516   
3517      while (!found && (ptr!=NULL))
3518      {
3519   
3520        // get index of wildcard wildcard, get index
3521        windx= strcspn(ptr->atomname,"*"); 
3522        if (windx == strlen(ptr->atomname))
3523        {
3524          // there is no wildcard here
3525          comp_code = strcasecmp(atomtype, ptr->atomname);   
3526        }
3527        else
3528        {
3529          comp_code = strncasecmp(atomtype, ptr->atomname, windx); 
3530        }  
3531
3532        if (comp_code == 0)
3533        {
3534          /*  Found a match!                                */
3535          atom_ptr->vdw_type=ptr->index;
3536          found=1;
3537          char errbuf[100];
3538          sprintf(errbuf,"VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
3539                         atomtype, ptr->atomname);
3540          int i;
3541          for(i=0; i<error_msgs.size(); i++) {
3542            if ( strcmp(errbuf,error_msgs[i]) == 0 ) break;
3543          }
3544          if ( i == error_msgs.size() ) {
3545            char *newbuf = new char[strlen(errbuf)+1];
3546            strcpy(newbuf,errbuf);
3547            error_msgs.add(newbuf);
3548            iout << iWARN << newbuf << "\n" << endi;
3549          }
3550        }
3551        else if (comp_code < 0)
3552        {
3553           /*  Go to the left                                */
3554                 ptr=ptr->left;
3555        }
3556        else
3557        {
3558          /*  Go to the right                                */
3559                 ptr=ptr->right;
3560        }
3561      
3562      }
3563                 
3564   }
3565   //****** END CHARMM/XPLOR type changes
3566
3567   /*  Make sure we found it          */
3568   if (!found)
3569   {
3570     char err_msg[100];
3571
3572     sprintf(err_msg, "DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s",
3573        atomtype);
3574     NAMD_die(err_msg);
3575   }
3576
3577   return;
3578 }
3579 /*      END OF FUNCTION assign_vdw_index    */
3580
3581 /************************************************************************
3582  * FUNCTION get_table_pair_params
3583  *
3584  * Inputs:
3585  * atom1 - atom type for atom 1
3586  * atom2 - atom type for atom 2
3587  * K - an integer value for the table type to populate
3588  *
3589  * Outputs:
3590  * If a match is found, K is populated and 1 is returned. Otherwise,
3591  * 0 is returned.
3592  *
3593  * This function finds the proper type index for tabulated nonbonded 
3594  * interactions between two atoms. If no such interactions are found,
3595  * the atoms are assumed to interact through standard VDW potentials.
3596  * 
3597  ************************************************************************/
3598
3599 int Parameters::get_table_pair_params(Index ind1, Index ind2, int *K) {
3600   IndexedTablePair *ptr;
3601   Index temp;
3602   int found=FALSE;
3603
3604   ptr=tab_pair_tree;
3605   //
3606   //  We need the smaller type in ind1, so if it isn't already that 
3607   //  way, switch them        */
3608   if (ind1 > ind2)
3609   {
3610     temp = ind1;
3611     ind1 = ind2;
3612     ind2 = temp;
3613   }
3614
3615   /*  While we haven't found a match and we're not at the end  */
3616   /*  of the tree, compare the bond passed in with the tree  */
3617   while (!found && (ptr!=NULL))
3618   {
3619 //    printf("Comparing %i with %i and %i with %i\n", ind1, ptr->ind1, ind2, ptr->ind2);
3620     if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
3621     {
3622        found = TRUE;
3623     }
3624     else if ( (ind1 < ptr->ind1) || 
3625         ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
3626     {
3627       /*  Go left          */
3628       ptr=ptr->left;
3629     }
3630     else
3631     {
3632       /*  Go right          */
3633       ptr=ptr->right;
3634     }
3635   }
3636
3637   /*  If we found a match, assign the values      */
3638   if (found)
3639   {
3640     *K = ptr->K;
3641     return(TRUE);
3642   }
3643   else
3644   {
3645     return(FALSE);
3646   }
3647 }
3648 /*      END OF FUNCTION get_table_pair_params    */
3649
3650 /************************************************************************/
3651 /*                  */
3652 /*      FUNCTION get_vdw_pair_params      */
3653 /*                  */
3654 /*   INPUTS:                */
3655 /*  atom1 - atom type for atom 1          */
3656 /*  atom2 - atom type for atom 2          */
3657 /*  A - A value to populate            */
3658 /*  B - B value to populate            */
3659 /*  A14 - A 1-4 value to populate          */
3660 /*  B14 - B 1-4 value to populate          */
3661 /*                  */
3662 /*   OUTPUTS:                */
3663 /*  If a match is found, A, B, A14, and B14 are all populated and a */
3664 /*   1 is returned.  Otherwise, a 0 is returned.      */
3665 /*                    */
3666 /*  This function finds a set of vdw_pair paramters.  It is given   */
3667 /*   the two types of atoms involved.  This is the only paramter for    */
3668 /*   which a match is NOT guaranteed.  There will only be a match if    */
3669 /*   there are specific van der waals parameters for the two atom types */
3670 /*   involved.                */
3671 /*                  */
3672 /************************************************************************/
3673
3674 int Parameters::get_vdw_pair_params(Index ind1, Index ind2, Real *A, 
3675         Real *B, Real *A14, Real *B14)
3676
3677 {
3678   IndexedVdwPair *ptr;    //  Current location in tree
3679   Index temp;      //  Temporary value for swithcing
3680           // values
3681   int found=FALSE;    //  Flag 1-> found a match
3682
3683   ptr=vdw_pair_tree;
3684
3685   //  We need the smaller type in ind1, so if it isn't already that 
3686   //  way, switch them        */
3687   if (ind1 > ind2)
3688   {
3689     temp = ind1;
3690     ind1 = ind2;
3691     ind2 = temp;
3692   }
3693
3694   /*  While we haven't found a match and we're not at the end  */
3695   /*  of the tree, compare the bond passed in with the tree  */
3696   while (!found && (ptr!=NULL))
3697   {
3698     if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
3699     {
3700        found = TRUE;
3701     }
3702     else if ( (ind1 < ptr->ind1) || 
3703         ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
3704     {
3705       /*  Go left          */
3706       ptr=ptr->left;
3707     }
3708     else
3709     {
3710       /*  Go right          */
3711       ptr=ptr->right;
3712     }
3713   }
3714
3715   /*  If we found a match, assign the values      */
3716   if (found)
3717   {
3718     *A = ptr->A;
3719     *B = ptr->B;
3720     *A14 = ptr->A14;
3721     *B14 = ptr->B14;
3722
3723     return(TRUE);
3724   }
3725   else
3726   {
3727     return(FALSE);
3728   }
3729 }
3730 /*      END OF FUNCTION get_vdw_pair_params    */
3731
3732
3733 /************************************************************************/
3734 /*                  */
3735 /*        FUNCTION assign_bond_index    */
3736 /*                  */
3737 /*   INPUTS:                */
3738 /*  atom1 - atom type for atom 1          */
3739 /*  atom2 - atom type for atom 2          */
3740 /*  bond_ptr - pointer to bond structure to populate    */
3741 /*                  */
3742 /*   OUTPUTS:                */
3743 /*  the structure pointed to by bond_ptr is populated    */
3744 /*                  */
3745 /*  This function finds a bond in the binary tree of bond values    */
3746 /*   and assigns its index.  If the bond is found, than the bond_type   */
3747 /*   field of the bond structure is populated.  If the parameter is     */
3748 /*   not found, NAMD will terminate.          */
3749 /*                  */
3750 /************************************************************************/
3751
3752 void Parameters::assign_bond_index(char *atom1, char *atom2, Bond *bond_ptr)
3753
3754 {
3755   struct bond_params *ptr;  //  Current location in tree
3756   int found=0;      //  Flag 1-> found a match
3757   int cmp_code;      //  return code from strcasecmp
3758   char tmp_name[15];    //  Temporary atom name
3759
3760   /*  Check to make sure the files have all been read    */
3761   if (!AllFilesRead)
3762   {
3763     NAMD_die("Tried to assign bond index before all parameter files were read");
3764   }
3765
3766   /*  We need atom1 < atom2, so if that's not the way they        */
3767   /*  were passed, flip them          */
3768   if (strcasecmp(atom1, atom2) > 0)
3769   {
3770     strcpy(tmp_name, atom1);
3771     strcpy(atom1, atom2);
3772     strcpy(atom2, tmp_name);
3773   }
3774
3775   /*  Start at the top            */
3776   ptr=bondp;
3777
3778   /*  While we haven't found a match and we're not at the end  */
3779   /*  of the tree, compare the bond passed in with the tree  */
3780   while (!found && (ptr!=NULL))
3781   {
3782     cmp_code=strcasecmp(atom1, ptr->atom1name);
3783
3784     if (cmp_code == 0)
3785     {
3786       cmp_code=strcasecmp(atom2, ptr->atom2name);
3787     }
3788
3789     if (cmp_code == 0)
3790     {
3791       /*  Found a match        */
3792       found=1;
3793       bond_ptr->bond_type = ptr->index;
3794     }
3795     else if (cmp_code < 0)
3796     {
3797       /*  Go left          */
3798       ptr=ptr->left;
3799     }
3800     else
3801     {
3802       /*  Go right          */
3803       ptr=ptr->right;
3804     }
3805   }
3806
3807   /*  Check to see if we found anything        */
3808   if (!found)
3809   {
3810     if ((strcmp(atom1, "DRUD")==0 || strcmp(atom2, "DRUD")==0)
3811         && (strcmp(atom1, "X")!=0 && strcmp(atom2, "X")!=0)) {
3812       /* try a wildcard DRUD X match for this Drude bond */
3813       char a1[8] = "DRUD", a2[8] = "X";
3814       return assign_bond_index(a1, a2, bond_ptr);  /* recursive call */
3815     }
3816     else {
3817       char err_msg[128];
3818
3819       sprintf(err_msg, "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
3820          atom1, atom2, bond_ptr->atom1+1, bond_ptr->atom2+1);
3821       NAMD_die(err_msg);
3822     }
3823   }
3824
3825   return;
3826 }
3827 /*      END OF FUNCTION assign_bond_index    */
3828
3829 /************************************************************************/
3830 /*                  */
3831 /*      FUNCTION assign_angle_index      */
3832 /*                  */
3833 /*   INPUTS:                */
3834 /*  atom1 - atom type for atom 1          */
3835 /*  atom2 - atom type for atom 2          */
3836 /*  atom3 - atom type for atom 3          */
3837 /*  angle_ptr - pointer to angle structure to populate    */
3838 /*                  */
3839 /*   OUTPUTS:                */
3840 /*  the structure pointed to by angle_ptr is populated    */
3841 /*                  */
3842 /*  This function assigns an angle index to a specific angle.  */
3843 /*   It searches the binary tree of angle parameters for the appropriate*/
3844 /*   values.  If they are found, the index is assigned.  If they are    */
3845 /*   not found, then NAMD will terminate.        */
3846 /*                  */
3847 /************************************************************************/
3848
3849 void Parameters::assign_angle_index(char *atom1, char *atom2, char*atom3,
3850           Angle *angle_ptr, int notFoundIndex)
3851
3852 {
3853   struct angle_params *ptr;  //  Current position in tree
3854   int comp_val;      //  value from strcasecmp
3855   int found=0;      //  flag 1->found a match
3856   char tmp_name[15];    //  Temporary atom name
3857
3858   /*  Check to make sure the files have all been read    */
3859   if (!AllFilesRead)
3860   {
3861     NAMD_die("Tried to assign angle index before all parameter files were read");
3862   }
3863
3864   /*  We need atom1 < atom3.  If that was not what we were   */
3865   /*  passed, switch them            */
3866   if (strcasecmp(atom1, atom3) > 0)
3867   {
3868     strcpy(tmp_name, atom1);
3869     strcpy(atom1, atom3);
3870     strcpy(atom3, tmp_name);
3871   }
3872
3873   /*  Start at the top            */
3874   ptr=anglep;
3875
3876   /*  While we don't have a match and we haven't reached the  */
3877   /*  bottom of the tree, compare values        */
3878   while (!found && (ptr != NULL))
3879   {
3880     comp_val = strcasecmp(atom1, ptr->atom1name);
3881
3882     if (comp_val == 0)
3883     {
3884       /*  Atom 1 matches, so compare atom 2    */
3885       comp_val = strcasecmp(atom2, ptr->atom2name);
3886       
3887       if (comp_val == 0)
3888       {
3889         /*  Atoms 1&2 match, try atom 3    */
3890         comp_val = strcasecmp(atom3, ptr->atom3name);
3891       }
3892     }
3893
3894     if (comp_val == 0)
3895     {
3896       /*  Found a match        */
3897       found = 1;
3898       angle_ptr->angle_type = ptr->index;
3899     }
3900     else if (comp_val < 0)
3901     {
3902       /*  Go left          */
3903       ptr=ptr->left;
3904     }
3905     else
3906     {
3907       /*  Go right          */
3908       ptr=ptr->right;
3909     }
3910   }
3911
3912   /*  Make sure we found a match          */
3913   if (!found)
3914   {
3915     char err_msg[128];
3916
3917     sprintf(err_msg, "UNABLE TO FIND ANGLE PARAMETERS FOR %s %s %s (ATOMS %i %i %i)",
3918        atom1, atom2, atom3, angle_ptr->atom1+1, angle_ptr->atom2+1, angle_ptr->atom3+1);
3919
3920     if ( notFoundIndex ) {
3921       angle_ptr->angle_type = notFoundIndex;
3922       iout << iWARN << err_msg << "\n" << endi;
3923       return;
3924     } else NAMD_die(err_msg);
3925   }
3926
3927   return;
3928 }
3929 /*      END OF FUNCTION assign_angle_index    */
3930
3931 /************************************************************************/
3932 /*                  */
3933 /*      FUNCTION assign_dihedral_index      */
3934 /*                  */
3935 /*   INPUTS:                */
3936 /*  atom1 - atom type for atom 1          */
3937 /*  atom2 - atom type for atom 2          */
3938 /*  atom3 - atom type for atom 3          */
3939 /*  atom4 - atom type for atom 4          */
3940 /*  dihedral_ptr - pointer to dihedral structure to populate  */
3941 /*  multiplicity - Multiplicity to assign to this bond    */
3942 /*                  */
3943 /*   OUTPUTS:                */
3944 /*  the structure pointed to by dihedral_ptr is populated    */
3945 /*                  */
3946 /*  This function searchs the linked list of dihedral parameters for*/
3947 /*   a given bond.  If a match is found, the dihedral type is assigned. */
3948 /*   If no match is found, NAMD terminates        */
3949 /*                  */
3950 /************************************************************************/
3951
3952 void Parameters::assign_dihedral_index(char *atom1, char *atom2, char *atom3,
3953         char *atom4, Dihedral *dihedral_ptr,
3954         int multiplicity, int notFoundIndex)
3955
3956 {
3957   struct dihedral_params *ptr;  //  Current position in list
3958   int found=0;      //  Flag 1->found a match
3959
3960   /*  Start at the begining of the list        */
3961   ptr=dihedralp;
3962
3963   /*  While we haven't found a match and we haven't reached       */
3964   /*  the end of the list, keep looking        */
3965   while (!found && (ptr!=NULL))
3966   {
3967     /*  Do a linear search through the linked list of   */
3968     /*  dihedral parameters.  Since the list is arranged    */
3969     /*  with wildcard paramters at the end of the list, we  */
3970     /*  can simply do a linear search and be guaranteed that*/
3971     /*  we will find exact matches before wildcard matches. */
3972     /*  Also, we must check for an exact match, and a match */
3973     /*  in reverse, since they are really the same          */
3974     /*  physically.            */
3975     if ( ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom1)==0) ) && 
3976          ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom2)==0) ) &&
3977          ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom3)==0) ) &&
3978          ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom4)==0) ) ) 
3979     {
3980       /*  Found an exact match      */
3981       found=1;
3982     }
3983     else if ( ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom1)==0) ) &&
3984               ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom2)==0) ) &&
3985               ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom3)==0) ) &&
3986               ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom4)==0) ) )
3987     {
3988       /*  Found a reverse match      */
3989       found=1;
3990     }
3991     else
3992     {
3993       /*  Didn't find a match, go to the next node  */
3994       ptr=ptr->next;
3995     }
3996   }
3997
3998   /*  Make sure we found a match          */
3999   if (!found)
4000   {
4001     char err_msg[128];
4002
4003     sprintf(err_msg, "UNABLE TO FIND DIHEDRAL PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
4004        atom1, atom2, atom3, atom4, dihedral_ptr->atom1+1, dihedral_ptr->atom2+1,
4005        dihedral_ptr->atom3+1, dihedral_ptr->atom4+1);
4006     
4007     if ( notFoundIndex ) {
4008       dihedral_ptr->dihedral_type = notFoundIndex;
4009       iout << iWARN << err_msg << "\n" << endi;
4010       return;
4011     } else NAMD_die(err_msg);
4012   }
4013
4014  if (paramType == paraXplor) {
4015   //  Check to make sure the number of multiples specified in the psf
4016   //  file doesn't exceed the number of parameters in the parameter
4017   //  files
4018   if (multiplicity > maxDihedralMults[ptr->index])
4019   {
4020     char err_msg[257];
4021
4022     sprintf(err_msg, "Multiplicity of Paramters for diehedral bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxDihedralMults[ptr->index]);
4023     NAMD_die(err_msg);
4024   }
4025
4026   //  If the multiplicity from the current bond is larger than that
4027   //  seen in the past, increase the multiplicity for this bond
4028   if (multiplicity > dihedral_array[ptr->index].multiplicity)
4029   {
4030     dihedral_array[ptr->index].multiplicity = multiplicity;
4031   }
4032  }
4033
4034   dihedral_ptr->dihedral_type = ptr->index;
4035
4036   return;
4037 }
4038 /*      END OF FUNCTION assign_dihedral_index    */
4039
4040 /************************************************************************/
4041 /*                  */
4042 /*      FUNCTION assign_improper_index      */
4043 /*                  */
4044 /*   INPUTS:                */
4045 /*  atom1 - atom type for atom 1          */
4046 /*  atom2 - atom type for atom 2          */
4047 /*  atom3 - atom type for atom 3          */
4048 /*  atom4 - atom type for atom 4          */
4049 /*  improper_ptr - pointer to improper structure to populate  */
4050 /*   multiplicity - Multiplicity to assign to this bond    */
4051 /*                  */
4052 /*   OUTPUTS:                */
4053 /*  the structure pointed to by improper_ptr is populated    */
4054 /*                  */
4055 /*  This function searchs the linked list of improper parameters for*/
4056 /*   a given bond.  If a match is found, the improper_type is assigned. */
4057 /*   If no match is found, NAMD will terminate.        */
4058 /*                  */
4059 /************************************************************************/
4060
4061 void Parameters::assign_improper_index(char *atom1, char *atom2, char *atom3,
4062         char *atom4, Improper *improper_ptr,
4063         int multiplicity)
4064
4065 {
4066   struct improper_params *ptr;  //  Current position in list
4067   int found=0;      //  Flag 1->found a match
4068
4069   /*  Start at the head of the list        */
4070   ptr=improperp;
4071
4072   /*  While we haven't fuond a match and haven't reached the end  */
4073   /*  of the list, keep looking          */
4074   while (!found && (ptr!=NULL))
4075   {
4076     /*  Do a linear search through the linked list of   */
4077     /*  improper parameters.  Since the list is arranged    */
4078     /*  with wildcard paramters at the end of the list, we  */
4079     /*  can simply do a linear search and be guaranteed that*/
4080     /*  we will find exact matches before wildcard matches. */
4081     /*  Also, we must check for an exact match, and a match */
4082     /*  in reverse, since they are really the same          */
4083     /*  physically.            */
4084     if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
4085            (strcasecmp(ptr->atom1name, "X")==0) ) &&
4086        ( (strcasecmp(ptr->atom2name, atom2)==0) || 
4087            (strcasecmp(ptr->atom2name, "X")==0) ) &&
4088        ( (strcasecmp(ptr->atom3name, atom3)==0) || 
4089            (strcasecmp(ptr->atom3name, "X")==0) ) &&
4090        ( (strcasecmp(ptr->atom4name, atom4)==0) || 
4091            (strcasecmp(ptr->atom4name, "X")==0) ) )
4092     {
4093       /*  Found an exact match      */
4094       found=1;
4095     }
4096     else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
4097            (strcasecmp(ptr->atom4name, "X")==0) ) &&
4098        ( (strcasecmp(ptr->atom3name, atom2)==0) || 
4099            (strcasecmp(ptr->atom3name, "X")==0) ) &&
4100        ( (strcasecmp(ptr->atom2name, atom3)==0) || 
4101            (strcasecmp(ptr->atom2name, "X")==0) ) &&
4102        ( (strcasecmp(ptr->atom1name, atom4)==0) || 
4103            (strcasecmp(ptr->atom1name, "X")==0) ) )
4104     {
4105       /*  Found a reverse match      */
4106       found=1;
4107     }
4108     else
4109     {
4110       /*  Didn't find a match, go to the next node  */
4111       ptr=ptr->next;
4112     }
4113   }
4114
4115   /*  Make sure we found a match          */
4116   if (!found)
4117   {
4118     char err_msg[128];
4119
4120     sprintf(err_msg, "UNABLE TO FIND IMPROPER PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
4121        atom1, atom2, atom3, atom4, improper_ptr->atom1+1, improper_ptr->atom2+1,
4122        improper_ptr->atom3+1, improper_ptr->atom4+1);
4123     
4124     NAMD_die(err_msg);
4125   }
4126
4127  if (paramType == paraXplor) {
4128   //  Check to make sure the number of multiples specified in the psf
4129   //  file doesn't exceed the number of parameters in the parameter
4130   //  files
4131   if (multiplicity > maxImproperMults[ptr->index])
4132   {
4133     char err_msg[257];
4134
4135     sprintf(err_msg, "Multiplicity of Paramters for improper bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxImproperMults[ptr->index]);
4136     NAMD_die(err_msg);
4137   }
4138
4139   //  If the multiplicity from the current bond is larger than that
4140   //  seen in the past, increase the multiplicity for this bond
4141   if (multiplicity > improper_array[ptr->index].multiplicity)
4142   {
4143     improper_array[ptr->index].multiplicity = multiplicity;
4144   }
4145  }
4146
4147   /*  Assign the constants          */
4148   improper_ptr->improper_type = ptr->index;
4149
4150   return;
4151 }
4152 /*      END OF FUNCTION assign_improper_index    */
4153
4154 /************************************************************************/
4155 /*                  */
4156 /*      FUNCTION assign_crossterm_index      */
4157 /*                  */
4158 /************************************************************************/
4159
4160 void Parameters::assign_crossterm_index(char *atom1, char *atom2, char *atom3,
4161         char *atom4, char *atom5, char *atom6, char *atom7,
4162         char *atom8, Crossterm *crossterm_ptr)
4163 {
4164   struct crossterm_params *ptr;  //  Current position in list
4165   int found=0;      //  Flag 1->found a match
4166
4167   /*  Start at the head of the list        */
4168   ptr=crosstermp;
4169
4170   /*  While we haven't fuond a match and haven't reached the end  */
4171   /*  of the list, keep looking          */
4172   while (!found && (ptr!=NULL))
4173   {
4174     /*  Do a linear search through the linked list of   */
4175     /*  crossterm parameters.  Since the list is arranged    */
4176     /*  with wildcard paramters at the end of the list, we  */
4177     /*  can simply do a linear search and be guaranteed that*/
4178     /*  we will find exact matches before wildcard matches. */
4179     /*  Also, we must check for an exact match, and a match */
4180     /*  in reverse, since they are really the same          */
4181     /*  physically.            */
4182     if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
4183            (strcasecmp(ptr->atom1name, "X")==0) ) &&
4184        ( (strcasecmp(ptr->atom2name, atom2)==0) || 
4185            (strcasecmp(ptr->atom2name, "X")==0) ) &&
4186        ( (strcasecmp(ptr->atom3name, atom3)==0) || 
4187            (strcasecmp(ptr->atom3name, "X")==0) ) &&
4188        ( (strcasecmp(ptr->atom4name, atom4)==0) || 
4189            (strcasecmp(ptr->atom4name, "X")==0) ) )
4190     {
4191       /*  Found an exact match      */
4192       found=1;
4193     }
4194     else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
4195            (strcasecmp(ptr->atom4name, "X")==0) ) &&
4196        ( (strcasecmp(ptr->atom3name, atom2)==0) || 
4197            (strcasecmp(ptr->atom3name, "X")==0) ) &&
4198        ( (strcasecmp(ptr->atom2name, atom3)==0) || 
4199            (strcasecmp(ptr->atom2name, "X")==0) ) &&
4200        ( (strcasecmp(ptr->atom1name, atom4)==0) || 
4201            (strcasecmp(ptr->atom1name, "X")==0) ) )
4202     {
4203       /*  Found a reverse match      */
4204       found=1;
4205     }
4206     if ( ! found ) {
4207       /*  Didn't find a match, go to the next node  */
4208       ptr=ptr->next;
4209       continue;
4210     }
4211     found = 0;
4212     if ( ( (strcasecmp(ptr->atom5name, atom5)==0) || 
4213            (strcasecmp(ptr->atom5name, "X")==0) ) &&
4214        ( (strcasecmp(ptr->atom6name, atom6)==0) || 
4215            (strcasecmp(ptr->atom6name, "X")==0) ) &&
4216        ( (strcasecmp(ptr->atom7name, atom7)==0) || 
4217            (strcasecmp(ptr->atom7name, "X")==0) ) &&
4218        ( (strcasecmp(ptr->atom8name, atom8)==0) || 
4219            (strcasecmp(ptr->atom8name, "X")==0) ) )
4220     {
4221       /*  Found an exact match      */
4222       found=1;
4223     }
4224     else if ( ( (strcasecmp(ptr->atom8name, atom5)==0) || 
4225            (strcasecmp(ptr->atom8name, "X")==0) ) &&
4226        ( (strcasecmp(ptr->atom7name, atom6)==0) || 
4227            (strcasecmp(ptr->atom7name, "X")==0) ) &&
4228        ( (strcasecmp(ptr->atom6name, atom7)==0) || 
4229            (strcasecmp(ptr->atom6name, "X")==0) ) &&
4230        ( (strcasecmp(ptr->atom5name, atom8)==0) || 
4231            (strcasecmp(ptr->atom5name, "X")==0) ) )
4232     {
4233       /*  Found a reverse match      */
4234       found=1;
4235     }
4236     if ( ! found ) {