NAMD: Parameters.C Source File

00001 
00007 /*
00008    The class Parameters is used to hold all of the parameters read
00009    in from the parameter files.  The class provides a routine to read in
00010    parameter files (as many parameter files as desired can be read in) and
00011    a series of routines that allow the parameters that have been read in
00012    to be queried.
00013 */
00014 
00015 #include <stdio.h>
00016 #include <stdlib.h>
00017 #include <string.h>
00018 #include <vector>
00019 #ifndef WIN32
00020 #include <strings.h>
00021 #endif
00022 #include "InfoStream.h"
00023 #include <charm++.h>
00024 #include "Parameters.h"
00025 #include "Communicate.h"
00026 #include "ConfigList.h"
00027 //****** BEGIN CHARMM/XPLOR type changes
00028 #include "SimParameters.h"
00029 //****** END CHARMM/XPLOR type changes
00030 
00031 #define MIN_DEBUG_LEVEL 3
00032 //#define DEBUGM
00033 #include "Debug.h"
00034 
00035 #define INDEX(ncols,i,j)  ((i)*ncols + (j))
00036 
00037 #define ENABLETABLES
00038 
00039 static char** table_types;
00040 
00041 //  struct bond_params is used to form a binary tree of bond parameters.
00042 //  The two atom names are used to determine the order of the nodes in the
00043 //  tree.  atom1name should ALWAYS be lexically before atom2name
00044 
00045 struct bond_params
00046 {
00047   char atom1name[11];
00048   char atom2name[11];
00049   Real forceconstant;
00050   Real distance;
00051   Index index;
00052   struct bond_params *left;
00053   struct bond_params *right;
00054 };
00055 
00056 //  struct angle_params is used to form a binary tree of bond parameters.
00057 //  The three atom names are used to determine the order of the nodes in
00058 //  the tree.  atom1name should ALWAYS be lexically before atom3name
00059 
00060 struct angle_params
00061 {
00062   char atom1name[11];
00063   char atom2name[11];
00064   char atom3name[11];
00065   Real forceconstant;
00066   int normal;
00067   Real angle;
00068   Real k_ub;
00069   Real r_ub;
00070   Index index;
00071   struct angle_params *left;
00072   struct angle_params *right;
00073 };
00074 
00075 //  struct dihedral_params is used to form a linked list of the dihedral
00076 //  parameters.  The linked list is arranged in such a way that any
00077 //  bonds with wildcards are at the end of the list so that a linear
00078 //  search can be done but we will still find exact matches before
00079 //  wildcard matches
00080 
00081 struct dihedral_params
00082 {
00083   char atom1name[11];
00084   char atom2name[11];
00085   char atom3name[11];
00086   char atom4name[11];
00087   char atom1wild;
00088   char atom2wild;
00089   char atom3wild;
00090   char atom4wild;
00091   int multiplicity;
00092   FourBodyConsts values[MAX_MULTIPLICITY];
00093   Index index;
00094   struct dihedral_params *next;
00095   dihedral_params() { memset(this, 0, sizeof(dihedral_params)); }
00096 };
00097 
00098 //  struct improper_params is used to form a linked list of the improper
00099 //  parameters.  The linked list is arranged in such a way that any
00100 //  bonds with wildcards are at the end of the list so that a linear
00101 //  search can be done but we will still find exact matches before
00102 //  wildcard matches
00103 
00104 struct improper_params
00105 {
00106   char atom1name[11];
00107   char atom2name[11];
00108   char atom3name[11];
00109   char atom4name[11];
00110   int multiplicity;
00111   FourBodyConsts values[MAX_MULTIPLICITY];
00112   Index index;
00113   struct improper_params *next;
00114 };
00115 
00116 struct crossterm_params
00117 {
00118   crossterm_params(int dim) : dimension(dim) {
00119     values = new double[dimension*dimension];
00120   }
00121   ~crossterm_params() {
00122     delete [] values;
00123   }
00124   char atom1name[11];
00125   char atom2name[11];
00126   char atom3name[11];
00127   char atom4name[11];
00128   char atom5name[11];
00129   char atom6name[11];
00130   char atom7name[11];
00131   char atom8name[11];
00132   int dimension;  // usually 24
00133   double *values;  // dimension * dimension data
00134   Index index;
00135   struct crossterm_params *next;
00136 };
00137 
00138 //  struct vdw_params is used to form a binary serach tree of the
00139 //  vdw paramters for a single atom.
00140 
00141 struct vdw_params
00142 {
00143   char atomname[11];
00144   Real sigma;
00145   Real epsilon;
00146   Real sigma14;
00147   Real epsilon14;
00148   Index index;
00149   struct vdw_params *left;
00150   struct vdw_params *right;
00151 };
00152 
00153 //  struct vdw_pair_params is used to form a linked list of the
00154 //  vdw parameters for a pair of atoms
00155 
00156 struct vdw_pair_params
00157 {
00158   char atom1name[11];
00159   char atom2name[11];
00160   Real A;
00161   Real B;
00162   Real A14;
00163   Real B14;
00164   struct vdw_pair_params *next;
00165 };
00166 
00167 struct table_pair_params
00168 {
00169   char atom1name[11];
00170   char atom2name[11];
00171   int K;
00172   struct table_pair_params *next;
00173 };
00174 
00175 Parameters::Parameters() {
00176   initialize();
00177 }
00178 
00179 void Parameters::initialize() {
00180 
00181   paramType = -1;
00182 
00183   /*  Set all the pointers to NULL        */
00184   atomTypeNames=NULL;
00185   bondp=NULL;
00186   anglep=NULL;
00187   improperp=NULL;
00188   dihedralp=NULL;
00189   crosstermp=NULL;
00190   vdwp=NULL;
00191   vdw_pairp=NULL;
00192   table_pairp=NULL;
00193   bond_array=NULL;
00194   angle_array=NULL;
00195   dihedral_array=NULL;
00196   improper_array=NULL;
00197   crossterm_array=NULL;
00198   vdw_array=NULL;
00199   vdw_pair_tree=NULL;
00200   tab_pair_tree=NULL;
00201   maxDihedralMults=NULL;
00202   maxImproperMults=NULL;
00203 
00204   /*  Set all the counts to 0          */
00205   NumBondParams=0;
00206   NumAngleParams=0;
00207   NumDihedralParams=0;
00208   NumImproperParams=0;
00209   NumCrosstermParams=0;
00210   NumVdwParams=0;
00211   NumVdwPairParams=0;
00212   NumTablePairParams=0;
00213   NumCosAngles=0;
00214 }
00215 
00216 /************************************************************************/
00217 /*                  */
00218 /*      FUNCTION Parameters        */
00219 /*                  */
00220 /*  This is the constructor for the class.  It simlpy sets the      */
00221 /*  pointers to the list and trees to NULL and the count of all the     */
00222 /*  parameters to 0.              */
00223 /*  The type (format) of the input parameters (Xplor,Charmm) is set here. */
00224 /*                  */
00225 /************************************************************************/
00226 
00227 Parameters::Parameters(SimParameters *simParams, StringList *f)
00228 {
00229   initialize();
00230 
00231   //****** BEGIN CHARMM/XPLOR type changes
00233   if (simParams->paraTypeXplorOn)
00234   {
00235     paramType = paraXplor;
00236   }
00237   else if (simParams->paraTypeCharmmOn)
00238   {
00239     paramType = paraCharmm;
00240   }
00241   //****** END CHARMM/XPLOR type changes
00242   //Test for cos-based angles
00243   if (simParams->cosAngles) {
00244     cosAngles = true;
00245   } else {
00246     cosAngles = false;
00247   }
00248 
00249   /* Set up AllFilesRead flag to FALSE.  Once all of the files    */
00250   /* have been read in, then this will be set to true and the     */
00251   /* arrays of parameters will be set up        */
00252   AllFilesRead = FALSE;
00253 
00254   numenerentries=0;
00255   table_ener = NULL;
00256   if (simParams->tabulatedEnergies) {
00257 
00258           fprintf(stdout,"Working on tables\n");
00259           read_ener_table(simParams);
00260   }
00261 
00262   if (NULL != f) 
00263   {
00264     do
00265     {
00266       //****** BEGIN CHARMM/XPLOR type changes
00267       if (paramType == paraXplor)
00268       {
00269         read_parameter_file(f->data);
00270       }
00271       else if (paramType == paraCharmm)
00272       {
00273         read_charmm_parameter_file(f->data);
00274       }
00275       //****** END CHARMM/XPLOR type changes
00276       f = f->next;
00277     } while ( f != NULL );
00278 
00279     done_reading_files();
00280   }
00281 
00282 }
00283 /*      END OF FUNCTION Parameters      */
00284 
00285 /************************************************************************/
00286 /*                  */
00287 /*      FUNCTION ~Parameters        */
00288 /*                  */
00289 /*  This is the destructor for this class.  It basically just       */
00290 /*  frees all of the memory allocated for the parameters.    */
00291 /*                  */
00292 /************************************************************************/
00293 
00294 Parameters::~Parameters()
00295 
00296 {
00297         if (atomTypeNames)
00298           delete [] atomTypeNames;
00299 
00300   if (bondp != NULL)
00301     free_bond_tree(bondp);
00302 
00303   if (anglep != NULL)
00304     free_angle_tree(anglep);
00305 
00306   if (dihedralp != NULL)
00307     free_dihedral_list(dihedralp);
00308 
00309   if (improperp != NULL)
00310     free_improper_list(improperp);
00311 
00312   if (crosstermp != NULL)
00313     free_crossterm_list(crosstermp);
00314 
00315   if (vdwp != NULL)
00316     free_vdw_tree(vdwp);
00317 
00318   if (vdw_pairp != NULL)
00319     free_vdw_pair_list();
00320 
00321   if (bond_array != NULL)
00322     delete [] bond_array;
00323 
00324   if (angle_array != NULL)
00325     delete [] angle_array;
00326 
00327   if (dihedral_array != NULL)
00328     delete [] dihedral_array;
00329 
00330   if (improper_array != NULL)
00331     delete [] improper_array;
00332 
00333   if (crossterm_array != NULL)
00334     delete [] crossterm_array;
00335 
00336   if (vdw_array != NULL)
00337     delete [] vdw_array;
00338   
00339   if (tab_pair_tree != NULL)
00340     free_table_pair_tree(tab_pair_tree);
00341 
00342   if (vdw_pair_tree != NULL)
00343     free_vdw_pair_tree(vdw_pair_tree);
00344 
00345   if (maxDihedralMults != NULL)
00346     delete [] maxDihedralMults;
00347 
00348   if (maxImproperMults != NULL)
00349     delete [] maxImproperMults;
00350 
00351   for( int i = 0; i < error_msgs.size(); ++i ) {
00352     delete [] error_msgs[i];
00353   }
00354   error_msgs.resize(0);
00355 }
00356 /*      END OF FUNCTION ~Parameters      */
00357 
00358 /************************************************************************/
00359 /*                  */
00360 /*      FUNCTION read_paramter_file      */
00361 /*                  */
00362 /*   INPUTS:                */
00363 /*  fname - name of the parameter file to read      */
00364 /*                  */
00365 /*  This function reads in a parameter file and adds the parameters */
00366 /*   from this file to the current group of parameters.  The basic      */
00367 /*   logic of the routine is to read in a line from the file, looks at  */
00368 /*   the first word of the line to determine what kind of parameter we  */
00369 /*   have, and then call the appropriate routine to add the parameter   */
00370 /*   to the parameters that we have.          */
00371 /*                  */
00372 /************************************************************************/
00373 
00374 void Parameters::read_parameter_file(char *fname)
00375 
00376 {
00377   char buffer[512];  //  Buffer to store each line of the file
00378   char first_word[512];  //  First word of the current line
00379   FILE *pfile;    //  File descriptor for the parameter file
00380 
00381   /*  Check to make sure that we haven't previously been told     */
00382   /*  that all the files were read        */
00383   if (AllFilesRead)
00384   {
00385     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00386   }
00387 
00388   /*  Try and open the file          */
00389   if ( (pfile = Fopen(fname, "r")) == NULL)
00390   {
00391     char err_msg[256];
00392 
00393     sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
00394     NAMD_die(err_msg);
00395   }
00396 
00397   /*  Keep reading in lines until we hit the EOF      */
00398   while (NAMD_read_line(pfile, buffer) != -1)
00399   {
00400     /*  Get the first word of the line      */
00401     NAMD_find_first_word(buffer, first_word);
00402 
00403     /*  First, screen out things that we ignore, such as    */
00404     /*  blank lines, lines that start with '!', lines that  */
00405     /*  start with "REMARK", lines that start with set",    */
00406     /*  and most of the HBOND parameters which include      */
00407     /*  AEXP, REXP, HAEX, AAEX, but not the HBOND statement */
00408     /*  which is parsed.                                    */
00409     if ((buffer[0] != '!') && 
00410         !NAMD_blank_string(buffer) &&
00411         (strncasecmp(first_word, "REMARK", 6) != 0) &&
00412         (strcasecmp(first_word, "set")!=0) &&
00413         (strncasecmp(first_word, "AEXP", 4) != 0) &&
00414         (strncasecmp(first_word, "REXP", 4) != 0) &&
00415         (strncasecmp(first_word, "HAEX", 4) != 0) &&
00416         (strncasecmp(first_word, "AAEX", 4) != 0) &&
00417         (strncasecmp(first_word, "NBOND", 5) != 0) &&
00418         (strncasecmp(first_word, "CUTNB", 5) != 0) &&
00419         (strncasecmp(first_word, "END", 3) != 0) &&
00420         (strncasecmp(first_word, "CTONN", 5) != 0) &&
00421         (strncasecmp(first_word, "EPS", 3) != 0) &&
00422         (strncasecmp(first_word, "VSWI", 4) != 0) &&
00423         (strncasecmp(first_word, "NBXM", 4) != 0) &&
00424         (strncasecmp(first_word, "INHI", 4) != 0) )
00425     {
00426       /*  Now, call the appropriate function based    */
00427       /*  on the type of parameter we have    */
00428       if (strncasecmp(first_word, "bond", 4)==0)
00429       {
00430         add_bond_param(buffer);
00431         NumBondParams++;
00432       }
00433       else if (strncasecmp(first_word, "angl", 4)==0)
00434       {
00435         add_angle_param(buffer);
00436         NumAngleParams++;
00437       }
00438       else if (strncasecmp(first_word, "dihe", 4)==0)
00439       {
00440         add_dihedral_param(buffer, pfile);
00441         NumDihedralParams++;
00442       }
00443       else if (strncasecmp(first_word, "impr", 4)==0)
00444       {
00445         add_improper_param(buffer, pfile);
00446         NumImproperParams++;
00447       }
00448       else if (strncasecmp(first_word, "nonb", 4)==0)
00449       {
00450         add_vdw_param(buffer);
00451         NumVdwParams++; 
00452       }
00453       else if (strncasecmp(first_word, "nbfi", 4)==0)
00454       {
00455         add_vdw_pair_param(buffer);
00456         NumVdwPairParams++; 
00457       }
00458       else if (strncasecmp(first_word, "nbta", 4)==0)
00459       {
00460 
00461         if (table_ener == NULL) {
00462           continue;
00463         }
00464 
00465         add_table_pair_param(buffer);
00466         NumTablePairParams++; 
00467       }
00468       else if (strncasecmp(first_word, "hbon", 4)==0)
00469       {
00470         add_hb_pair_param(buffer);
00471       }
00472       else
00473       {
00474         /*  This is an unknown paramter.        */
00475         /*  This is BAD        */
00476         char err_msg[512];
00477 
00478         sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
00479            fname, buffer);
00480         NAMD_die(err_msg);
00481       }
00482     }
00483   }
00484 
00485   /*  Close the file            */
00486   Fclose(pfile);
00487 
00488   return;
00489 }
00490 /*      END OF FUNCTION read_paramter_file    */
00491 
00492 //****** BEGIN CHARMM/XPLOR type changes
00493 /************************************************************************/
00494 /*                                                                        */
00495 /*                        FUNCTION read_charmm_paramter_file                */
00496 /*                                                                        */
00497 /*   INPUTS:                                                                */
00498 /*        fname - name of the parameter file to read                        */
00499 /*                                                                        */
00500 /*        This function reads in a CAHRMM parameter file and adds the     */ 
00501 /*   parameters from this file to the current group of parameters.      */
00502 /*   The basic logic of the routine is to first find out what type of   */
00503 /*   parameter we have in the file. Then look at each line in turn      */
00504 /*   and call the appropriate routine to add the parameters until we hit*/
00505 /*   a new type of parameter or EOF.                                    */
00506 /*                                                                        */
00507 /************************************************************************/
00508 
00509 void Parameters::read_charmm_parameter_file(char *fname)
00510 
00511 {
00512   int  par_type=0;         //  What type of parameter are we currently
00513                            //  dealing with? (vide infra)
00514   int  skipline;           //  skip this line?
00515   int  skipall = 0;        //  skip rest of file;
00516   char buffer[512];           //  Buffer to store each line of the file
00517   char first_word[512];           //  First word of the current line
00518   FILE *pfile;                   //  File descriptor for the parameter file
00519 
00520   /*  Check to make sure that we haven't previously been told     */
00521   /*  that all the files were read                                */
00522   if (AllFilesRead)
00523   {
00524     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00525   }
00526 
00527   /*  Try and open the file                                        */
00528   if ( (pfile = fopen(fname, "r")) == NULL)
00529   {
00530     char err_msg[256];
00531 
00532     sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
00533     NAMD_die(err_msg);
00534   }
00535 
00536   /*  Keep reading in lines until we hit the EOF                        */
00537   while (NAMD_read_line(pfile, buffer) != -1)
00538   {
00539     /*  Get the first word of the line                        */
00540     NAMD_find_first_word(buffer, first_word);
00541     skipline=0;
00542 
00543     /*  First, screen out things that we ignore.                   */   
00544     /*  blank lines, lines that start with '!' or '*', lines that  */
00545     /*  start with "END".                                          */
00546     if (!NAMD_blank_string(buffer) &&
00547         (strncmp(first_word, "!", 1) != 0) &&
00548          (strncmp(first_word, "*", 1) != 0) &&
00549         (strncasecmp(first_word, "END", 3) != 0))
00550     {
00551       if ( skipall ) {
00552         iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
00553         break;
00554       }
00555       /*  Now, determine the apropriate parameter type.   */
00556       if (strncasecmp(first_word, "bond", 4)==0)
00557       {
00558         par_type=1; skipline=1;
00559       }
00560       else if (strncasecmp(first_word, "angl", 4)==0)
00561       {
00562         par_type=2; skipline=1;
00563       }
00564       else if (strncasecmp(first_word, "thet", 4)==0)
00565       {
00566         par_type=2; skipline=1;
00567       }
00568       else if (strncasecmp(first_word, "dihe", 4)==0)
00569       {
00570         par_type=3; skipline=1;
00571       }
00572       else if (strncasecmp(first_word, "phi", 3)==0)
00573       {
00574         par_type=3; skipline=1;
00575       }
00576       else if (strncasecmp(first_word, "impr", 4)==0)
00577       {
00578         par_type=4; skipline=1;
00579       }
00580       else if (strncasecmp(first_word, "imph", 4)==0)
00581       {
00582         par_type=4; skipline=1;
00583       }
00584       else if (strncasecmp(first_word, "nbon", 4)==0)
00585       {
00586         par_type=5; skipline=1;
00587       }
00588       else if (strncasecmp(first_word, "nonb", 4)==0)
00589       {
00590         par_type=5; skipline=1;
00591       }
00592       else if (strncasecmp(first_word, "nbfi", 4)==0)
00593       {
00594         par_type=6; skipline=1;
00595       }
00596       else if (strncasecmp(first_word, "hbon", 4)==0)
00597       {
00598         par_type=7; skipline=1;
00599       }
00600       else if (strncasecmp(first_word, "cmap", 4)==0)
00601       {
00602         par_type=8; skipline=1;
00603       }
00604       else if (strncasecmp(first_word, "nbta", 4)==0)
00605       {
00606         par_type=9; skipline=1;
00607       }
00608       else if (strncasecmp(first_word, "read", 4)==0)
00609       {
00610         skip_stream_read(buffer, pfile);  skipline=1;
00611       }
00612       else if (strncasecmp(first_word, "return", 4)==0)
00613       {
00614         skipall=1;  skipline=1;
00615       }
00616       else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
00617                (strncasecmp(first_word, "grou", 4) == 0) ||
00618                (strncasecmp(first_word, "cdie", 4) == 0) ||
00619                (strncasecmp(first_word, "shif", 4) == 0) ||
00620                (strncasecmp(first_word, "vgro", 4) == 0) ||
00621                (strncasecmp(first_word, "vdis", 4) == 0) ||
00622                (strncasecmp(first_word, "vswi", 4) == 0) ||
00623                (strncasecmp(first_word, "cutn", 4) == 0) ||
00624                (strncasecmp(first_word, "ctof", 4) == 0) ||
00625                (strncasecmp(first_word, "cton", 4) == 0) ||
00626                (strncasecmp(first_word, "eps", 3) == 0) ||
00627                (strncasecmp(first_word, "e14f", 4) == 0) ||
00628                (strncasecmp(first_word, "wmin", 4) == 0) ||
00629                (strncasecmp(first_word, "aexp", 4) == 0) ||
00630                (strncasecmp(first_word, "rexp", 4) == 0) ||
00631                (strncasecmp(first_word, "haex", 4) == 0) ||
00632                (strncasecmp(first_word, "aaex", 4) == 0) ||
00633                (strncasecmp(first_word, "noac", 4) == 0) ||
00634                (strncasecmp(first_word, "hbno", 4) == 0) ||
00635                (strncasecmp(first_word, "cuth", 4) == 0) ||
00636                (strncasecmp(first_word, "ctof", 4) == 0) ||
00637                (strncasecmp(first_word, "cton", 4) == 0) ||
00638                (strncasecmp(first_word, "cuth", 4) == 0) ||
00639                (strncasecmp(first_word, "ctof", 4) == 0) ||
00640                (strncasecmp(first_word, "cton", 4) == 0) ) 
00641       {
00642         if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
00643         {
00644           char err_msg[512];
00645 
00646           sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00647           NAMD_die(err_msg);
00648         }
00649         else 
00650         {
00651           skipline = 1;
00652         }
00653       }        
00654       else if (par_type == 0)
00655       {
00656         /*  This is an unknown paramter.        */
00657         /*  This is BAD                                */
00658         char err_msg[512];
00659 
00660         sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00661         NAMD_die(err_msg);
00662       }
00663     }
00664     else
00665     {
00666       skipline=1;
00667     }
00668 
00669     if ( (par_type != 0) && (!skipline) )
00670     {
00671       /*  Now, call the appropriate function based    */
00672       /*  on the type of parameter we have                */
00673       /*  I know, this should really be a switch ...  */
00674       if (par_type == 1)
00675       {
00676         add_bond_param(buffer);
00677         NumBondParams++;
00678       }
00679       else if (par_type == 2)
00680       {
00681         add_angle_param(buffer);
00682         NumAngleParams++;
00683       }
00684       else if (par_type == 3)
00685       {
00686         add_dihedral_param(buffer, pfile);
00687         NumDihedralParams++;
00688       }
00689       else if (par_type == 4)
00690       {
00691         add_improper_param(buffer, pfile);
00692         NumImproperParams++;
00693       }
00694       else if (par_type == 5)
00695       {
00696         add_vdw_param(buffer);
00697         NumVdwParams++;
00698       }
00699       else if (par_type == 6)
00700       {
00701         add_vdw_pair_param(buffer);
00702         NumVdwPairParams++; 
00703       }
00704       else if (par_type == 7)
00705       {
00706         add_hb_pair_param(buffer);                  
00707       }
00708       else if (par_type == 8)
00709       {
00710         add_crossterm_param(buffer, pfile);                  
00711         NumCrosstermParams++;
00712       }
00713       else if (par_type == 9)
00714       {
00715 
00716         if (table_ener == NULL) {
00717           continue;
00718         }
00719 
00720         add_table_pair_param(buffer);                  
00721         NumTablePairParams++;
00722       }
00723       else
00724       {
00725         /*  This really should not occour!      */
00726         /*  This is an internal error.          */
00727         /*  This is VERY BAD                        */
00728         char err_msg[512];
00729 
00730         sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00731         NAMD_die(err_msg);
00732       }
00733     }
00734   }
00735 
00736   /*  Close the file                                                */
00737   fclose(pfile);
00738 
00739   return;
00740 }
00741 /*                        END OF FUNCTION read_charmm_paramter_file                */
00742 //****** END CHARMM/XPLOR type changes
00743 
00744 
00745 void Parameters::skip_stream_read(char *buf, FILE *fd) {
00746 
00747   char buffer[513];
00748   char first_word[513];
00749   char s1[128];
00750   char s2[128];
00751   int rval = sscanf(buf, "%s %s", s1, s2);
00752   if (rval != 2) {
00753         char err_msg[512];
00754         sprintf(err_msg, "BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
00755         NAMD_die(err_msg);
00756   }
00757   if ( ! strncasecmp(s2,"PARA",4) ) return;  // read parameters
00758   
00759   iout << iINFO << "SKIPPING " << s2 << " SECTION IN STREAM FILE\n" << endi;
00760 
00761   while (NAMD_read_line(fd, buffer) != -1)
00762   {
00763     // read until we find "END"
00764     NAMD_find_first_word(buffer, first_word);
00765     if (!NAMD_blank_string(buffer) &&
00766         (strncmp(first_word, "!", 1) != 0) &&
00767          (strncmp(first_word, "*", 1) != 0) &&
00768          (strncasecmp(first_word, "END", 3) == 0) ) return;
00769   }
00770 
00771 }
00772 
00773 
00774 /************************************************************************/
00775 /*                  */
00776 /*      FUNCTION add_bond_param        */
00777 /*                  */
00778 /*   INPUTS:                */
00779 /*  buf - Line from parameter file containing bond parameters  */
00780 /*                  */
00781 /*  This function adds a new bond paramter to the binary tree of    */
00782 /*   angle paramters that we have.  If a duplicate is found, a warning  */
00783 /*   message is printed and the new parameters are used.    */
00784 /*                  */
00785 /************************************************************************/
00786 
00787 void Parameters::add_bond_param(char *buf)
00788 
00789 {
00790   char atom1name[11];    //  Atom type for atom 1
00791   char atom2name[11];    //  Atom type for atom 2
00792   Real forceconstant;    //  Force constant for bond
00793   Real distance;      //  Rest distance for bond
00794   int read_count;      //  Count from sscanf
00795   struct bond_params *new_node;  //  New node in tree
00796 
00797   //****** BEGIN CHARMM/XPLOR type changes
00798   /*  Use sscanf to parse up the input line      */
00799   if (paramType == paraXplor)
00800   {
00801     /* read XPLOR format */
00802     read_count=sscanf(buf, "%*s %s %s %f %f\n", atom1name, atom2name, 
00803        &forceconstant, &distance);
00804   }
00805   else if (paramType == paraCharmm)
00806   {
00807     /* read CHARMM format */
00808     read_count=sscanf(buf, "%s %s %f %f\n", atom1name, atom2name, 
00809        &forceconstant, &distance);
00810   }
00811   //****** END CHARMM/XPLOR type changes
00812 
00813   /*  Check to make sure we found everything we expeceted    */
00814   if (read_count != 4)
00815   {
00816     char err_msg[512];
00817 
00818     if (paramType == paraXplor)
00819     {
00820       sprintf(err_msg, "BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
00821     }
00822     else if (paramType == paraCharmm)
00823     {
00824       sprintf(err_msg, "BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
00825     }
00826     NAMD_die(err_msg);
00827   }
00828 
00829   /*  Allocate a new node            */
00830   new_node = new bond_params;
00831 
00832   if (new_node == NULL)
00833   {
00834     NAMD_die("memory allocation failed in Parameters::add_bond_param\n");
00835   }
00836 
00837   /*  Order the atoms so that the atom that comes alphabetically  */
00838   /*  first is atom 1.  Since the bond is symmetric, it doesn't   */
00839   /*  matter physically which atom is first.  And this allows the */
00840   /*  search of the binary tree to be done in a logical manner    */
00841   if (strcasecmp(atom1name, atom2name) < 0)
00842   {
00843     strcpy(new_node->atom1name, atom1name);
00844     strcpy(new_node->atom2name, atom2name);
00845   }
00846   else
00847   {
00848     strcpy(new_node->atom2name, atom1name);
00849     strcpy(new_node->atom1name, atom2name);
00850   }
00851 
00852   /*  Assign force constant and distance        */
00853   new_node->forceconstant = forceconstant;
00854   new_node->distance = distance;
00855 
00856   /*  Set pointers to null          */
00857   new_node->left = NULL;
00858   new_node->right = NULL;
00859 
00860   /*  Make call to recursive call to actually add the node to the */
00861   /*  tree              */
00862   bondp=add_to_bond_tree(new_node, bondp);
00863 
00864   return;
00865 }
00866 /*      END OF FUNCTION add_bond_param      */
00867 
00868 /************************************************************************/
00869 /*                  */
00870 /*      FUNCTION add_to_bond_tree      */
00871 /*                  */
00872 /*   INPUTS:                */
00873 /*  new_node - Node to add to the tree        */
00874 /*  tree - tree to add the node to          */
00875 /*                  */
00876 /*   OUTPUTS:                */
00877 /*  ths function returns a pointer to the new tree with the node    */
00878 /*   added to it.  Most of the time it will be the same pointer as was  */
00879 /*   passed in, but not if the current tree is empty.      */
00880 /*                  */
00881 /*  this is a receursive function that adds a node to the binary    */
00882 /*   tree used to store bond parameters.        */
00883 /*                  */
00884 /************************************************************************/
00885 
00886 struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
00887              struct bond_params *tree)
00888 
00889 {
00890   int compare_code;  //  Results from strcasecmp
00891 
00892   /*  If the tree is currently empty, then the new tree consists  */
00893   /*  only of the new node          */
00894   if (tree == NULL)
00895     return(new_node);
00896 
00897   /*  Compare the atom1 name from the new node and the head of    */
00898   /*  the tree              */
00899   compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
00900 
00901   /*  Check to see if they are the same        */
00902   if (compare_code == 0)
00903   {
00904     /*  The atom 1 names are the same, compare atom 2  */
00905     compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
00906 
00907     /*  If atom 1 AND atom 2 are the same, we have a duplicate */
00908     if (compare_code == 0)
00909     {
00910       /*  We have a duplicate.  So print out a warning*/
00911       /*  message.  Then assign the new values to the */
00912       /*  tree and free the new_node      */
00913       //****** BEGIN CHARMM/XPLOR type changes
00914       /* we do not care about identical replacement */
00915       if ((tree->forceconstant != new_node->forceconstant) || 
00916           (tree->distance != new_node->distance))
00917       {
00918         iout << "\n" << iWARN << "DUPLICATE BOND ENTRY FOR "
00919           << new_node->atom1name << "-"
00920           << new_node->atom2name
00921           << "\nPREVIOUS VALUES  k=" << tree->forceconstant
00922           << "  x0=" << tree->distance
00923           << "\n   USING VALUES  k=" << new_node->forceconstant
00924           << "  x0=" << new_node->distance
00925           << "\n" << endi;
00926 
00927         tree->forceconstant=new_node->forceconstant;
00928         tree->distance=new_node->distance;
00929       }
00930       //****** END CHARMM/XPLOR type changes
00931 
00932       delete new_node;
00933 
00934       return(tree);
00935     }
00936   }
00937 
00938   /*  We don't have a duplicate, so if the new value is less      */
00939   /*  than the head of the tree, add it to the left child,   */
00940   /*  otherwise add it to the right child        */
00941   if (compare_code < 0)
00942   {
00943     tree->left = add_to_bond_tree(new_node, tree->left);
00944   }
00945   else
00946   {
00947     tree->right = add_to_bond_tree(new_node, tree->right);
00948   }
00949 
00950   return(tree);
00951 }
00952 /*    END OF FUNCTION add_to_bond_tree      */
00953 
00954 /************************************************************************/
00955 /*                  */
00956 /*      FUNCTION add_angle_param      */
00957 /*                  */
00958 /*   INPUTS:                */
00959 /*  buf - line from paramter file with angle parameters    */
00960 /*                  */
00961 /*  this function adds an angle parameter.  It parses up the input  */
00962 /*   line and then adds it to the binary tree used to store the angle   */
00963 /*   parameters.              */
00964 /*                  */
00965 /************************************************************************/
00966 
00967 void Parameters::add_angle_param(char *buf)
00968 
00969 {
00970   char atom1name[11];    // Type for atom 1
00971   char atom2name[11];    // Type for atom 2
00972   char atom3name[11];    // Type for atom 3
00973   char norm[4]="xxx";
00974   Real forceconstant;    // Force constant
00975   Real angle;      // Theta 0
00976   Real k_ub;      // Urey-Bradley force constant
00977   Real r_ub;      // Urey-Bradley distance
00978   int read_count;      // count from sscanf
00979   struct angle_params *new_node;  // new node in tree
00980 
00981   //****** BEGIN CHARMM/XPLOR type changes
00982   /*  parse up the input line with sscanf        */
00983   if (paramType == paraXplor)
00984   {
00985     /* read XPLOR format */
00986     read_count=sscanf(buf, "%*s %s %s %s %f %f UB %f %f\n", 
00987        atom1name, atom2name, atom3name, &forceconstant, &angle,
00988        &k_ub, &r_ub);
00989   }
00990   else if ((paramType == paraCharmm) && cosAngles) {
00991     read_count=sscanf(buf, "%s %s %s %f %f %3s %f %f\n", 
00992        atom1name, atom2name, atom3name, &forceconstant, &angle, norm,
00993        &k_ub, &r_ub);
00994 //    printf("%s\n", buf);
00995 //    printf("Data: %s %s %s %f %f %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle, norm, k_ub, r_ub);
00996   }  
00997   else if (paramType == paraCharmm)
00998   {
00999     /* read CHARMM format */
01000     read_count=sscanf(buf, "%s %s %s %f %f %f %f\n", 
01001        atom1name, atom2name, atom3name, &forceconstant, &angle,
01002        &k_ub, &r_ub);
01003 //    printf("%s\n", buf);
01004 //    printf("Data: %s %s %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle);
01005   }
01006   //****** END CHARMM/XPLOR type changes
01007 
01008   /*  Check to make sure we got what we expected      */
01009   if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
01010   {
01011     char err_msg[512];
01012 
01013     if (paramType == paraXplor)
01014     {
01015       sprintf(err_msg, "BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01016     }
01017     else if (paramType == paraCharmm)
01018     {
01019       sprintf(err_msg, "BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
01020     }
01021     NAMD_die(err_msg);
01022   }
01023 
01024   /*  Allocate the new node          */
01025   new_node = new angle_params;
01026 
01027   if (new_node == NULL)
01028   {
01029     NAMD_die("memory allocation failed in Parameters::add_angle_param");
01030   }
01031 
01032   /*  As with the bond, we want the atom type is comes first  */
01033   /*  alphbetically first between atom 1 and atom 3 to be in      */
01034   /*  atom 1 so that we can search the tree reliably.    */
01035   if (strcasecmp(atom1name, atom3name) < 0)
01036   {
01037     strcpy(new_node->atom1name, atom1name);
01038     strcpy(new_node->atom2name, atom2name);
01039     strcpy(new_node->atom3name, atom3name);
01040   }
01041   else
01042   {
01043     strcpy(new_node->atom3name, atom1name);
01044     strcpy(new_node->atom2name, atom2name);
01045     strcpy(new_node->atom1name, atom3name);
01046   }
01047 
01048   /*  Assign the constants and pointer values      */
01049   new_node->forceconstant = forceconstant;
01050   new_node->angle = angle;
01051 
01052   if (cosAngles) {
01053     if (strcmp("cos",norm)==0) {
01054 //      iout << "Info: Using cos mode for angle " << buf << endl;
01055       NumCosAngles++;
01056       new_node->normal = 0;
01057     } else {
01058 //      iout << "Info: Using x^2 mode for angle " << buf << endl;
01059       new_node->normal = 1;
01060     }
01061   } else {
01062     new_node->normal = 1;
01063   }
01064 
01065   if (read_count == 7)
01066   {
01067     //  Urey-Bradley constants
01068     if (new_node->normal == 0) {
01069       NAMD_die("ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
01070     }
01071     new_node->k_ub = k_ub;
01072     new_node->r_ub = r_ub;
01073   }
01074   else
01075   {
01076     new_node->k_ub = 0.0;
01077     new_node->r_ub = 0.0;
01078   }
01079 
01080   new_node->left = NULL;
01081   new_node->right = NULL;
01082 
01083   /*  Insert it into the tree          */
01084   anglep = add_to_angle_tree(new_node, anglep);
01085 
01086   return;
01087 }
01088 /*      END OF FUNCTION add_angle_param      */
01089 
01090 /************************************************************************/
01091 /*                  */
01092 /*      FUNCTION add_to_angle_tree      */
01093 /*                  */
01094 /*   INPUTS:                */
01095 /*  new_node - new node to add to the angle tree      */
01096 /*  tree - tree to add the node to          */
01097 /*                  */
01098 /*   OUTPUTS:                */
01099 /*  the function returns a pointer to the new tree with the node    */
01100 /*   added.  Most of the time, this will be the same as the value passed*/
01101 /*   in, but not in the case where the tree is empty.      */
01102 /*                  */
01103 /*  this is a recursive function that adds an angle parameter  */
01104 /*   to the binary tree storing the angle parameters.  If a duplicate   */
01105 /*   is found, a warning message is printed, the current values in the  */
01106 /*   tree are replaced with the new values, and the new node is free'd  */
01107 /*                  */
01108 /************************************************************************/
01109 
01110 struct angle_params *Parameters::add_to_angle_tree(struct angle_params *new_node,
01111              struct angle_params *tree)
01112 
01113 {
01114   int compare_code;  //  Return code from strcasecmp
01115 
01116   /*  If the tree is empty, then the new_node is the tree    */
01117   if (tree == NULL)
01118     return(new_node);
01119 
01120   /*  Compare atom 1 from the new node and the head of the tree   */
01121   compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
01122 
01123   if (compare_code == 0)
01124   {
01125     /*  Atom 1 is the same, compare atom 2      */
01126     compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
01127 
01128     if (compare_code == 0)
01129     {
01130       /*  Atoms 1 & 2 are the same, compare atom 3  */
01131       compare_code = strcasecmp(new_node->atom3name, 
01132             tree->atom3name);
01133 
01134       if (compare_code == 0)
01135       {
01136         /*  All three atoms were the same, this */
01137         /*  is a duplicate.  Print a warning    */
01138         /*  message, replace the current values,*/
01139         /*  and free the new node    */
01140         //****** BEGIN CHARMM/XPLOR type changes
01141         /* we do not care about identical replacement */
01142         if ((tree->forceconstant != new_node->forceconstant) ||
01143             (tree->angle != new_node->angle) ||
01144             (tree->k_ub != new_node->k_ub) ||
01145             (tree->r_ub != new_node->r_ub) || (tree->normal != new_node->normal))
01146         {
01147           iout << "\n" << iWARN << "DUPLICATE ANGLE ENTRY FOR "
01148             << new_node->atom1name << "-"
01149             << new_node->atom2name << "-"
01150             << new_node->atom3name
01151             << "\nPREVIOUS VALUES  k="
01152             << tree->forceconstant << "  theta0="
01153             << tree->angle << " k_ub="
01154             << tree->k_ub << " r_ub="
01155             << tree->r_ub
01156             << "\n   USING VALUES  k="
01157             << new_node->forceconstant << "  theta0="
01158             << new_node->angle << " k_ub="
01159             << new_node->k_ub << " r_ub=" << new_node->r_ub 
01160             << "\n" << endi;
01161 
01162           tree->forceconstant=new_node->forceconstant;
01163           tree->angle=new_node->angle;
01164           tree->k_ub=new_node->k_ub;
01165           tree->r_ub=new_node->r_ub;
01166           tree->normal=new_node->normal;
01167         }
01168         //****** END CHARMM/XPLOR type changes
01169 
01170         delete new_node;
01171 
01172         return(tree);
01173       }
01174     }
01175   }
01176 
01177   /*  Didn't find a duplicate, so if the new_node is smaller  */
01178   /*  than the current head, add it to the left child.  Otherwise */
01179   /*  add it to the right child.          */
01180   if (compare_code < 0)
01181   {
01182     tree->left = add_to_angle_tree(new_node, tree->left);
01183   }
01184   else
01185   {
01186     tree->right = add_to_angle_tree(new_node, tree->right);
01187   }
01188 
01189   return(tree);
01190 }
01191 /*      END OF FUNCTION add_to_angle_tree    */
01192 
01193 /************************************************************************/
01194 /*                  */
01195 /*      FUNCTION add_dihedral_param      */
01196 /*                  */
01197 /*   INPUTS:                */
01198 /*  buf - line from paramter file with dihedral parameters    */
01199 /*                  */
01200 /*  this function adds an dihedral parameter.  It parses up the     */
01201 /*   input line and then adds it to the binary tree used to store the   */
01202 /*   dihedral parameters.            */
01203 /*                  */
01204 /************************************************************************/
01205 
01206 void Parameters::add_dihedral_param(char *buf, FILE *fd)
01207 
01208 {
01209   char atom1name[11];       //  Type of atom 1
01210   char atom2name[11];       //  Type of atom 2
01211   char atom3name[11];       //  Type of atom 3
01212   char atom4name[11];       //  Type of atom 4
01213   Real forceconstant;       //  Force constant
01214   int periodicity;       //  Periodicity
01215   Real phase_shift;       //  Phase shift
01216   int read_count;         //  Count from sscanf
01217   struct dihedral_params *new_node;  //  New node
01218   int multiplicity;       //  Multiplicity for bonds
01219   int i;           //  Loop counter
01220   char buffer[513];       //  Buffer for new line
01221   int ret_code;         //  Return code
01222 
01223   //****** BEGIN CHARMM/XPLOR type changes
01224   /*  Parse up the input line using sscanf      */
01225   if (paramType == paraXplor)
01226   {
01227     /* read XPLOR format */
01228     read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n", 
01229        atom1name, atom2name, atom3name, atom4name, &multiplicity,
01230        &forceconstant, &periodicity, &phase_shift);
01231   }
01232   else if (paramType == paraCharmm)
01233   {
01234     /* read CHARMM format */
01235     read_count=sscanf(buf, "%s %s %s %s %f %d %f\n", 
01236        atom1name, atom2name, atom3name, atom4name,
01237        &forceconstant, &periodicity, &phase_shift);
01238     multiplicity=1; 
01239   }
01240 
01241   if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
01242   {
01243     char err_msg[512];
01244 
01245     sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01246     NAMD_die(err_msg);
01247   }
01248   else if ( (read_count != 7) && (paramType == paraCharmm) )
01249   {
01250     char err_msg[512];
01251 
01252     sprintf(err_msg, "BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
01253     NAMD_die(err_msg);
01254   }
01255 
01256   if ( (read_count == 4) && (paramType == paraXplor) )
01257   //****** END CHARMM/XPLOR type changes
01258   {
01259     read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n", 
01260           &forceconstant, &periodicity, &phase_shift);
01261 
01262     /*  Check to make sure we got what we expected    */
01263     if (read_count != 3)
01264     {
01265       char err_msg[512];
01266 
01267       sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01268       NAMD_die(err_msg);
01269     }
01270 
01271     multiplicity = 1;
01272   }
01273 
01274   if (multiplicity > MAX_MULTIPLICITY)
01275   {
01276     char err_msg[181];
01277 
01278     sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
01279        multiplicity, MAX_MULTIPLICITY);
01280     NAMD_die(err_msg);
01281   }
01282 
01283   /*  Allocate new node            */
01284   new_node = new dihedral_params;
01285 
01286   if (new_node == NULL)
01287   {
01288     NAMD_die("memory allocation failed in Parameters::add_dihedral_param\n");
01289   }
01290 
01291   /*  Assign all of the values for this node.  Notice that since  */
01292   /*  the dihedrals and impropers are implemented with a linked   */
01293   /*  list rather than a binary tree, we don't really care about  */
01294   /*  the order of the atoms any more        */
01295   strcpy(new_node->atom1name, atom1name);
01296   strcpy(new_node->atom2name, atom2name);
01297   strcpy(new_node->atom3name, atom3name);
01298   strcpy(new_node->atom4name, atom4name);
01299   new_node->atom1wild = ! strcasecmp(atom1name, "X");
01300   new_node->atom2wild = ! strcasecmp(atom2name, "X");
01301   new_node->atom3wild = ! strcasecmp(atom3name, "X");
01302   new_node->atom4wild = ! strcasecmp(atom4name, "X");
01303   new_node->multiplicity = multiplicity;
01304   if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01305   new_node->values[0].k = forceconstant;
01306   new_node->values[0].n = periodicity;
01307   new_node->values[0].delta = phase_shift;
01308 
01309   new_node->next = NULL;
01310 
01311   //  If the multiplicity is greater than 1, then read in other parameters
01312   if (multiplicity > 1)
01313   {
01314     for (i=1; i<multiplicity; i++)
01315     {
01316       ret_code = NAMD_read_line(fd, buffer);
01317 
01318       //  Get rid of comments at the end of a line
01319       if (ret_code == 0)
01320       {
01321         NAMD_remove_comment(buffer);
01322       }
01323 
01324       //  Keep reading lines until we get one that isn't blank
01325       while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
01326       {
01327         ret_code = NAMD_read_line(fd, buffer);
01328       }
01329 
01330       if (ret_code != 0)
01331       {
01332         NAMD_die("EOF encoutner in middle of multiple dihedral");
01333       }
01334 
01335       read_count=sscanf(buffer, "%f %d %f\n", 
01336             &forceconstant, &periodicity, &phase_shift);
01337 
01338       if (read_count != 3)
01339       {
01340         char err_msg[512];
01341 
01342         sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
01343         NAMD_die(err_msg);
01344       }
01345 
01346       if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01347       new_node->values[i].k = forceconstant;
01348       new_node->values[i].n = periodicity;
01349       new_node->values[i].delta = phase_shift;
01350     }
01351   }
01352 
01353   //****** BEGIN CHARMM/XPLOR type changes
01354   /*  Add this node to the list          */
01355   if (paramType == paraXplor)
01356   {
01357     add_to_dihedral_list(new_node); // XPLOR
01358   }
01359   else if (paramType == paraCharmm)
01360   {
01361     add_to_charmm_dihedral_list(new_node); // CHARMM
01362   }
01363  //****** END CHARMM/XPLOR type changes
01364 
01365   return;
01366 }
01367 /*      END OF FUNCTION add_dihedral_param    */
01368 
01369 /************************************************************************/
01370 /*                  */
01371 /*      FUNCTION add_to_dihedral_list      */
01372 /*                  */
01373 /*   INPUTS:                */
01374 /*  new_node - node that is to be added to dihedral_list    */
01375 /*                  */
01376 /*  this function adds a new dihedral parameter to the linked list  */
01377 /*   of dihedral parameters.  First, it checks for duplicates.  If a    */
01378 /*   duplicate is found, a warning message is printed, the old values   */
01379 /*   are replaced with the new values, and the new node is freed.  If   */
01380 /*   Otherwise, the node is added to the list.  This list is arranged   */
01381 /*   so that bods with wildcards are placed at the tail of the list.    */
01382 /*   This will guarantee that if we just do a linear search, we will    */
01383 /*   always find an exact match before a wildcard match.    */
01384 /*                  */
01385 /************************************************************************/
01386 
01387 void Parameters::add_to_dihedral_list(
01388         struct dihedral_params *new_node)
01389 
01390 {
01391   static struct dihedral_params *ptr;   //  position within list
01392   static struct dihedral_params *tail;  //  Pointer to the end of 
01393                 //  the list so we can add
01394                 //  entries to the end of the
01395                 //  list in constant time
01396   int i;              //  Loop counter
01397 
01398   /*  If the list is currently empty, then the new node is the list*/
01399   if (dihedralp == NULL)
01400   {
01401     dihedralp=new_node;
01402     tail=new_node;
01403 
01404     return;
01405   }
01406 
01407   /*  The list isn't empty, so check for a duplicate    */
01408   ptr=dihedralp;
01409 
01410   while (ptr != NULL)
01411   {
01412     if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01413            (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01414            (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01415            (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01416          ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01417            (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01418            (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01419            (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
01420     {
01421       /*  Found a duplicate        */
01422       //****** BEGIN CHARMM/XPLOR type changes
01423       /* we do not care about identical replacement */
01424       int echoWarn=0;  // echo warning messages ?
01425 
01426       if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
01427       
01428       if (!echoWarn)
01429       {
01430         for (i=0; i<ptr->multiplicity; i++)
01431         {
01432           if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
01433           if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
01434           if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
01435         }
01436       }
01437 
01438       if (echoWarn)
01439       {
01440         iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
01441           << ptr->atom1name << "-"
01442           << ptr->atom2name << "-"
01443           << ptr->atom3name << "-"
01444           << ptr->atom4name
01445           << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
01446         
01447         for (i=0; i<ptr->multiplicity; i++)
01448         {
01449           iout     << "  k=" << ptr->values[i].k
01450                    << "  n=" << ptr->values[i].n
01451                    << "  delta=" << ptr->values[i].delta;
01452         }
01453 
01454         iout << "\nUSING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
01455 
01456         for (i=0; i<new_node->multiplicity; i++)
01457         {
01458           iout <<     "  k=" << new_node->values[i].k
01459                    << "  n=" << new_node->values[i].n
01460                    << "  delta=" << new_node->values[i].delta;
01461         }
01462 
01463         iout << endi;
01464 
01465         ptr->multiplicity = new_node->multiplicity;
01466 
01467         for (i=0; i<new_node->multiplicity; i++)
01468         {
01469           ptr->values[i].k = new_node->values[i].k;
01470           ptr->values[i].n = new_node->values[i].n;
01471           ptr->values[i].delta = new_node->values[i].delta;
01472         }
01473 
01474       }
01475       //****** END CHARMM/XPLOR type changes
01476 
01477       delete new_node;
01478 
01479       return;
01480     }
01481 
01482     ptr=ptr->next;
01483   }
01484 
01485   /*  Check to see if we have any wildcards.  Since specific  */
01486   /*  entries are to take precedence, we'll put anything without  */
01487   /*  wildcards at the begining of the list and anything with     */
01488   /*  wildcards at the end of the list.  Then, we can just do a   */
01489   /*  linear search for a bond and be guaranteed to have specific */
01490   /*  entries take precendence over over wildcards          */
01491   if ( new_node->atom1wild ||
01492        new_node->atom2wild ||
01493        new_node->atom3wild ||
01494        new_node->atom4wild )
01495   {
01496     /*  add to the end of the list        */
01497     tail->next=new_node;
01498     tail=new_node;
01499 
01500     return;
01501   }
01502   else
01503   {
01504     /*  add to the head of the list        */
01505     new_node->next=dihedralp;
01506     dihedralp=new_node;
01507 
01508     return;
01509   }
01510 
01511 }
01512 /*    END OF FUNCTION add_to_dihedral_list      */
01513 
01514 //****** BEGIN CHARMM/XPLOR type changes
01515 /************************************************************************/
01516 /*                                                                        */
01517 /*                        FUNCTION add_to_charmm_dihedral_list                */
01518 /*                                                                        */
01519 /*   INPUTS:                                                                */
01520 /*        new_node - node that is to be added to dihedral_list                */
01521 /*                                                                        */
01522 /*        this function adds a new dihedral parameter to the linked list  */
01523 /*   of dihedral parameters in CHARMM format.                           */
01524 /*   First, it checks for duplicates.  If a duplicate is found, a       */
01525 /*   warning message is printed. If the periodicity is the same as of   */
01526 /*   a previous dihedral the old values are replaced with the new       */
01527 /*   values, otherwise, the dihedral is added and the multiplicity is   */
01528 /*   increased.                                                         */
01529 /*   Otherwise, the node is added to the list.  This list is arranged   */
01530 /*   so that bonds with wildcards are placed at the tail of the list.   */
01531 /*   This will guarantee that if we just do a linear search, we will    */
01532 /*   always find an exact match before a wildcard match.                */
01533 /*                                                                        */
01534 /************************************************************************/
01535 
01536 void Parameters::add_to_charmm_dihedral_list(
01537                                 struct dihedral_params *new_node)
01538 
01539 {
01540         static struct dihedral_params *ptr;   //  position within list
01541         static struct dihedral_params *tail;  //  Pointer to the end of 
01542                                               //  the list so we can add
01543                                               //  entries to the end of the
01544                                               //  list in constant time
01545         int i;                                      //  Loop counter
01546         int replace;                          //  replace values?
01547 
01548         // keep track of the last dihedral param read to avoid spurious
01549         // error messages.
01550         static struct dihedral_params last_dihedral; 
01551 
01552         /*  If the list is currently empty, then the new node is the list*/
01553         if (dihedralp == NULL)
01554         {
01555                 dihedralp=new_node;
01556                 tail=new_node;
01557                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01558 
01559                 return;
01560         }
01561 
01562         /*  The list isn't empty, so check for a duplicate                */
01563         ptr=dihedralp;
01564 
01565         while (ptr != NULL)
01566         {
01567                 int same_as_last = 0;
01568                 if (  ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01569                        (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01570                        (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01571                        (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01572                      ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01573                        (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01574                        (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01575                        (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) )
01576                        )
01577                 {
01578                         /*  Found a duplicate                                */
01579                         
01580                         // check for same_as_last.  Note: don't believe the echoWarn crap; it controls
01581                         // not just whether we print warning messages, but whether we actually change
01582                         // values or not!  
01583 
01584                         if ( ( !strcmp(ptr->atom1name, last_dihedral.atom1name) && 
01585                                !strcmp(ptr->atom2name, last_dihedral.atom2name) &&
01586                                !strcmp(ptr->atom3name, last_dihedral.atom3name) &&
01587                                !strcmp(ptr->atom4name, last_dihedral.atom4name)))
01588                           same_as_last = 1;
01589 
01590                         //****** BEGIN CHARMM/XPLOR type changes
01591                         /* we do not care about identical replacement */
01592                         int echoWarn=1;  // echo warning messages ?
01593 
01594                         // ptr->multiplicity will always be >= new_node->multiplicity
01595                         for (i=0; i<ptr->multiplicity; i++)
01596                         {
01597                           if ((ptr->values[i].k == new_node->values[0].k) && 
01598                               (ptr->values[i].n == new_node->values[0].n) &&
01599                               (ptr->values[i].delta == new_node->values[0].delta)) 
01600                           {
01601                             // found an identical replacement
01602                             echoWarn=0; 
01603                             break;
01604                           }
01605 
01606                         }
01607                   
01608                         if (echoWarn)
01609                         {
01610                           if (!same_as_last) {
01611                             iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
01612                                  << ptr->atom1name << "-"
01613                                  << ptr->atom2name << "-"
01614                                  << ptr->atom3name << "-"
01615                                  << ptr->atom4name
01616                                  << "\nPREVIOUS VALUES MULTIPLICITY: " << ptr->multiplicity << "\n";
01617                           }
01618                           replace=0;
01619                           
01620                           for (i=0; i<ptr->multiplicity; i++)
01621                           {
01622                             if (!same_as_last) {
01623                               iout << "  k=" << ptr->values[i].k
01624                                    << "  n=" << ptr->values[i].n
01625                                    << "  delta=" << ptr->values[i].delta << "\n";
01626                             }
01627                             if (ptr->values[i].n == new_node->values[0].n)
01628                             {
01629                               iout << iWARN << "IDENTICAL PERIODICITY! REPLACING OLD VALUES BY: \n";
01630                               ptr->values[i].k = new_node->values[0].k;
01631                               ptr->values[i].delta = new_node->values[0].delta;
01632                               iout << "  k=" << ptr->values[i].k
01633                                    << "  n=" << ptr->values[i].n
01634                                    << "  delta=" << ptr->values[i].delta<< "\n";
01635                               replace=1;
01636                               break;
01637                             }
01638                           }
01639 
01640                           if (!replace)
01641                           {
01642                             ptr->multiplicity += 1;
01643 
01644                             if (ptr->multiplicity > MAX_MULTIPLICITY)
01645                             {
01646                               char err_msg[181];
01647 
01648                               sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
01649                                       ptr->multiplicity, MAX_MULTIPLICITY);
01650                               NAMD_die(err_msg);
01651                             }
01652                             if (!same_as_last) 
01653                               iout << "INCREASING MULTIPLICITY TO: " << ptr->multiplicity << "\n";
01654 
01655                             i= ptr->multiplicity - 1; 
01656                             ptr->values[i].k = new_node->values[0].k;
01657                             ptr->values[i].n = new_node->values[0].n;
01658                             ptr->values[i].delta = new_node->values[0].delta;
01659 
01660                             if (!same_as_last) 
01661                               iout << "  k=" << ptr->values[i].k
01662                                    << "  n=" << ptr->values[i].n
01663                                    << "  delta=" << ptr->values[i].delta<< "\n";
01664                           }
01665                         
01666                           iout << endi;
01667                         } 
01668                         //****** END CHARMM/XPLOR type changes
01669 
01670                         memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01671                         delete new_node;
01672 
01673                         return;
01674                 }
01675 
01676                 ptr=ptr->next;
01677         }
01678 
01679         /*  CHARMM and XPLOR wildcards for dihedrals are luckily the same */
01680         /*  Check to see if we have any wildcards.  Since specific        */
01681         /*  entries are to take precedence, we'll put anything without  */
01682         /*  wildcards at the begining of the list and anything with     */
01683         /*  wildcards at the end of the list.  Then, we can just do a   */
01684         /*  linear search for a bond and be guaranteed to have specific */
01685         /*  entries take precendence over over wildcards                */
01686         if ( new_node->atom1wild ||
01687              new_node->atom2wild ||
01688              new_node->atom3wild ||
01689              new_node->atom4wild )
01690         {
01691                 /*  add to the end of the list                                */
01692                 tail->next=new_node;
01693                 tail=new_node;
01694 
01695                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01696                 return;
01697         }
01698         else
01699         {
01700                 /*  add to the head of the list                                */
01701                 new_node->next=dihedralp;
01702                 dihedralp=new_node;
01703 
01704                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01705                 return;
01706         }
01707 
01708 }
01709 /*                END OF FUNCTION add_to_charmm_dihedral_list                */
01710 //****** END CHARMM/XPLOR type changes
01711 
01712 /************************************************************************/
01713 /*                  */
01714 /*      FUNCTION add_improper_param      */
01715 /*                  */
01716 /*   INPUTS:                */
01717 /*  buf - line from paramter file with improper parameters    */
01718 /*                  */
01719 /*  this function adds an improper parameter.  It parses up the     */
01720 /*   input line and then adds it to the binary tree used to store the   */
01721 /*   improper parameters.            */
01722 /*                  */
01723 /************************************************************************/
01724 
01725 void Parameters::add_improper_param(char *buf, FILE *fd)
01726 
01727 {
01728   char atom1name[11];       //  Atom 1 type
01729   char atom2name[11];       //  Atom 2 type
01730   char atom3name[11];       //  Atom 3 type
01731   char atom4name[11];       //  Atom 4 type
01732   Real forceconstant;       //  Force constant 
01733   int periodicity;       //  Periodicity
01734   Real phase_shift;       //  Phase shift
01735   int read_count;         //  Count from sscanf
01736   struct improper_params *new_node;  //  New node
01737   int multiplicity;       //  Multiplicity for bonds
01738   int i;           //  Loop counter
01739   char buffer[513];       //  Buffer for new line
01740   int ret_code;         //  Return code
01741 
01742   //****** BEGIN CHARMM/XPLOR type changes
01743   /*  Parse up the line with sscanf                                */
01744   if (paramType == paraXplor)
01745   {
01746     /* read XPLOR format */
01747     read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n", 
01748        atom1name, atom2name, atom3name, atom4name, &multiplicity, 
01749        &forceconstant, &periodicity, &phase_shift);
01750   }
01751   else if (paramType == paraCharmm)
01752   {
01753     /* read CHARMM format */
01754     read_count=sscanf(buf, "%s %s %s %s %f %d %f\n", 
01755        atom1name, atom2name, atom3name, atom4name,  
01756        &forceconstant, &periodicity, &phase_shift); 
01757     multiplicity=1;      
01758   }
01759 
01760   if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
01761   {
01762     char err_msg[512];
01763 
01764     sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
01765     NAMD_die(err_msg);
01766   }
01767   else if ( (read_count != 7) && (paramType == paraCharmm) )
01768   {
01769     char err_msg[512];
01770 
01771     sprintf(err_msg, "BAD IMPROPER FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
01772     NAMD_die(err_msg);
01773   }
01774 
01775   if ( (read_count == 4) && (paramType == paraXplor) )
01776   //****** END CHARMM/XPLOR type changes
01777   {
01778     read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n", 
01779           &forceconstant, &periodicity, &phase_shift);
01780 
01781     /*  Check to make sure we got what we expected    */
01782     if (read_count != 3)
01783     {
01784       char err_msg[512];
01785 
01786       sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01787       NAMD_die(err_msg);
01788     }
01789 
01790     multiplicity = 1;
01791   }
01792 
01793   if (multiplicity > MAX_MULTIPLICITY)
01794   {
01795     char err_msg[181];
01796 
01797     sprintf(err_msg, "Multiple improper with multiplicity of %d greater than max of %d",
01798        multiplicity, MAX_MULTIPLICITY);
01799     NAMD_die(err_msg);
01800   }
01801 
01802   /*  Allocate a new node            */
01803   new_node = new improper_params;
01804 
01805   if (new_node == NULL)
01806   {
01807     NAMD_die("memory allocation failed in Parameters::add_improper_param");
01808   }
01809 
01810   /*  Assign the values for this bond.  As with the dihedrals,    */
01811   /*  the atom order doesn't matter        */
01812   strcpy(new_node->atom1name, atom1name);
01813   strcpy(new_node->atom2name, atom2name);
01814   strcpy(new_node->atom3name, atom3name);
01815   strcpy(new_node->atom4name, atom4name);
01816   new_node->multiplicity = multiplicity;
01817   if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01818   new_node->values[0].k = forceconstant;
01819   new_node->values[0].n = periodicity;
01820   new_node->values[0].delta = phase_shift;
01821 
01822   new_node->next = NULL;
01823 
01824   //  Check to see if this improper has multiple values
01825   if (multiplicity > 1)
01826   {
01827     //  Loop through and read the other values
01828     for (i=1; i<multiplicity; i++)
01829     {
01830       ret_code = NAMD_read_line(fd, buffer);
01831 
01832       //  Strip off comments at the end of the line
01833       if (ret_code == 0)
01834       {
01835         NAMD_remove_comment(buffer);
01836       }
01837 
01838       //  Skip blank lines
01839       while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
01840       {
01841         ret_code = NAMD_read_line(fd, buffer);
01842       }
01843 
01844       if (ret_code != 0)
01845       {
01846         NAMD_die("EOF encoutner in middle of multiple improper");
01847       }
01848 
01849       //  Get the values from the line
01850       read_count=sscanf(buffer, "%f %d %f\n", 
01851             &forceconstant, &periodicity, &phase_shift);
01852 
01853       if (read_count != 3)
01854       {
01855         char err_msg[512];
01856 
01857         sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
01858         NAMD_die(err_msg);
01859       }
01860 
01861       if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01862       new_node->values[i].k = forceconstant;
01863       new_node->values[i].n = periodicity;
01864       new_node->values[i].delta = phase_shift;
01865     }
01866   }
01867 
01868   /*  Add the paramter to the list        */
01869   add_to_improper_list(new_node);  // works for both XPLOR & CHARMM
01870 
01871   return;
01872 }
01873 /*      END OF FUNCTION add_improper_param    */
01874 
01875 /************************************************************************/
01876 /*                  */
01877 /*      FUNCTION add_to_improper_list      */
01878 /*                  */
01879 /*   INPUTS:                */
01880 /*  new_node - node that is to be added to imporper_list    */
01881 /*                  */
01882 /*  this function adds a new dihedral parameter to the linked list  */
01883 /*   of improper parameters.  First, it checks for duplicates.  If a    */
01884 /*   duplicate is found, a warning message is printed, the old values   */
01885 /*   are replaced with the new values, and the new node is freed.  If   */
01886 /*   Otherwise, the node is added to the list.  This list is arranged   */
01887 /*   so that bods with wildcards are placed at the tail of the list.    */
01888 /*   This will guarantee that if we just do a linear search, we will    */
01889 /*   always find an exact match before a wildcard match.    */
01890 /*                  */
01891 /************************************************************************/
01892 
01893 void Parameters::add_to_improper_list(struct improper_params *new_node)
01894 
01895 {
01896   int i;              //  Loop counter
01897   static struct improper_params *ptr;   //  position within list
01898   static struct improper_params *tail;  //  Pointer to the end of 
01899                 //  the list so we can add
01900                 //  entries to the end of the
01901                 //  list in constant time
01902 
01903   /*  If the list is currently empty, then the new node is the list*/
01904   if (improperp == NULL)
01905   {
01906     improperp=new_node;
01907     tail=new_node;
01908 
01909     return;
01910   }
01911 
01912   /*  The list isn't empty, so check for a duplicate    */
01913   ptr=improperp;
01914 
01915   while (ptr != NULL)
01916   {
01917     if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01918            (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01919            (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01920            (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01921          ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01922            (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01923            (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01924            (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
01925     {
01926       /*  Found a duplicate        */
01927       //****** BEGIN CHARMM/XPLOR type changes
01928       /* we do not care about identical replacement */
01929       int echoWarn=0;  // echo warning messages ?
01930 
01931       if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
01932       
01933       if (!echoWarn)
01934       {
01935         for (i=0; i<ptr->multiplicity; i++)
01936         {
01937           if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
01938           if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
01939           if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
01940         }
01941       }
01942 
01943       if (echoWarn)
01944       {
01945         iout << "\n" << iWARN << "DUPLICATE IMPROPER DIHEDRAL ENTRY FOR "
01946           << ptr->atom1name << "-"
01947           << ptr->atom2name << "-"
01948           << ptr->atom3name << "-"
01949           << ptr->atom4name
01950           << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
01951         
01952         for (i=0; i<ptr->multiplicity; i++)
01953         {
01954           iout <<     "  k=" << ptr->values[i].k
01955                    << "  n=" << ptr->values[i].n
01956                    << "  delta=" << ptr->values[i].delta;
01957         }
01958 
01959         iout << "\n" << "USING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
01960 
01961         for (i=0; i<new_node->multiplicity; i++)
01962         {
01963           iout <<     "  k=" << new_node->values[i].k
01964                    << "  n=" << new_node->values[i].n
01965                    << "  delta=" << new_node->values[i].delta;
01966         }
01967 
01968         iout << endi;
01969 
01970         ptr->multiplicity = new_node->multiplicity;
01971 
01972         for (i=0; i<new_node->multiplicity; i++)
01973         {
01974           ptr->values[i].k = new_node->values[i].k;
01975           ptr->values[i].n = new_node->values[i].n;
01976           ptr->values[i].delta = new_node->values[i].delta;
01977         }
01978       }
01979       //****** END CHARMM/XPLOR type changes
01980 
01981       delete new_node;
01982 
01983       return;
01984     }
01985 
01986     ptr=ptr->next;
01987   }
01988 
01989   /*  Check to see if we have any wildcards.  Since specific  */
01990   /*  entries are to take precedence, we'll put anything without  */
01991   /*  wildcards at the begining of the list and anything with     */
01992   /*  wildcards at the end of the list.  Then, we can just do a   */
01993   /*  linear search for a bond and be guaranteed to have specific */
01994   /*  entries take precendence over over wildcards          */
01995   if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
01996        (strcasecmp(new_node->atom2name, "X") == 0) ||
01997        (strcasecmp(new_node->atom3name, "X") == 0) ||
01998        (strcasecmp(new_node->atom4name, "X") == 0) )
01999   {
02000     /*  add to the end of the list        */
02001     tail->next=new_node;
02002     tail=new_node;
02003 
02004     return;
02005   }
02006   else
02007   {
02008     /*  add to the head of the list        */
02009     new_node->next=improperp;
02010     improperp=new_node;
02011 
02012     return;
02013   }
02014 }
02015 /*    END OF FUNCTION add_to_improper_list      */
02016 
02017 /************************************************************************/
02018 /*                  */
02019 /*      FUNCTION add_crossterm_param      */
02020 /*                  */
02021 /*   INPUTS:                */
02022 /*  buf - line from paramter file with crossterm parameters    */
02023 /*                  */
02024 /*  this function adds an crossterm parameter.  It parses up the     */
02025 /*   input line and then adds it to the binary tree used to store the   */
02026 /*   crossterm parameters.            */
02027 /*                  */
02028 /************************************************************************/
02029 
02030 void Parameters::add_crossterm_param(char *buf, FILE *fd)
02031 
02032 {
02033   char atom1name[11];       //  Atom 1 type
02034   char atom2name[11];       //  Atom 2 type
02035   char atom3name[11];       //  Atom 3 type
02036   char atom4name[11];       //  Atom 4 type
02037   char atom5name[11];       //  Atom 1 type
02038   char atom6name[11];       //  Atom 2 type
02039   char atom7name[11];       //  Atom 3 type
02040   char atom8name[11];       //  Atom 4 type
02041   int dimension;
02042   int read_count;         //  Count from sscanf
02043   struct crossterm_params *new_node;  //  New node
02044   char buffer[513];       //  Buffer for new line
02045   int ret_code;         //  Return code
02046 
02047   /* read CHARMM format */
02048   read_count=sscanf(buf, "%s %s %s %s %s %s %s %s %d\n", 
02049      atom1name, atom2name, atom3name, atom4name,  
02050      atom5name, atom6name, atom7name, atom8name,  
02051      &dimension); 
02052 
02053   if ( (read_count != 9) || dimension < 1 || dimension > 1000 )
02054   {
02055     char err_msg[512];
02056 
02057     sprintf(err_msg, "BAD CMAP FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
02058     NAMD_die(err_msg);
02059   }
02060 
02061   /*  Allocate a new node            */
02062   new_node = new crossterm_params(dimension);
02063 
02064   /*  Assign the values for this bond.  As with the dihedrals,    */
02065   /*  the atom order doesn't matter        */
02066   strcpy(new_node->atom1name, atom1name);
02067   strcpy(new_node->atom2name, atom2name);
02068   strcpy(new_node->atom3name, atom3name);
02069   strcpy(new_node->atom4name, atom4name);
02070   strcpy(new_node->atom5name, atom5name);
02071   strcpy(new_node->atom6name, atom6name);
02072   strcpy(new_node->atom7name, atom7name);
02073   strcpy(new_node->atom8name, atom8name);
02074 
02075   new_node->next = NULL;
02076 
02077   int nterms = dimension * dimension;
02078   int nread = 0;
02079 
02080   //  Loop through and read the other values
02081   while ( nread < nterms ) {
02082     ret_code = NAMD_read_line(fd, buffer);
02083 
02084     //  Strip off comments at the end of the line
02085     if (ret_code == 0) {
02086       NAMD_remove_comment(buffer);
02087     }
02088 
02089     //  Skip blank lines
02090     while ( (ret_code == 0) && (NAMD_blank_string(buffer)) ) {
02091       ret_code = NAMD_read_line(fd, buffer);
02092       if (ret_code == 0) {
02093         NAMD_remove_comment(buffer);
02094       }
02095     }
02096 
02097     if (ret_code != 0) {
02098       NAMD_die("EOF encoutner in middle of CMAP");
02099     }
02100 
02101     //  Get the values from the line
02102     read_count=sscanf(buffer, "%lf %lf %lf %lf %lf %lf %lf %lf\n",
02103         new_node->values + nread,
02104         new_node->values + nread+1,
02105         new_node->values + nread+2,
02106         new_node->values + nread+3,
02107         new_node->values + nread+4,
02108         new_node->values + nread+5,
02109         new_node->values + nread+6,
02110         new_node->values + nread+7);
02111 
02112     nread += read_count;
02113 
02114     if (read_count == 0 || nread > nterms) {
02115       char err_msg[512];
02116 
02117       sprintf(err_msg, "BAD CMAP FORMAT IN PARAMETER FILE\nLINE=*%s*\n", buffer);
02118       NAMD_die(err_msg);
02119     }
02120   }
02121 
02122   /*  Add the paramter to the list        */
02123   add_to_crossterm_list(new_node);
02124 
02125   return;
02126 }
02127 /*      END OF FUNCTION add_crossterm_param    */
02128 
02129 /************************************************************************/
02130 /*                  */
02131 /*      FUNCTION add_to_crossterm_list      */
02132 /*                  */
02133 /*   INPUTS:                */
02134 /*  new_node - node that is to be added to crossterm_list    */
02135 /*                  */
02136 /*  this function adds a new crossterm parameter to the linked list  */
02137 /*   of crossterm parameters.  First, it checks for duplicates.  If a    */
02138 /*   duplicate is found, a warning message is printed, the old values   */
02139 /*   are replaced with the new values, and the new node is freed.  If   */
02140 /*   Otherwise, the node is added to the list.  This list is arranged   */
02141 /*   so that bods with wildcards are placed at the tail of the list.    */
02142 /*   This will guarantee that if we just do a linear search, we will    */
02143 /*   always find an exact match before a wildcard match.    */
02144 /*                  */
02145 /************************************************************************/
02146 
02147 void Parameters::add_to_crossterm_list(struct crossterm_params *new_node)
02148 
02149 {
02150   int i;              //  Loop counter
02151   static struct crossterm_params *ptr;   //  position within list
02152   static struct crossterm_params *tail;  //  Pointer to the end of 
02153                 //  the list so we can add
02154                 //  entries to the end of the
02155                 //  list in constant time
02156 
02157   /*  If the list is currently empty, then the new node is the list*/
02158   if (crosstermp == NULL)
02159   {
02160     crosstermp=new_node;
02161     tail=new_node;
02162 
02163     return;
02164   }
02165 
02166   /*  The list isn't empty, so check for a duplicate    */
02167   ptr=crosstermp;
02168 
02169   while (ptr != NULL)
02170   {
02171     if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
02172            (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
02173            (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
02174            (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) &&
02175            (strcasecmp(new_node->atom5name, ptr->atom5name) == 0) &&
02176            (strcasecmp(new_node->atom6name, ptr->atom6name) == 0) &&
02177            (strcasecmp(new_node->atom7name, ptr->atom7name) == 0) &&
02178            (strcasecmp(new_node->atom8name, ptr->atom8name) == 0) ) )
02179     {
02180       /*  Found a duplicate        */
02181       /* we do not care about identical replacement */
02182       int echoWarn=0;  // echo warning messages ?
02183 
02184       if (ptr->dimension != new_node->dimension) {echoWarn=1;}
02185       
02186       if (!echoWarn)
02187       {
02188         int nvals = ptr->dimension * ptr->dimension;
02189         for (i=0; i<nvals; i++)
02190         {
02191           if (ptr->values[i] != new_node->values[i]) {echoWarn=1; break;}
02192         }
02193       }
02194 
02195       if (echoWarn)
02196       {
02197         iout << "\n" << iWARN << "DUPLICATE CMAP ENTRY FOR "
02198           << ptr->atom1name << "-"
02199           << ptr->atom2name << "-"
02200           << ptr->atom3name << "-"
02201           << ptr->atom4name << " "
02202           << ptr->atom5name << "-"
02203           << ptr->atom6name << "-"
02204           << ptr->atom7name << "-"
02205           << ptr->atom8name << ", USING NEW VALUES\n";
02206         
02207         iout << endi;
02208 
02209         ptr->dimension = new_node->dimension;
02210 
02211         BigReal *tmpvalues = ptr->values;
02212         ptr->values = new_node->values;
02213         new_node->values = tmpvalues;
02214       }
02215 
02216       delete new_node;
02217 
02218       return;
02219     }
02220 
02221     ptr=ptr->next;
02222   }
02223 
02224   /*  Check to see if we have any wildcards.  Since specific  */
02225   /*  entries are to take precedence, we'll put anything without  */
02226   /*  wildcards at the begining of the list and anything with     */
02227   /*  wildcards at the end of the list.  Then, we can just do a   */
02228   /*  linear search for a bond and be guaranteed to have specific */
02229   /*  entries take precendence over over wildcards          */
02230   if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
02231        (strcasecmp(new_node->atom2name, "X") == 0) ||
02232        (strcasecmp(new_node->atom3name, "X") == 0) ||
02233        (strcasecmp(new_node->atom4name, "X") == 0) ||
02234        (strcasecmp(new_node->atom5name, "X") == 0) ||
02235        (strcasecmp(new_node->atom6name, "X") == 0) ||
02236        (strcasecmp(new_node->atom7name, "X") == 0) ||
02237        (strcasecmp(new_node->atom8name, "X") == 0) )
02238   {
02239     /*  add to the end of the list        */
02240     tail->next=new_node;
02241     tail=new_node;
02242 
02243     return;
02244   }
02245   else
02246   {
02247     /*  add to the head of the list        */
02248     new_node->next=crosstermp;
02249     crosstermp=new_node;
02250 
02251     return;
02252   }
02253 }
02254 /*    END OF FUNCTION add_to_crossterm_list      */
02255 
02256 /************************************************************************/
02257 /*                  */
02258 /*      FUNCTION add_vdw_param        */
02259 /*                  */
02260 /*  INPUTS:                */
02261 /*  buf - line containing the vdw information      */
02262 /*                  */
02263 /*  add_vdw_param adds a vdw parameter for an atom to the current   */
02264 /*  binary tree of values.            */
02265 /*                  */
02266 /************************************************************************/
02267 
02268 void Parameters::add_vdw_param(char *buf)
02269 
02270 {
02271   char atomname[11];    //  atom type of paramter
02272   Real sigma;      //  sigma value for this atom
02273   Real epsilon;      //  epsilon value for this atom
02274   Real sigma14;      //  sigma value for 1-4 interactions
02275   Real epsilon14;      //  epsilon value for 1-4 interactions
02276   Real sqrt26;         //  2^(1/6)
02277   int read_count;      //  count returned by sscanf
02278   struct vdw_params *new_node;  //  new node for tree
02279 
02280   //****** BEGIN CHARMM/XPLOR type changes
02281   /*  Parse up the line with sscanf        */
02282   if (paramType == paraXplor)
02283   {
02284     /* read XPLOR format */
02285     read_count=sscanf(buf, "%*s %s %f %f %f %f\n", atomname, 
02286        &epsilon, &sigma, &epsilon14, &sigma14);
02287   }
02288   else if (paramType == paraCharmm)
02289   {
02290     /* read CHARMM format */
02291     read_count=sscanf(buf, "%s %*f %f %f %*f %f %f\n", atomname, 
02292        &epsilon, &sigma, &epsilon14, &sigma14);
02293   }
02294 
02295   /*  Check to make sure we got what we expected      */
02296   if ((read_count != 5) && (paramType == paraXplor))
02297   {
02298     char err_msg[512];
02299 
02300     sprintf(err_msg, "BAD vdW FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
02301     NAMD_die(err_msg);
02302   }
02303   else if ( ((read_count != 5) && (read_count != 3)) && (paramType == paraCharmm))
02304   {
02305     char err_msg[512];
02306 
02307     sprintf(err_msg, "BAD vdW FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
02308     NAMD_die(err_msg);
02309   }
02310 
02311   if (paramType == paraCharmm)
02312   {
02313     // convert CHARMM to XPLOR format
02314     epsilon*=-1.;
02315     sqrt26=pow(2.,(1./6.));
02316     sigma=2.*sigma/sqrt26; 
02317 
02318     if (read_count == 3)
02319     {
02320       epsilon14=epsilon;
02321       sigma14=sigma;
02322     }
02323     else
02324     {
02325       epsilon14*=-1.;
02326       sigma14=2.*sigma14/sqrt26; 
02327     }
02328   }
02329   //****** END CHARMM/XPLOR type changes
02330 
02331   if ( epsilon < 0. || epsilon14 < 0. ) {
02332     iout << iWARN << "Ignoring VDW parameter with negative epsilon:\n"
02333         << buf << "\n" << endi;
02334     return;
02335   }
02336 
02337   /*  Allocate a new node            */
02338   new_node = new vdw_params;
02339 
02340   if (new_node == NULL)
02341   {
02342     NAMD_die("memory allocation failed in Parameters::add_vdw_param");
02343   }
02344 
02345   /*  Assign the values to the new node        */
02346   strcpy(new_node->atomname, atomname);
02347   new_node->sigma = sigma;
02348   new_node->sigma14 = sigma14;
02349   new_node->epsilon = epsilon;
02350   new_node->epsilon14 = epsilon14;
02351 
02352   new_node->left = NULL;
02353   new_node->right = NULL;
02354 
02355   /*  Add the new node into the tree        */
02356   vdwp=add_to_vdw_tree(new_node, vdwp);
02357 
02358   return;
02359 }
02360 /*      END OF FUNCTION add_vdw_param      */
02361 
02362 /************************************************************************/
02363 /*                  */
02364 /*      FUNCTION add_to_vdw_tree      */
02365 /*                  */
02366 /*   INPUTS:                */
02367 /*  new_node - node to add to tree          */
02368 /*  tree - tree to add the node to          */
02369 /*                  */
02370 /*   OUTPUTS:                */
02371 /*  the function returns a pointer to the tree with the node added  */
02372 /*                  */
02373 /*  this function adds a vdw to the binary tree containing the      */
02374 /*   parameters.              */
02375 /*                  */
02376 /************************************************************************/
02377 
02378 struct vdw_params *Parameters::add_to_vdw_tree(struct vdw_params *new_node,
02379              struct vdw_params *tree)
02380 
02381 {
02382   int compare_code;  //  Return code from strcasecmp
02383 
02384   /*  If the tree is currently empty, the new node is the tree    */
02385   if (tree == NULL)
02386     return(new_node);
02387 
02388   compare_code = strcasecmp(new_node->atomname, tree->atomname);
02389 
02390   /*  Check to see if we have a duplicate        */
02391   if (compare_code==0)
02392   {
02393     /*  We have a duplicate.  So print out a warning   */
02394     /*  message, copy the new values into the current node  */
02395     /*  of the tree, and then free the new_node    */
02396     if ((tree->sigma != new_node->sigma) || 
02397         (tree->epsilon != new_node->epsilon) ||
02398         (tree->sigma14 != new_node->sigma14) ||
02399         (tree->epsilon14 != new_node->epsilon14))
02400     {
02401       iout << iWARN << "DUPLICATE vdW ENTRY FOR " << tree->atomname
02402         << "\nPREVIOUS VALUES  sigma=" << tree->sigma
02403         << " epsilon=" << tree->epsilon
02404         << " sigma14=" << tree->sigma14
02405         << " epsilon14=" << tree->epsilon14
02406         << "\n   USING VALUES  sigma=" << new_node->sigma
02407         << " epsilon=" << new_node->epsilon
02408         << " sigma14=" << new_node->sigma14
02409         << " epsilon14=" << new_node->epsilon14
02410         << "\n" << endi;
02411 
02412       tree->sigma=new_node->sigma;
02413       tree->epsilon=new_node->epsilon;
02414       tree->sigma14=new_node->sigma14;
02415       tree->epsilon14=new_node->epsilon14;
02416     }
02417 
02418     delete new_node;
02419 
02420     return(tree);
02421   }
02422 
02423   /*  Otherwise, if the new node is less than the head of    */
02424   /*  the tree, add it to the left child, and if it is greater  */
02425   /*  add it to the right child          */
02426   if (compare_code < 0)
02427   {
02428     tree->left = add_to_vdw_tree(new_node, tree->left);
02429   }
02430   else
02431   {
02432     tree->right = add_to_vdw_tree(new_node, tree->right);
02433   }
02434 
02435   return(tree);
02436 }
02437 /*      END OF FUNCTION add_to_vdw_tree      */
02438 
02439 /************************************************************************/
02440 /*                  */
02441 /*      FUNCTION add_table_pair_param      */
02442 /*                  */
02443 /*   INPUTS:                */
02444 /*  buf - line containing the table_pair information      */
02445 /*                  */
02446 /*  this function adds a table_pair parameter to the current          */
02447 /*   parameters.              */
02448 /*                  */
02449 /************************************************************************/
02450 
02451 void Parameters::add_table_pair_param(char *buf)
02452 
02453 {
02454   char atom1name[11];      //  Atom 1 name
02455   char atom2name[11];      //  Atom 2 name
02456   char tabletype[11];      // Name of pair interaction
02457   int K;           // Table entry type for pair
02458   int read_count;        //  count from sscanf
02459   struct table_pair_params *new_node;  //  new node
02460 
02461   /*  Parse up the input line using sscanf      */
02462   read_count=sscanf(buf, "%s %s %s\n", atom1name, 
02463   atom2name, tabletype);
02464 
02465   /*  Check to make sure we got what we expected      */
02466   if ((read_count != 3))
02467   {
02468     char err_msg[512];
02469 
02470     sprintf(err_msg, "BAD TABLE PAIR FORMAT IN PARAMETER FILE\nLINE=*%s*", buf);
02471     NAMD_die(err_msg);
02472   }
02473 
02474   /*  Allocate a new node            */
02475   new_node = new table_pair_params;
02476 
02477   if (new_node == NULL)
02478   {
02479     NAMD_die("memory allocation failed in Parameters::add_table_pair_param\n");
02480   }
02481 
02482   strcpy(new_node->atom1name, atom1name);
02483   strcpy(new_node->atom2name, atom2name);
02484 
02485   /* Find the proper table type for this pairing */
02486   K = get_int_table_type(tabletype);
02487 //  printf("Looking for type %s; got %i\n", tabletype, K);
02488   if (K < 0) {
02489     char err_msg[512];
02490     sprintf(err_msg, "Couldn't find table parameters for table interaction %s!\n", tabletype);
02491     NAMD_die(err_msg);
02492   }
02493 
02494   /*  Assign values to this node          */
02495   new_node->K = K;
02496 
02497   new_node->next = NULL;
02498 
02499   /*  Add this node to the tree          */
02500 //  printf("Adding pair parameter with index %i\n", K);
02501   add_to_table_pair_list(new_node);
02502 
02503   return;
02504 }
02505 /*      END OF FUNCTION add_vdw_par_param    */
02506 
02507 /************************************************************************/
02508 /*                  */
02509 /*      FUNCTION add_vdw_pair_param      */
02510 /*                  */
02511 /*   INPUTS:                */
02512 /*  buf - line containing the vdw_pair information      */
02513 /*                  */
02514 /*  this function adds a vdw_pair parameter to the current          */
02515 /*   parameters.              */
02516 /*                  */
02517 /************************************************************************/
02518 
02519 void Parameters::add_vdw_pair_param(char *buf)
02520 
02521 {
02522   char atom1name[11];      //  Atom 1 name
02523   char atom2name[11];      //  Atom 2 name
02524   Real A;          //  A value for pair
02525   Real B;          //  B value for pair
02526   Real A14;        //  A value for 1-4 ints
02527   Real B14;        //  B value for 1-4 ints
02528   Real sqrt26;     //  2^(1/6)
02529   Real expo;       //  just for pow
02530   int read_count;        //  count from sscanf
02531   struct vdw_pair_params *new_node;  //  new node
02532 
02533   /*  Parse up the input line using sscanf      */
02534   if (paramType == paraXplor)
02535   {
02536     /* read XPLOR format */
02537     read_count=sscanf(buf, "%*s %s %s %f %f %f %f\n", atom1name, 
02538        atom2name, &A, &B, &A14, &B14);
02539   }
02540   else if (paramType == paraCharmm)
02541   {
02542     // XXX CHARMM CAN HAVE 1-4 PARAMETERS TOO!!!
02543     /* read CHARMM format */
02544     read_count=sscanf(buf, "%s %s %f %f\n", atom1name, 
02545        atom2name, &A, &B);
02546     // convert to XPLOR format and use A14, B14 as dummies
02547     /* XXX doesn't this just do the following?
02548        A = -eps*pow(sig,6.);
02549        B = -2*eps*pow(sig,12.);
02550        Why do we need to do this in such a complicated way?
02551        -pgrayson */
02552     A14=-A;
02553     sqrt26=pow(2.,(1./6.));
02554     B14=B/sqrt26;
02555     expo=12.;
02556     A=pow(B14,expo);
02557     A=A*4.*A14;
02558     expo=6.;
02559     B=pow(B14,expo);
02560     B=B*4.*A14;
02561     A14=A;
02562     B14=B;
02563   }
02564 
02565   /*  Check to make sure we got what we expected      */
02566   if ((read_count != 6) && (paramType == paraXplor))
02567   {
02568     char err_msg[512];
02569 
02570     sprintf(err_msg, "BAD vdW PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
02571     NAMD_die(err_msg);
02572   }
02573   if ((read_count != 4) && (paramType == paraCharmm))
02574   {
02575     char err_msg[512];
02576 
02577     sprintf(err_msg, "BAD vdW PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
02578     NAMD_die(err_msg);
02579   }
02580 
02581 
02582   /*  Allocate a new node            */
02583   new_node = new vdw_pair_params;
02584 
02585   if (new_node == NULL)
02586   {
02587     NAMD_die("memory allocation failed in Parameters::add_vdw_pair_param\n");
02588   }
02589 
02590   strcpy(new_node->atom1name, atom1name);
02591   strcpy(new_node->atom2name, atom2name);
02592 
02593   /*  Assign values to this node          */
02594   new_node->A = A;
02595   new_node->A14 = A14;
02596   new_node->B = B;
02597   new_node->B14 = B14;
02598 
02599   new_node->next = NULL;
02600 
02601   /*  Add this node to the tree          */
02602   add_to_vdw_pair_list(new_node);
02603 
02604   return;
02605 }
02606 /*      END OF FUNCTION add_vdw_par_param    */
02607 
02608 /************************************************************************/
02609 /*                  */
02610 /*      FUNCTION add_hb_pair_param      */
02611 /*                  */
02612 /*   INPUTS:                */
02613 /*  buf - line containing the hydrogen bond information    */
02614 /*                  */
02615 /*  this function adds data for a hydrogen bond interaction pair    */
02616 /*   to the hbondParams object.                                         */
02617 /*                  */
02618 /************************************************************************/
02619 
02620 void Parameters::add_hb_pair_param(char *buf)
02621 
02622 {
02623 #if 0
02624   char a1n[11];      //  Atom 1 name
02625   char a2n[11];      //  Atom 2 name
02626   Real A, B;      //  A, B value for pair
02627 
02628   //****** BEGIN CHARMM/XPLOR type changes
02630   /*  Parse up the input line using sscanf      */
02631   if (paramType == paraXplor) {
02632     if (sscanf(buf, "%*s %s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
02633       char err_msg[512];
02634       sprintf(err_msg, "BAD HBOND PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
02635       NAMD_die(err_msg);
02636     }
02637   }
02638   else if (paramType == paraCharmm) {
02639     if (sscanf(buf, "%s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
02640       char err_msg[512];
02641       sprintf(err_msg, "BAD HBOND PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
02642       NAMD_die(err_msg);
02643     }
02644   }
02645   //****** END CHARMM/XPLOR type changes
02646 
02647   /*  add data */
02648   if (hbondParams.add_hbond_pair(a1n, a2n, A, B) == FALSE) {
02649     iout << "\n" << iWARN << "Duplicate HBOND parameters for types " << a1n
02650     << " and " << a2n << " found; using latest values." << "\n" << endi;
02651   }
02652 #endif
02653 }
02654 /*      END OF FUNCTION add_hb_par_param    */
02655 
02656 /************************************************************************/
02657 /*                  */
02658 /*      FUNCTION add_to_table_pair_list      */
02659 /*                  */
02660 /*   INPUTS:                */
02661 /*  new_node - node to be added to list        */
02662 /*                  */
02663 /*  This function adds a link to the end of the table_pair_list list  */
02664 /*                  */
02665 /************************************************************************/
02666 
02667 void Parameters::add_to_table_pair_list(struct table_pair_params *new_node)
02668 
02669 {
02670      static struct table_pair_params *tail=NULL;
02671   struct table_pair_params *ptr;
02672   int compare_code;
02673   
02674 
02675   //  If the list was empty, then just make the new node the list
02676   if (table_pairp == NULL)
02677   {
02678      table_pairp = new_node;
02679      tail = new_node;
02680      return;
02681   }
02682   
02683   ptr = table_pairp;
02684 
02685   //  Now check the list to see if we have a duplicate entry
02686   while (ptr!=NULL)
02687   {
02688       /*  Compare atom 1            */
02689       compare_code = strncasecmp(new_node->atom1name, ptr->atom1name, 5);
02690       
02691       if (compare_code == 0)
02692       {
02693     /*  Atom 1 is the same, compare atom 2      */
02694     compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
02695 
02696     if (compare_code==0)
02697     {
02698       /*  Found a duplicate.  Print out a warning   */
02699       /*  message, assign the values to the current   */
02700       /*  node in the tree, and then free the new_node*/
02701       iout << iWARN << "DUPLICATE TABLE PAIR ENTRY FOR "
02702         << new_node->atom1name << "-"
02703         << new_node->atom2name
02704         << "\n" << endi;
02705 
02706         ptr->K=new_node->K;
02707 
02708       delete new_node;
02709 
02710       return;
02711     }
02712       }
02713       
02714       ptr = ptr->next;
02715   }
02716 
02717   //  We didn't find a duplicate, so add this node to the end
02718   //  of the list
02719   tail->next = new_node;
02720   tail = new_node;
02721 }
02722 /*      END OF FUNCTION add_to_vdw_pair_list    */
02723 
02724 /************************************************************************/
02725 /*                  */
02726 /*      FUNCTION add_to_vdw_pair_list      */
02727 /*                  */
02728 /*   INPUTS:                */
02729 /*  new_node - node to be added to list        */
02730 /*                  */
02731 /*  This function adds a link to the end of the vdw_pair_list list  */
02732 /*                  */
02733 /************************************************************************/
02734 
02735 void Parameters::add_to_vdw_pair_list(struct vdw_pair_params *new_node)
02736 
02737 {
02738      static struct vdw_pair_params *tail=NULL;
02739   struct vdw_pair_params *ptr;
02740   int compare_code;
02741   
02742 
02743   //  If the list was empty, then just make the new node the list
02744   if (vdw_pairp == NULL)
02745   {
02746      vdw_pairp = new_node;
02747      tail = new_node;
02748      return;
02749   }
02750   
02751   ptr = vdw_pairp;
02752 
02753   //  Now check the list to see if we have a duplicate entry
02754   while (ptr!=NULL)
02755   {
02756       /*  Compare atom 1            */
02757       compare_code = strcasecmp(new_node->atom1name, ptr->atom1name);
02758       
02759       if (compare_code == 0)
02760       {
02761     /*  Atom 1 is the same, compare atom 2      */
02762     compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
02763 
02764     if (compare_code==0)
02765     {
02766       /*  Found a duplicate.  Print out a warning   */
02767       /*  message, assign the values to the current   */
02768       /*  node in the tree, and then free the new_node*/
02769       if ((ptr->A != new_node->A) ||
02770           (ptr->B != new_node->B) ||
02771           (ptr->A14 != new_node->A14) ||
02772           (ptr->B14 != new_node->B14))
02773       {
02774         iout << iWARN << "DUPLICATE vdW PAIR ENTRY FOR "
02775           << new_node->atom1name << "-"
02776           << new_node->atom2name
02777           << "\nPREVIOUS VALUES  A=" << ptr->A
02778           << " B=" << ptr->B
02779           << " A14=" << ptr->A14
02780           << " B14" << ptr->B14
02781           << "\n   USING VALUES  A=" << new_node->A
02782           << " B=" << new_node->B
02783           << " A14=" << new_node->A14
02784           << " B14" << new_node->B14
02785           << "\n" << endi;
02786 
02787         ptr->A=new_node->A;
02788         ptr->B=new_node->B;
02789         ptr->A14=new_node->A14;
02790         ptr->B14=new_node->B14;
02791       }
02792 
02793       delete new_node;
02794 
02795       return;
02796     }
02797       }
02798       
02799       ptr = ptr->next;
02800   }
02801 
02802   //  We didn't find a duplicate, so add this node to the end
02803   //  of the list
02804   tail->next = new_node;
02805   tail = new_node;
02806 }
02807 /*      END OF FUNCTION add_to_vdw_pair_list    */
02808 
02809 /************************************************************************/
02810 /*                  */
02811 /*      FUNCTION done_reading_files      */
02812 /*                  */
02813 /*  This function is used to signal the Parameters object that all  */
02814 /*  of the parameter files have been read.  Once the object knows this, */
02815 /*  it can set un indexes for all the parameters and transfer the values*/
02816 /*  to linear arrays.  This will allow constant time access from this   */
02817 /*  point on.                */
02818 /*                  */
02819 /************************************************************************/
02820 
02821 void Parameters::done_reading_files()
02822 
02823 {
02824   AllFilesRead = TRUE;
02825 
02826   //  Allocate space for all of the arrays
02827   if (NumBondParams)
02828   {
02829     bond_array = new BondValue[NumBondParams];
02830 
02831     if (bond_array == NULL)
02832     {
02833       NAMD_die("memory allocation of bond_array failed!");
02834     }
02835   }
02836 
02837   if (NumAngleParams)
02838   {
02839     angle_array = new AngleValue[NumAngleParams];
02840 
02841     if (angle_array == NULL)
02842     {
02843       NAMD_die("memory allocation of angle_array failed!");
02844     }
02845   }
02846 
02847   if (NumDihedralParams)
02848   {
02849     dihedral_array = new DihedralValue[NumDihedralParams];
02850 
02851     if (dihedral_array == NULL)
02852     {
02853       NAMD_die("memory allocation of dihedral_array failed!");
02854     }
02855     memset(dihedral_array, 0, NumDihedralParams*sizeof(DihedralValue));
02856   }
02857 
02858   if (NumImproperParams)
02859   {
02860     improper_array = new ImproperValue[NumImproperParams];
02861 
02862     if (improper_array == NULL)
02863     {
02864       NAMD_die("memory allocation of improper_array failed!");
02865     }
02866     memset(improper_array, 0, NumImproperParams*sizeof(ImproperValue));
02867   }
02868 
02869   if (NumCrosstermParams)
02870   {
02871     crossterm_array = new CrosstermValue[NumCrosstermParams];
02872     memset(crossterm_array, 0, NumCrosstermParams*sizeof(CrosstermValue));
02873   }
02874 
02875   if (NumVdwParams)
02876   {
02877           atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
02878     vdw_array = new VdwValue[NumVdwParams];
02879     
02880     if (vdw_array == NULL)
02881     {
02882       NAMD_die("memory allocation of vdw_array failed!");
02883     }
02884   }
02885 
02886   //  Assign indexes to each of the parameters and populate the
02887   //  arrays using the binary trees and linked lists that we have
02888   //  already read in
02889   index_bonds(bondp, 0);
02890   index_angles(anglep, 0);
02891   NumVdwParamsAssigned = index_vdw(vdwp, 0);
02892   index_dihedrals();
02893   index_impropers();
02894   index_crossterms();
02895   
02896   //  Convert the vdw pairs
02897   convert_vdw_pairs();
02898   convert_table_pairs();
02899 }
02900 /*      END OF FUNCTION done_reading_files    */
02901 
02902 /************************************************************************/
02903 /*                  */
02904 /*      FUNCTION index_bonds        */
02905 /*                  */
02906 /*   INPUTS:                */
02907 /*  tree - The tree that is to be indexed        */
02908 /*  index - index to start with          */
02909 /*                  */
02910 /*  This is a recursive routine that will traverse the binary tree  */
02911 /*   of bond parameters, assigning an index to each one, and copying    */
02912 /*   the data from the binary tree to the array that will be used from  */
02913 /*   here on.                */
02914 /*                  */
02915 /************************************************************************/
02916 
02917 Index Parameters::index_bonds(struct bond_params *tree, Index index)
02918 
02919 {
02920   //  Tree is empty, do nothing
02921   if (tree==NULL)
02922     return(index);
02923 
02924   //  If I have a left subtree, index it first
02925   if (tree->left != NULL)
02926   {
02927     index=index_bonds(tree->left, index);
02928   }
02929 
02930   //  Now assign an index to top node and populate array
02931   tree->index = index;
02932   bond_array[index].k = tree->forceconstant;
02933   bond_array[index].x0 = tree->distance;
02934   index++;
02935 
02936   //  If I have a right subtree, index it
02937   if (tree->right != NULL)
02938   {
02939     index=index_bonds(tree->right, index);
02940   }
02941 
02942   return(index);
02943 }
02944 /*      END OF FUNCTION index_bonds      */
02945 
02946 /************************************************************************/
02947 /*                  */
02948 /*      FUNCTION index_angles        */
02949 /*                  */
02950 /*   INPUTS:                */
02951 /*  tree - The tree that is to be indexed        */
02952 /*  index - index to start with          */
02953 /*                  */
02954 /*  This is a recursive routine that will traverse the binary tree  */
02955 /*   of angle parameters, assigning an index to each one, and copying   */
02956 /*   the data from the binary tree to the array that will be used from  */
02957 /*   here on.                */
02958 /*                  */
02959 /************************************************************************/
02960 
02961 Index Parameters::index_angles(struct angle_params *tree, Index index)
02962 
02963 {
02964   //  Tree is empty, do nothing
02965   if (tree==NULL)
02966     return(index);
02967 
02968   //  If I have a left subtree, index it first
02969   if (tree->left != NULL)
02970   {
02971     index=index_angles(tree->left, index);
02972   }
02973 
02974   //  Now assign an index to top node and populate array
02975   tree->index = index;
02976 
02977   angle_array[index].k = tree->forceconstant;
02978   angle_array[index].k_ub = tree->k_ub;
02979   angle_array[index].r_ub = tree->r_ub;
02980   angle_array[index].normal = tree->normal;
02981 
02982   //  Convert the angle to radians before storing it
02983   angle_array[index].theta0 = (tree->angle*PI)/180.0;
02984   index++;
02985 
02986   //  If I have a right subtree, index it
02987   if (tree->right != NULL)
02988   {
02989     index=index_angles(tree->right, index);
02990   }
02991 
02992   return(index);
02993 }
02994 /*      END OF FUNCTION index_angles      */
02995 
02996 /************************************************************************/
02997 /*                  */
02998 /*      FUNCTION index_dihedrals      */
02999 /*                  */
03000 /*  This function walks down the linked list of dihedral parameters */
03001 /*  and assigns an index to each one.  It also copies the data from this*/
03002 /*  linked list to the arrays that will be used from here on out  */
03003 /*                  */
03004 /************************************************************************/
03005 
03006 void Parameters::index_dihedrals()
03007 
03008 {
03009   struct dihedral_params *ptr;  //  Current location in list
03010   Index index=0;      //  Current index value
03011   int i;        //  Loop counter
03012 
03013   //  Allocate an array to hold the multiplicity present in the
03014   //  parameter file for each bond.  This will be used to check
03015   //  the multiplicities that are detected in the psf file
03016 
03017   //  This is kind of ugly, but necessary because of the way that
03018   //  X-PLOR psf files deal with Charmm22 parameters.  The way
03019   //  that multiple periodicities are specified is by having
03020   //  the bonds appear multiple times in the psf file.  This even
03021   //  if a bond type has multiple parameters defined, they
03022   //  will be used if the bond appears multiple times in the
03023   //  psf file.  So we need to store the number of parameters
03024   //  we have to make sure the psf file doesn't ask for more
03025   //  parameters than we really have, and we also need to track
03026   //  how many times the bond appears in the psf file so that
03027   //  we can decide how many parameters to actually use.
03028   //  This is different for CHARMM parameter files as stated below!
03029   maxDihedralMults = new int[NumDihedralParams];
03030 
03031   if (maxDihedralMults == NULL)
03032   {
03033     NAMD_die("memory allocation failed in Parameters::index_dihedrals()");
03034   }
03035   
03036   //  Start at the head
03037   ptr = dihedralp;
03038 
03039   while (ptr != NULL)
03040   {
03041     //  Copy data to array and assign index
03042 
03043     //  Save the multiplicity in another array
03044     maxDihedralMults[index] = ptr->multiplicity;
03045 
03046 
03047     //****** BEGIN CHARMM/XPLOR type changes
03048     if (paramType == paraXplor)
03049     {
03050       //  Assign the multiplicity in the actual structure a bogus value
03051       //  that we will update in assign_dihedral_index
03052       dihedral_array[index].multiplicity = -1;
03053     }
03054     else if (paramType == paraCharmm)
03055     {
03056       // In a CHARMM psf file each dihedral will be only listed once
03057       // even if it has multiple terms. There is no point in comparing
03058       // to the psf information
03059       dihedral_array[index].multiplicity = ptr->multiplicity;
03060     } 
03061     //****** END CHARMM/XPLOR type changes
03062 
03063     for (i=0; i<ptr->multiplicity; i++)
03064     {
03065       dihedral_array[index].values[i].k = ptr->values[i].k;
03066       dihedral_array[index].values[i].n = ptr->values[i].n;
03067 
03068       //  Convert the angle to radians before storing it
03069       dihedral_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
03070     }
03071 
03072     ptr->index = index;
03073 
03074     index++;
03075     ptr=ptr->next;
03076   }
03077 }
03078 /*      END OF FUNCTION index_dihedrals      */
03079 
03080 /************************************************************************/
03081 /*                  */
03082 /*      FUNCTION index_impropers      */
03083 /*                  */
03084 /*  This function walks down the linked list of improper parameters */
03085 /*  and assigns an index to each one.  It also copies the data from this*/
03086 /*  linked list to the arrays that will be used from here on out  */
03087 /*                  */
03088 /************************************************************************/
03089 
03090 void Parameters::index_impropers()
03091 
03092 {
03093   struct improper_params *ptr;  //  Current place in list
03094   Index index=0;      //  Current index value
03095   int i;        //  Loop counter
03096 
03097   //  Allocate an array to hold the multiplicity present in the
03098   //  parameter file for each bond.  This will be used to check
03099   //  the multiplicities that are detected in the psf file
03100 
03101   //  This is kind of ugly, but necessary because of the way that
03102   //  X-PLOR psf files deal with Charmm22 parameters.  The way
03103   //  that multiple periodicities are specified is by having
03104   //  the bonds appear multiple times in the psf file.  This even
03105   //  if a bond type has multiple parameters defined, they
03106   //  will be used if the bond appears multiple times in the
03107   //  psf file.  So we need to store the number of parameters
03108   //  we have to make sure the psf file doesn't ask for more
03109   //  parameters than we really have, and we also need to track
03110   //  how many times the bond appears in the psf file so that
03111   //  we can decide how many parameters to actually use.
03112   maxImproperMults = new int[NumImproperParams];
03113 
03114   if (maxImproperMults == NULL)
03115   {
03116     NAMD_die("memory allocation failed in Parameters::index_impropers()");
03117   }
03118   
03119   //  Start at the head
03120   ptr = improperp;
03121 
03122   while (ptr != NULL)
03123   {
03124     //  Copy data to array and assign index
03125 
03126     //  Save the multiplicity in another array
03127     maxImproperMults[index] = ptr->multiplicity;
03128 
03129     //  Assign the multiplicity in the actual structure a bogus value
03130     //  that we will update in assign_dihedral_index
03131     improper_array[index].multiplicity = -1;
03132 
03133     for (i=0; i<ptr->multiplicity; i++)
03134     {
03135       improper_array[index].values[i].k = ptr->values[i].k;
03136       improper_array[index].values[i].n = ptr->values[i].n;
03137 
03138       //  Convert the angle to radians before storing it
03139       improper_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
03140     }
03141 
03142     ptr->index=index;
03143 
03144     index++;
03145     ptr=ptr->next;
03146   }
03147 }
03148 /*      END OF FUNCTION index_impropers      */
03149 
03150 
03151 /************************************************************************/
03152 /*                  */
03153 /*      FUNCTION index_crossterms      */
03154 /*                  */
03155 /*  This function walks down the linked list of crossterm parameters */
03156 /*  and assigns an index to each one.  It also copies the data from this*/
03157 /*  linked list to the arrays that will be used from here on out  */
03158 /*                  */
03159 /************************************************************************/
03160 
03161 void crossterm_setup(CrosstermData *);
03162 
03163 void Parameters::index_crossterms()
03164 
03165 {
03166   struct crossterm_params *ptr;  //  Current place in list
03167   Index index=0;      //  Current index value
03168   int i,j,k;        //  Loop counter
03169 
03170   //  Start at the head
03171   ptr = crosstermp;
03172 
03173   while (ptr != NULL)
03174   {
03175     //  Copy data to array and assign index
03176 
03177     int N = CrosstermValue::dim - 1;
03178 
03179     if ( ptr->dimension != N ) {
03180       NAMD_die("Sorry, only CMAP dimension of 24 is supported");
03181     }
03182 
03183     k = 0;
03184     for (i=0; i<N; i++) {
03185       for (j=0; j<N; j++) {
03186         crossterm_array[index].c[i][j].d00 = ptr->values[k];
03187         ++k;
03188       }
03189     }
03190     for (i=0; i<N; i++) {
03191         crossterm_array[index].c[i][N].d00 = 
03192                                 crossterm_array[index].c[i][0].d00;
03193         crossterm_array[index].c[N][i].d00 = 
03194                                 crossterm_array[index].c[0][i].d00;
03195     }
03196     crossterm_array[index].c[N][N].d00 = 
03197                                 crossterm_array[index].c[0][0].d00;
03198 
03199     crossterm_setup(&crossterm_array[index].c[0][0]);
03200 
03201     ptr->index=index;
03202 
03203     index++;
03204     ptr=ptr->next;
03205   }
03206 }
03207 /*      END OF FUNCTION index_crossterms      */
03208 
03209 /************************************************************************/
03210 /*                  */
03211 /*      FUNCTION index_vdw        */
03212 /*                  */
03213 /*   INPUTS:                */
03214 /*  tree - The tree that is to be indexed        */
03215 /*  index - index to start with          */
03216 /*                  */
03217 /*  This is a recursive routine that will traverse the binary tree  */
03218 /*   of vdw parameters, assigning an index to each one, and copying     */
03219 /*   the data from the binary tree to the array that will be used from  */
03220 /*   here on.                */
03221 /*                  */
03222 /************************************************************************/
03223 
03224 Index Parameters::index_vdw(struct vdw_params *tree, Index index)
03225 
03226 {
03227   //  If the tree is empty, do nothing
03228   if (tree==NULL)
03229     return(index);
03230 
03231   //  If I have a left subtree, populate it first
03232   if (tree->left != NULL)
03233   {
03234     index=index_vdw(tree->left, index);
03235   }
03236 
03237   //  Assign the index and copy the data to the array
03238   tree->index = index;
03239 
03240   vdw_array[index].sigma = tree->sigma;
03241   vdw_array[index].epsilon = tree->epsilon;
03242   vdw_array[index].sigma14 = tree->sigma14;
03243   vdw_array[index].epsilon14 = tree->epsilon14;
03244 
03245   char *nameloc = atom_type_name(index);
03246   strncpy(nameloc, tree->atomname, MAX_ATOMTYPE_CHARS);
03247   nameloc[MAX_ATOMTYPE_CHARS] = '\0';
03248 
03249 //  iout << iWARN << "Parameters: Stored name for type " << index << ": '";
03250 //      iout << iWARN << nameloc << "'" << "\n" << endi;
03251 
03252   index++;
03253 
03254   //  If I have a right subtree, index it
03255   if (tree->right != NULL)
03256   {
03257     index=index_vdw(tree->right, index);
03258   }
03259 
03260   return(index);
03261 }
03262 /*      END OF FUNCTION index_vdw      */
03263 
03264 /************************************************************************/
03265 /*                  */
03266 /*      FUNCTION assign_vdw_index      */
03267 /*                  */
03268 /*   INPUTS:                */
03269 /*  atomtype - atom type to find          */
03270 /*  atom_ptr - pointer to the atom structure to find vdw paramters  */
03271 /*       for              */
03272 /*                  */
03273 /*   OUTPUTS:                */
03274 /*  the vdw_index field of the atom structure is populated    */
03275 /*                  */
03276 /*  This function searches the binary tree of vdw parameters so     */
03277 /*   that an index can be assigned to this atom.  If the parameter is   */
03278 /*   is found, then the index is assigned.  If the parameter is not     */
03279 /*   found, then NAMD terminates.          */
03280 /*                  */
03281 /************************************************************************/
03282 #ifdef MEM_OPT_VERSION
03283 void Parameters::assign_vdw_index(char *atomtype, AtomCstInfo *atom_ptr)
03284 #else    
03285 void Parameters::assign_vdw_index(char *atomtype, Atom *atom_ptr)
03286 #endif
03287 {
03288   struct vdw_params *ptr;    //  Current position in trees
03289   int found=0;      //  Flag 1->found match
03290   int comp_code;      //  return code from strcasecmp
03291 
03292   /*  Check to make sure the files have all been read    */
03293   if (!AllFilesRead)
03294   {
03295     NAMD_die("Tried to assign vdw index before all parameter files were read");
03296   }
03297 
03298   /*  Start at the top            */
03299   ptr=vdwp;
03300 
03301   /*  While we haven't found a match, and we haven't reached      */
03302   /*  the bottom of the tree, compare the atom passed in with     */
03303   /*  the current value and decide if we have a match, or if not, */
03304   /*  which way to go            */
03305   while (!found && (ptr!=NULL))
03306   {
03307     comp_code = strcasecmp(atomtype, ptr->atomname);
03308 
03309     if (comp_code == 0)
03310     {
03311       /*  Found a match!        */
03312       atom_ptr->vdw_type=ptr->index;
03313       found=1;
03314     }
03315     else if (comp_code < 0)
03316     {
03317       /*  Go to the left        */
03318       ptr=ptr->left;
03319     }
03320     else
03321     {
03322       /*  Go to the right        */
03323       ptr=ptr->right;
03324     }
03325   }
03326 
03327   //****** BEGIN CHARMM/XPLOR type changes
03328   if (!found)
03329   {
03330     // since CHARMM allows wildcards "*" in vdw typenames
03331     // we have to look again if necessary, this way, if
03332     // we already had an exact match, this is never executed
03333     size_t windx;                      //  wildcard index
03334 
03335     /*  Start again at the top                                */
03336     ptr=vdwp;
03337   
03338      while (!found && (ptr!=NULL))
03339      {
03340   
03341        // get index of wildcard wildcard, get index
03342        windx= strcspn(ptr->atomname,"*"); 
03343        if (windx == strlen(ptr->atomname))
03344        {
03345          // there is no wildcard here
03346          comp_code = strcasecmp(atomtype, ptr->atomname);   
03347        }
03348        else
03349        {
03350          comp_code = strncasecmp(atomtype, ptr->atomname, windx); 
03351        }  
03352 
03353        if (comp_code == 0)
03354        {
03355          /*  Found a match!                                */
03356          atom_ptr->vdw_type=ptr->index;
03357          found=1;
03358          char errbuf[100];
03359          sprintf(errbuf,"VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
03360                         atomtype, ptr->atomname);
03361          int i;
03362          for(i=0; i<error_msgs.size(); i++) {
03363            if ( strcmp(errbuf,error_msgs[i]) == 0 ) break;
03364          }
03365          if ( i == error_msgs.size() ) {
03366            char *newbuf = new char[strlen(errbuf)+1];
03367            strcpy(newbuf,errbuf);
03368            error_msgs.add(newbuf);
03369            iout << iWARN << newbuf << "\n" << endi;
03370          }
03371        }
03372        else if (comp_code < 0)
03373        {
03374           /*  Go to the left                                */
03375                 ptr=ptr->left;
03376        }
03377        else
03378        {
03379          /*  Go to the right                                */
03380                 ptr=ptr->right;
03381        }
03382      
03383      }
03384                 
03385   }
03386   //****** END CHARMM/XPLOR type changes
03387 
03388   /*  Make sure we found it          */
03389   if (!found)
03390   {
03391     char err_msg[100];
03392 
03393     sprintf(err_msg, "DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s",
03394        atomtype);
03395     NAMD_die(err_msg);
03396   }
03397 
03398   return;
03399 }
03400 /*      END OF FUNCTION assign_vdw_index    */
03401 
03402 /************************************************************************
03403  * FUNCTION get_table_pair_params
03404  *
03405  * Inputs:
03406  * atom1 - atom type for atom 1
03407  * atom2 - atom type for atom 2
03408  * K - an integer value for the table type to populate
03409  *
03410  * Outputs:
03411  * If a match is found, K is populated and 1 is returned. Otherwise,
03412  * 0 is returned.
03413  *
03414  * This function finds the proper type index for tabulated nonbonded 
03415  * interactions between two atoms. If no such interactions are found,
03416  * the atoms are assumed to interact through standard VDW potentials.
03417  * 
03418  ************************************************************************/
03419 
03420 int Parameters::get_table_pair_params(Index ind1, Index ind2, int *K) {
03421   IndexedTablePair *ptr;
03422   Index temp;
03423   int found=FALSE;
03424 
03425   ptr=tab_pair_tree;
03426   //
03427   //  We need the smaller type in ind1, so if it isn't already that 
03428   //  way, switch them        */
03429   if (ind1 > ind2)
03430   {
03431     temp = ind1;
03432     ind1 = ind2;
03433     ind2 = temp;
03434   }
03435 
03436   /*  While we haven't found a match and we're not at the end  */
03437   /*  of the tree, compare the bond passed in with the tree  */
03438   while (!found && (ptr!=NULL))
03439   {
03440 //    printf("Comparing %i with %i and %i with %i\n", ind1, ptr->ind1, ind2, ptr->ind2);
03441     if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
03442     {
03443        found = TRUE;
03444     }
03445     else if ( (ind1 < ptr->ind1) || 
03446         ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
03447     {
03448       /*  Go left          */
03449       ptr=ptr->left;
03450     }
03451     else
03452     {
03453       /*  Go right          */
03454       ptr=ptr->right;
03455     }
03456   }
03457 
03458   /*  If we found a match, assign the values      */
03459   if (found)
03460   {
03461     *K = ptr->K;
03462     return(TRUE);
03463   }
03464   else
03465   {
03466     return(FALSE);
03467   }
03468 }
03469 /*      END OF FUNCTION get_table_pair_params    */
03470 
03471 /************************************************************************/
03472 /*                  */
03473 /*      FUNCTION get_vdw_pair_params      */
03474 /*                  */
03475 /*   INPUTS:                */
03476 /*  atom1 - atom type for atom 1          */
03477 /*  atom2 - atom type for atom 2          */
03478 /*  A - A value to populate            */
03479 /*  B - B value to populate            */
03480 /*  A14 - A 1-4 value to populate          */
03481 /*  B14 - B 1-4 value to populate          */
03482 /*                  */
03483 /*   OUTPUTS:                */
03484 /*  If a match is found, A, B, A14, and B14 are all populated and a */
03485 /*   1 is returned.  Otherwise, a 0 is returned.      */
03486 /*                    */
03487 /*  This function finds a set of vdw_pair paramters.  It is given   */
03488 /*   the two types of atoms involved.  This is the only paramter for    */
03489 /*   which a match is NOT guaranteed.  There will only be a match if    */
03490 /*   there are specific van der waals parameters for the two atom types */
03491 /*   involved.                */
03492 /*                  */
03493 /************************************************************************/
03494 
03495 int Parameters::get_vdw_pair_params(Index ind1, Index ind2, Real *A, 
03496         Real *B, Real *A14, Real *B14)
03497 
03498 {
03499   IndexedVdwPair *ptr;    //  Current location in tree
03500   Index temp;      //  Temporary value for swithcing
03501           // values
03502   int found=FALSE;    //  Flag 1-> found a match
03503 
03504   ptr=vdw_pair_tree;
03505 
03506   //  We need the smaller type in ind1, so if it isn't already that 
03507   //  way, switch them        */
03508   if (ind1 > ind2)
03509   {
03510     temp = ind1;
03511     ind1 = ind2;
03512     ind2 = temp;
03513   }
03514 
03515   /*  While we haven't found a match and we're not at the end  */
03516   /*  of the tree, compare the bond passed in with the tree  */
03517   while (!found && (ptr!=NULL))
03518   {
03519     if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
03520     {
03521        found = TRUE;
03522     }
03523     else if ( (ind1 < ptr->ind1) || 
03524         ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
03525     {
03526       /*  Go left          */
03527       ptr=ptr->left;
03528     }
03529     else
03530     {
03531       /*  Go right          */
03532       ptr=ptr->right;
03533     }
03534   }
03535 
03536   /*  If we found a match, assign the values      */
03537   if (found)
03538   {
03539     *A = ptr->A;
03540     *B = ptr->B;
03541     *A14 = ptr->A14;
03542     *B14 = ptr->B14;
03543 
03544     return(TRUE);
03545   }
03546   else
03547   {
03548     return(FALSE);
03549   }
03550 }
03551 /*      END OF FUNCTION get_vdw_pair_params    */
03552 
03553 
03554 /************************************************************************/
03555 /*                  */
03556 /*        FUNCTION assign_bond_index    */
03557 /*                  */
03558 /*   INPUTS:                */
03559 /*  atom1 - atom type for atom 1          */
03560 /*  atom2 - atom type for atom 2          */
03561 /*  bond_ptr - pointer to bond structure to populate    */
03562 /*                  */
03563 /*   OUTPUTS:                */
03564 /*  the structure pointed to by bond_ptr is populated    */
03565 /*                  */
03566 /*  This function finds a bond in the binary tree of bond values    */
03567 /*   and assigns its index.  If the bond is found, than the bond_type   */
03568 /*   field of the bond structure is populated.  If the parameter is     */
03569 /*   not found, NAMD will terminate.          */
03570 /*                  */
03571 /************************************************************************/
03572 
03573 void Parameters::assign_bond_index(char *atom1, char *atom2, Bond *bond_ptr)
03574 
03575 {
03576   struct bond_params *ptr;  //  Current location in tree
03577   int found=0;      //  Flag 1-> found a match
03578   int cmp_code;      //  return code from strcasecmp
03579   char tmp_name[15];    //  Temporary atom name
03580 
03581   /*  Check to make sure the files have all been read    */
03582   if (!AllFilesRead)
03583   {
03584     NAMD_die("Tried to assign bond index before all parameter files were read");
03585   }
03586 
03587   /*  We need atom1 < atom2, so if that's not the way they        */
03588   /*  were passed, flip them          */
03589   if (strcasecmp(atom1, atom2) > 0)
03590   {
03591     strcpy(tmp_name, atom1);
03592     strcpy(atom1, atom2);
03593     strcpy(atom2, tmp_name);
03594   }
03595 
03596   /*  Start at the top            */
03597   ptr=bondp;
03598 
03599   /*  While we haven't found a match and we're not at the end  */
03600   /*  of the tree, compare the bond passed in with the tree  */
03601   while (!found && (ptr!=NULL))
03602   {
03603     cmp_code=strcasecmp(atom1, ptr->atom1name);
03604 
03605     if (cmp_code == 0)
03606     {
03607       cmp_code=strcasecmp(atom2, ptr->atom2name);
03608     }
03609 
03610     if (cmp_code == 0)
03611     {
03612       /*  Found a match        */
03613       found=1;
03614       bond_ptr->bond_type = ptr->index;
03615     }
03616     else if (cmp_code < 0)
03617     {
03618       /*  Go left          */
03619       ptr=ptr->left;
03620     }
03621     else
03622     {
03623       /*  Go right          */
03624       ptr=ptr->right;
03625     }
03626   }
03627 
03628   /*  Check to see if we found anything        */
03629   if (!found)
03630   {
03631     char err_msg[128];
03632 
03633     sprintf(err_msg, "CAN'T FIND BOND PARAMETERS FOR BOND %s - %s IN PARAMETER FILES", atom1, atom2);
03634     NAMD_die(err_msg);
03635   }
03636 
03637   return;
03638 }
03639 /*      END OF FUNCTION assign_bond_index    */
03640 
03641 /************************************************************************/
03642 /*                  */
03643 /*      FUNCTION assign_angle_index      */
03644 /*                  */
03645 /*   INPUTS:                */
03646 /*  atom1 - atom type for atom 1          */
03647 /*  atom2 - atom type for atom 2          */
03648 /*  atom3 - atom type for atom 3          */
03649 /*  angle_ptr - pointer to angle structure to populate    */
03650 /*                  */
03651 /*   OUTPUTS:                */
03652 /*  the structure pointed to by angle_ptr is populated    */
03653 /*                  */
03654 /*  This function assigns an angle index to a specific angle.  */
03655 /*   It searches the binary tree of angle parameters for the appropriate*/
03656 /*   values.  If they are found, the index is assigned.  If they are    */
03657 /*   not found, then NAMD will terminate.        */
03658 /*                  */
03659 /************************************************************************/
03660 
03661 void Parameters::assign_angle_index(char *atom1, char *atom2, char*atom3,
03662           Angle *angle_ptr)
03663 
03664 {
03665   struct angle_params *ptr;  //  Current position in tree
03666   int comp_val;      //  value from strcasecmp
03667   int found=0;      //  flag 1->found a match
03668   char tmp_name[15];    //  Temporary atom name
03669 
03670   /*  Check to make sure the files have all been read    */
03671   if (!AllFilesRead)
03672   {
03673     NAMD_die("Tried to assign angle index before all parameter files were read");
03674   }
03675 
03676   /*  We need atom1 < atom3.  If that was not what we were   */
03677   /*  passed, switch them            */
03678   if (strcasecmp(atom1, atom3) > 0)
03679   {
03680     strcpy(tmp_name, atom1);
03681     strcpy(atom1, atom3);
03682     strcpy(atom3, tmp_name);
03683   }
03684 
03685   /*  Start at the top            */
03686   ptr=anglep;
03687 
03688   /*  While we don't have a match and we haven't reached the  */
03689   /*  bottom of the tree, compare values        */
03690   while (!found && (ptr != NULL))
03691   {
03692     comp_val = strcasecmp(atom1, ptr->atom1name);
03693 
03694     if (comp_val == 0)
03695     {
03696       /*  Atom 1 matches, so compare atom 2    */
03697       comp_val = strcasecmp(atom2, ptr->atom2name);
03698       
03699       if (comp_val == 0)
03700       {
03701         /*  Atoms 1&2 match, try atom 3    */
03702         comp_val = strcasecmp(atom3, ptr->atom3name);
03703       }
03704     }
03705 
03706     if (comp_val == 0)
03707     {
03708       /*  Found a match        */
03709       found = 1;
03710       angle_ptr->angle_type = ptr->index;
03711     }
03712     else if (comp_val < 0)
03713     {
03714       /*  Go left          */
03715       ptr=ptr->left;
03716     }
03717     else
03718     {
03719       /*  Go right          */
03720       ptr=ptr->right;
03721     }
03722   }
03723 
03724   /*  Make sure we found a match          */
03725   if (!found)
03726   {
03727     char err_msg[128];
03728 
03729     sprintf(err_msg, "UNABLE TO FIND ANGLE PARAMETERS FOR %s %s %s",
03730        atom1, atom2, atom3);
03731     NAMD_die(err_msg);
03732   }
03733 
03734   return;
03735 }
03736 /*      END OF FUNCTION assign_angle_index    */
03737 
03738 /************************************************************************/
03739 /*                  */
03740 /*      FUNCTION assign_dihedral_index      */
03741 /*                  */
03742 /*   INPUTS:                */
03743 /*  atom1 - atom type for atom 1          */
03744 /*  atom2 - atom type for atom 2          */
03745 /*  atom3 - atom type for atom 3          */
03746 /*  atom4 - atom type for atom 4          */
03747 /*  dihedral_ptr - pointer to dihedral structure to populate  */
03748 /*  multiplicity - Multiplicity to assign to this bond    */
03749 /*                  */
03750 /*   OUTPUTS:                */
03751 /*  the structure pointed to by dihedral_ptr is populated    */
03752 /*                  */
03753 /*  This function searchs the linked list of dihedral parameters for*/
03754 /*   a given bond.  If a match is found, the dihedral type is assigned. */
03755 /*   If no match is found, NAMD terminates        */
03756 /*                  */
03757 /************************************************************************/
03758 
03759 void Parameters::assign_dihedral_index(char *atom1, char *atom2, char *atom3,
03760         char *atom4, Dihedral *dihedral_ptr,
03761         int multiplicity)
03762 
03763 {
03764   struct dihedral_params *ptr;  //  Current position in list
03765   int found=0;      //  Flag 1->found a match
03766 
03767   /*  Start at the begining of the list        */
03768   ptr=dihedralp;
03769 
03770   /*  While we haven't found a match and we haven't reached       */
03771   /*  the end of the list, keep looking        */
03772   while (!found && (ptr!=NULL))
03773   {
03774     /*  Do a linear search through the linked list of   */
03775     /*  dihedral parameters.  Since the list is arranged    */
03776     /*  with wildcard paramters at the end of the list, we  */
03777     /*  can simply do a linear search and be guaranteed that*/
03778     /*  we will find exact matches before wildcard matches. */
03779     /*  Also, we must check for an exact match, and a match */
03780     /*  in reverse, since they are really the same          */
03781     /*  physically.            */
03782     if ( ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom1)==0) ) && 
03783          ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom2)==0) ) &&
03784          ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom3)==0) ) &&
03785          ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom4)==0) ) ) 
03786     {
03787       /*  Found an exact match      */
03788       found=1;
03789     }
03790     else if ( ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom1)==0) ) &&
03791               ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom2)==0) ) &&
03792               ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom3)==0) ) &&
03793               ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom4)==0) ) )
03794     {
03795       /*  Found a reverse match      */
03796       found=1;
03797     }
03798     else
03799     {
03800       /*  Didn't find a match, go to the next node  */
03801       ptr=ptr->next;
03802     }
03803   }
03804 
03805   /*  Make sure we found a match          */
03806   if (!found)
03807   {
03808     char err_msg[128];
03809 
03810     sprintf(err_msg, "CAN'T FIND DIHEDRAL PARAMETERS FOR %s  %s  %s  %s",
03811        atom1, atom2, atom3, atom4);
03812     
03813     NAMD_die(err_msg);
03814   }
03815 
03816   //  Check to make sure the number of multiples specified in the psf
03817   //  file doesn't exceed the number of parameters in the parameter
03818   //  files
03819   if (multiplicity > maxDihedralMults[ptr->index])
03820   {
03821     char err_msg[257];
03822 
03823     sprintf(err_msg, "Multiplicity of Paramters for diehedral bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxDihedralMults[ptr->index]);
03824     NAMD_die(err_msg);
03825   }
03826 
03827   //  If the multiplicity from the current bond is larger than that
03828   //  seen in the past, increase the multiplicity for this bond
03829   if (multiplicity > dihedral_array[ptr->index].multiplicity)
03830   {
03831     dihedral_array[ptr->index].multiplicity = multiplicity;
03832   }
03833 
03834   dihedral_ptr->dihedral_type = ptr->index;
03835 
03836   return;
03837 }
03838 /*      END OF FUNCTION assign_dihedral_index    */
03839 
03840 /************************************************************************/
03841 /*                  */
03842 /*      FUNCTION assign_improper_index      */
03843 /*                  */
03844 /*   INPUTS:                */
03845 /*  atom1 - atom type for atom 1          */
03846 /*  atom2 - atom type for atom 2          */
03847 /*  atom3 - atom type for atom 3          */
03848 /*  atom4 - atom type for atom 4          */
03849 /*  improper_ptr - pointer to improper structure to populate  */
03850 /*   multiplicity - Multiplicity to assign to this bond    */
03851 /*                  */
03852 /*   OUTPUTS:                */
03853 /*  the structure pointed to by improper_ptr is populated    */
03854 /*                  */
03855 /*  This function searchs the linked list of improper parameters for*/
03856 /*   a given bond.  If a match is found, the improper_type is assigned. */
03857 /*   If no match is found, NAMD will terminate.        */
03858 /*                  */
03859 /************************************************************************/
03860 
03861 void Parameters::assign_improper_index(char *atom1, char *atom2, char *atom3,
03862         char *atom4, Improper *improper_ptr,
03863         int multiplicity)
03864 
03865 {
03866   struct improper_params *ptr;  //  Current position in list
03867   int found=0;      //  Flag 1->found a match
03868 
03869   /*  Start at the head of the list        */
03870   ptr=improperp;
03871 
03872   /*  While we haven't fuond a match and haven't reached the end  */
03873   /*  of the list, keep looking          */
03874   while (!found && (ptr!=NULL))
03875   {
03876     /*  Do a linear search through the linked list of   */
03877     /*  improper parameters.  Since the list is arranged    */
03878     /*  with wildcard paramters at the end of the list, we  */
03879     /*  can simply do a linear search and be guaranteed that*/
03880     /*  we will find exact matches before wildcard matches. */
03881     /*  Also, we must check for an exact match, and a match */
03882     /*  in reverse, since they are really the same          */
03883     /*  physically.            */
03884     if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
03885            (strcasecmp(ptr->atom1name, "X")==0) ) &&
03886        ( (strcasecmp(ptr->atom2name, atom2)==0) || 
03887            (strcasecmp(ptr->atom2name, "X")==0) ) &&
03888        ( (strcasecmp(ptr->atom3name, atom3)==0) || 
03889            (strcasecmp(ptr->atom3name, "X")==0) ) &&
03890        ( (strcasecmp(ptr->atom4name, atom4)==0) || 
03891            (strcasecmp(ptr->atom4name, "X")==0) ) )
03892     {
03893       /*  Found an exact match      */
03894       found=1;
03895     }
03896     else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
03897            (strcasecmp(ptr->atom4name, "X")==0) ) &&
03898        ( (strcasecmp(ptr->atom3name, atom2)==0) || 
03899            (strcasecmp(ptr->atom3name, "X")==0) ) &&
03900        ( (strcasecmp(ptr->atom2name, atom3)==0) || 
03901            (strcasecmp(ptr->atom2name, "X")==0) ) &&
03902        ( (strcasecmp(ptr->atom1name, atom4)==0) || 
03903            (strcasecmp(ptr->atom1name, "X")==0) ) )
03904     {
03905       /*  Found a reverse match      */
03906       found=1;
03907     }
03908     else
03909     {
03910       /*  Didn't find a match, go to the next node  */
03911       ptr=ptr->next;
03912     }
03913   }
03914 
03915   /*  Make sure we found a match          */
03916   if (!found)
03917   {
03918     char err_msg[128];
03919 
03920     sprintf(err_msg, "CAN'T FIND IMPROPER PARAMETERS FOR %s  %s  %s  %s",
03921        atom1, atom2, atom3, atom4);
03922     
03923     NAMD_die(err_msg);
03924   }
03925 
03926   //  Check to make sure the number of multiples specified in the psf
03927   //  file doesn't exceed the number of parameters in the parameter
03928   //  files
03929   if (multiplicity > maxImproperMults[ptr->index])
03930   {
03931     char err_msg[257];
03932 
03933     sprintf(err_msg, "Multiplicity of Paramters for improper bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxImproperMults[ptr->index]);
03934     NAMD_die(err_msg);
03935   }
03936 
03937   //  If the multiplicity from the current bond is larger than that
03938   //  seen in the past, increase the multiplicity for this bond
03939   if (multiplicity > improper_array[ptr->index].multiplicity)
03940   {
03941     improper_array[ptr->index].multiplicity = multiplicity;
03942   }
03943 
03944   /*  Assign the constants          */
03945   improper_ptr->improper_type = ptr->index;
03946 
03947   return;
03948 }
03949 /*      END OF FUNCTION assign_improper_index    */
03950 
03951 /************************************************************************/
03952 /*                  */
03953 /*      FUNCTION assign_crossterm_index      */
03954 /*                  */
03955 /************************************************************************/
03956 
03957 void Parameters::assign_crossterm_index(char *atom1, char *atom2, char *atom3,
03958         char *atom4, char *atom5, char *atom6, char *atom7,
03959         char *atom8, Crossterm *crossterm_ptr)
03960 {
03961   struct crossterm_params *ptr;  //  Current position in list
03962   int found=0;      //  Flag 1->found a match
03963 
03964   /*  Start at the head of the list        */
03965   ptr=crosstermp;
03966 
03967   /*  While we haven't fuond a match and haven't reached the end  */
03968   /*  of the list, keep looking          */
03969   while (!found && (ptr!=NULL))
03970   {
03971     /*  Do a linear search through the linked list of   */
03972     /*  crossterm parameters.  Since the list is arranged    */
03973     /*  with wildcard paramters at the end of the list, we  */
03974     /*  can simply do a linear search and be guaranteed that*/
03975     /*  we will find exact matches before wildcard matches. */
03976     /*  Also, we must check for an exact match, and a match */
03977     /*  in reverse, since they are really the same          */
03978     /*  physically.            */
03979     if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
03980            (strcasecmp(ptr->atom1name, "X")==0) ) &&
03981        ( (strcasecmp(ptr->atom2name, atom2)==0) || 
03982            (strcasecmp(ptr->atom2name, "X")==0) ) &&
03983        ( (strcasecmp(ptr->atom3name, atom3)==0) || 
03984            (strcasecmp(ptr->atom3name, "X")==0) ) &&
03985        ( (strcasecmp(ptr->atom4name, atom4)==0) || 
03986            (strcasecmp(ptr->atom4name, "X")==0) ) )
03987     {
03988       /*  Found an exact match      */
03989       found=1;
03990     }
03991     else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
03992            (strcasecmp(ptr->atom4name, "X")==0) ) &&
03993        ( (strcasecmp(ptr->atom3name, atom2)==0) || 
03994            (strcasecmp(ptr->atom3name, "X")==0) ) &&
03995        ( (strcasecmp(ptr->atom2name, atom3)==0) || 
03996            (strcasecmp(ptr->atom2name, "X")==0) ) &&
03997        ( (strcasecmp(ptr->atom1name, atom4)==0) || 
03998            (strcasecmp(ptr->atom1name, "X")==0) ) )
03999     {
04000       /*  Found a reverse match      */
04001       found=1;
04002     }
04003     if ( ! found ) {
04004       /*  Didn't find a match, go to the next node  */
04005       ptr=ptr->next;
04006       continue;
04007     }
04008     found = 0;
04009     if ( ( (strcasecmp(ptr->atom5name, atom5)==0) || 
04010            (strcasecmp(ptr->atom5name, "X")==0) ) &&
04011        ( (strcasecmp(ptr->atom6name, atom6)==0) || 
04012            (strcasecmp(ptr->atom6name, "X")==0) ) &&
04013        ( (strcasecmp(ptr->atom7name, atom7)==0) || 
04014            (strcasecmp(ptr->atom7name, "X")==0) ) &&
04015        ( (strcasecmp(ptr->atom8name, atom8)==0) || 
04016            (strcasecmp(ptr->atom8name, "X")==0) ) )
04017     {
04018       /*  Found an exact match      */
04019       found=1;
04020     }
04021     else if ( ( (strcasecmp(ptr->atom8name, atom5)==0) || 
04022            (strcasecmp(ptr->atom8name, "X")==0) ) &&
04023        ( (strcasecmp(ptr->atom7name, atom6)==0) || 
04024            (strcasecmp(ptr->atom7name, "X")==0) ) &&
04025        ( (strcasecmp(ptr->atom6name, atom7)==0) || 
04026            (strcasecmp(ptr->atom6name, "X")==0) ) &&
04027        ( (strcasecmp(ptr->atom5name, atom8)==0) || 
04028            (strcasecmp(ptr->atom5name, "X")==0) ) )
04029     {
04030       /*  Found a reverse match      */
04031       found=1;
04032     }
04033     if ( ! found ) {
04034       /*  Didn't find a match, go to the next node  */
04035       ptr=ptr->next;
04036     }
04037   }
04038 
04039   /*  Make sure we found a match          */
04040   if (!found)
04041   {
04042     char err_msg[128];
04043 
04044     sprintf(err_msg, "CAN'T FIND CROSSTERM PARAMETERS FOR %s  %s  %s  %s  %s  %s  %s  %s",
04045        atom1, atom2, atom3, atom4, atom5, atom6, atom7, atom8);
04046     
04047     NAMD_die(err_msg);
04048   }
04049 
04050   /*  Assign the constants          */
04051   crossterm_ptr->crossterm_type = ptr->index;
04052 
04053   return;
04054 }
04055 /*      END OF FUNCTION assign_improper_index    */
04056 
04057 /************************************************************************/
04058 /*                  */
04059 /*      FUNCTION free_bond_tree        */
04060 /*                  */
04061 /*   INPUTS:                */
04062 /*  bond_ptr - pointer to bond tree to free        */
04063 /*                  */
04064 /*  this is a recursive function that is used to free the memory    */
04065 /*   allocated for a bond paramter tree.  It makes recursive calls to   */
04066 /*   free the left an right subtress, and then frees the head.  It is   */
04067 /*   only called by the destructor          */
04068 /*                  */
04069 /************************************************************************/
04070 
04071 void Parameters::free_bond_tree(struct bond_params *bond_ptr)
04072 
04073 {
04074   if (bond_ptr->left != NULL)
04075   {
04076     free_bond_tree(bond_ptr->left);
04077   }
04078 
04079   if (bond_ptr->right != NULL)
04080   {
04081     free_bond_tree(bond_ptr->right);
04082   }
04083 
04084   delete bond_ptr;
04085 
04086   return;
04087 }
04088 /*      END OF FUNCTION free_bond_tree      */
04089 
04090 /************************************************************************/
04091 /*                  */
04092 /*      FUNCTION free_angle_tree      */
04093 /*                  */
04094 /*   INPUTS:                */
04095 /*  angle_ptr - pointer to angle tree to free      */
04096 /*                  */
04097 /*  this is a recursive function that is used to free the memory    */
04098 /*   allocated for a angle paramter tree.  It makes recursive calls to  */
04099 /*   free the left an right subtress, and then frees the head.  It is   */
04100 /*   only called by the destructor          */
04101 /*                  */
04102 /************************************************************************/
04103 
04104 void Parameters::free_angle_tree(struct angle_params *angle_ptr)
04105 
04106 {
04107   if (angle_ptr->left != NULL)
04108   {
04109     free_angle_tree(angle_ptr->left);
04110   }
04111 
04112   if (angle_ptr->right != NULL)
04113   {
04114     free_angle_tree(angle_ptr->right);
04115   }
04116 
04117   delete angle_ptr;
04118 
04119   return;
04120 }
04121 /*      END OF FUNCTION free_angle_tree      */
04122 
04123 /************************************************************************/
04124 /*                  */
04125 /*      FUNCTION free_dihedral_list      */
04126 /*                  */
04127 /*   INPUTS:                */
04128 /*  dih_ptr - pointer to the list to free        */
04129 /*                  */
04130 /*  this function frees a linked list of dihedral parameters.  It   */
04131 /*   is only called by the destructor.          */
04132 /*                  */
04133 /************************************************************************/
04134 
04135 void Parameters::free_dihedral_list(struct dihedral_params *dih_ptr)
04136 
04137 {
04138   struct dihedral_params *ptr;  //  Current position in list
04139   struct dihedral_params *next; //  Next position in list
04140 
04141   ptr=dih_ptr;
04142 
04143   while (ptr != NULL)
04144   {
04145     next=ptr->next;
04146     delete ptr;
04147     ptr=next;
04148   }
04149 
04150   return;
04151 }
04152 /*      END OF FUNCTION free_dihedral_list    */
04153 
04154 /************************************************************************/
04155 /*                  */
04156 /*      FUNCTION free_improper_list      */
04157 /*                  */
04158 /*   INPUTS:                */
04159 /*  imp_ptr - pointer to the list to free        */
04160 /*                  */
04161 /*  this function frees a linked list of improper parameters.  It   */
04162 /*   is only called by the destructor.          */
04163 /*                  */
04164 /************************************************************************/
04165 
04166 void Parameters::free_improper_list(struct improper_params *imp_ptr)
04167 
04168 {
04169   struct improper_params *ptr;  //  Current position in list
04170   struct improper_params *next; //  Next position in list
04171 
04172   ptr=imp_ptr;
04173 
04174   while (ptr != NULL)
04175   {
04176     next=ptr->next;
04177     delete ptr;
04178     ptr=next;
04179   }
04180 
04181   return;
04182 }
04183 /*      END OF FUNCTION free_improper_list    */
04184     
04185 /************************************************************************/
04186 /*                  */
04187 /*      FUNCTION free_crossterm_list      */
04188 /*                  */
04189 /*   INPUTS:                */
04190 /*  imp_ptr - pointer to the list to free        */
04191 /*                  */
04192 /*  this function frees a linked list of crossterm parameters.  It   */
04193 /*   is only called by the destructor.          */
04194 /*                  */
04195 /************************************************************************/
04196 
04197 void Parameters::free_crossterm_list(struct crossterm_params *imp_ptr)
04198 
04199 {
04200   struct crossterm_params *ptr;  //  Current position in list
04201   struct crossterm_params *next; //  Next position in list
04202 
04203   ptr=imp_ptr;
04204 
04205   while (ptr != NULL)
04206   {
04207     next=ptr->next;
04208     delete ptr;
04209     ptr=next;
04210   }
04211 
04212   return;
04213 }
04214 /*      END OF FUNCTION free_crossterm_list    */
04215     
04216 /************************************************************************/
04217 /*                  */
04218 /*      FUNCTION free_vdw_tree        */
04219 /*                  */
04220 /*   INPUTS:                */
04221 /*  vdw_ptr - pointer to vdw tree to free        */
04222 /*                  */
04223 /*  this is a recursive function that is used to free the memory    */
04224 /*   allocated for a vdw paramter tree.  It makes recursive calls to    */
04225 /*   free the left an right subtress, and then frees the head.  It is   */
04226 /*   only called by the destructor          */
04227 /*                  */
04228 /************************************************************************/
04229 
04230 void Parameters::free_vdw_tree(struct vdw_params *vdw_ptr)
04231 
04232 {
04233   if (vdw_ptr->left != NULL)
04234   {
04235     free_vdw_tree(vdw_ptr->left);
04236   }
04237 
04238   if (vdw_ptr->right != NULL)
04239   {
04240     free_vdw_tree(vdw_ptr->right);
04241   }
04242 
04243   delete vdw_ptr;
04244 
04245   return;
04246 }
04247 /*      END OF FUNCTION free_vdw_tree      */
04248 
04249 /************************************************************************/
04250 /*                  */
04251 /*      FUNCTION free_vdw_pair_list      */
04252 /*                  */
04253 /*  This function frees the vdw_pair_list        */
04254 /*                  */
04255 /************************************************************************/
04256 
04257 void Parameters::free_vdw_pair_list()
04258 {
04259    struct vdw_pair_params *ptr, *next;
04260    
04261    ptr=vdw_pairp;
04262    
04263    while (ptr != NULL)
04264    {
04265       next = ptr->next;
04266       
04267       delete ptr;
04268       
04269       ptr = next;
04270    }
04271    
04272    vdw_pairp = NULL;
04273 }
04274 /*      END OF FUNCTION free_vdw_pair_list    */
04275 
04276 /************************************************************************
04277  * FUNCTION free_table_pair_tree
04278  *
04279  * Free a table_pair_tree given a pointer to its head
04280  * **********************************************************************/
04281 
04282 void Parameters::free_table_pair_tree(IndexedTablePair *table_pair_ptr) {
04283   if (table_pair_ptr->left != NULL)
04284   {
04285     free_table_pair_tree(table_pair_ptr->left);
04286   }
04287 
04288   if (table_pair_ptr->right != NULL)
04289   {
04290     free_table_pair_tree(table_pair_ptr->right);
04291   }
04292 
04293   delete table_pair_ptr;
04294 
04295   return;
04296 }
04297 
04298 
04299 /************************************************************************/
04300 /*                  */
04301 /*      FUNCTION free_vdw_pair_tree      */
04302 /*                  */
04303 /*   INPUTS:                */
04304 /*  vdw_pair_ptr - pointer to vdw_pair tree to free      */
04305 /*                  */
04306 /*  this is a recursive function that is used to free the memory    */
04307 /*   allocated for a vdw_pair paramter tree.  It makes recursive calls  */
04308 /*   to free the left an right subtress, and then frees the head.  It is*/
04309 /*   only called by the destructor          */
04310 /*                  */
04311 /************************************************************************/
04312 
04313 void Parameters::free_vdw_pair_tree(IndexedVdwPair *vdw_pair_ptr)
04314 
04315 {
04316   if (vdw_pair_ptr->left != NULL)
04317   {
04318     free_vdw_pair_tree(vdw_pair_ptr->left);
04319   }
04320 
04321   if (vdw_pair_ptr->right != NULL)
04322   {
04323     free_vdw_pair_tree(vdw_pair_ptr->right);
04324   }
04325 
04326   delete vdw_pair_ptr;
04327 
04328   return;
04329 }
04330 /*      END OF FUNCTION free_vdw_pair_tree    */
04331 
04332 /************************************************************************/
04333 /*                  */
04334 /*      FUNCTION traverse_bond_params      */
04335 /*                  */
04336 /*   INPUTS:                */
04337 /*  tree - the bond binary tree to traverse        */
04338 /*                  */
04339 /*  This is a recursive call used for debugging purposes that  */
04340 /*   prints out all the bond paramters in the bond parameter binary  */
04341 /*   search tree. It is only called by print_bond_params    */
04342 /*                  */
04343 /************************************************************************/
04344 
04345 void Parameters::traverse_bond_params(struct bond_params *tree)
04346 
04347 {
04348   if (tree==NULL)
04349     return;
04350 
04351   if (tree->left != NULL)
04352   {
04353     traverse_bond_params(tree->left);
04354   }
04355 
04356   DebugM(3,"BOND " <<  tree->atom1name << "  " << tree->atom2name \
04357       << " index=" << tree->index << " k=" << tree->forceconstant \
04358       << " x0=" << tree->distance);
04359 
04360   if (tree->right != NULL)
04361   {
04362     traverse_bond_params(tree->right);
04363   }
04364 }
04365 /*      END OF FUNCTION traverse_bond_params    */
04366 
04367 /************************************************************************/
04368 /*                  */
04369 /*      FUNCTION traverse_angle_params      */
04370 /*                  */
04371 /*   INPUTS:                */
04372 /*  tree - the angle binary tree to traverse      */
04373 /*                  */
04374 /*  This is a recursive call used for debugging purposes that  */
04375 /*   prints out all the angle paramters in the angle parameter binary  */
04376 /*   search tree. It is only called by print_angle_params    */
04377 /*                  */
04378 /************************************************************************/
04379 
04380 void Parameters::traverse_angle_params(struct angle_params *tree)
04381 
04382 {
04383   if (tree==NULL)
04384     return;
04385 
04386   if (tree->left != NULL)
04387   {
04388     traverse_angle_params(tree->left);
04389   }
04390   DebugM(3,"ANGLE " << tree->atom1name << "  " << tree->atom2name \
04391       << "  " << tree->atom3name << " index=" << tree->index \
04392       << " k=" << tree->forceconstant << " theta0=" << tree->angle \
04393       );
04394 
04395   if (tree->right != NULL)
04396   {
04397     traverse_angle_params(tree->right);
04398   }
04399 }
04400 /*      END OF FUNCTION traverse_angle_params    */
04401 
04402 /************************************************************************/
04403 /*                  */
04404 /*      FUNCTION traverse_dihedral_params    */
04405 /*                  */
04406 /*   INPUTS:                */
04407 /*  list - the dihedral linked list to traverse      */
04408 /*                  */
04409 /*  This is a call used for debugging purposes that prints out all  */
04410 /*   the bond paramters in the dihedral parameter linked list. It is    */
04411 /*   only called by print_dihedral_params.        */
04412 /*                  */
04413 /************************************************************************/
04414 
04415 void Parameters::traverse_dihedral_params(struct dihedral_params *list)
04416 
04417 {
04418   int i;
04419 
04420   while (list != NULL)
04421   {
04422     DebugM(3,"DIHEDRAL  " << list->atom1name << "  " \
04423         << list->atom2name << "  " << list->atom3name \
04424         << "  " << list->atom4name << " index=" \
04425         << list->index \
04426         << " multiplicity=" << list->multiplicity << "\n");
04427         
04428     for (i=0; i<list->multiplicity; i++)
04429     {
04430       DebugM(3,"k=" << list->values[i].k \
04431           << " n=" << list->values[i].n  \
04432           << " delta=" << list->values[i].delta);
04433     }
04434 
04435     list=list->next;
04436   }
04437 }
04438 /*      END OF FUNCTION traverse_dihedral_params  */
04439 
04440 /************************************************************************/
04441 /*                  */
04442 /*      FUNCTION traverse_improper_params    */
04443 /*                  */
04444 /*   INPUTS:                */
04445 /*  list - the improper linked list to traverse      */
04446 /*                  */
04447 /*  This is a call used for debugging purposes that prints out all  */
04448 /*   the improper paramters in the improper parameter linked list. It is*/
04449 /*   only called by print_improper_params.        */
04450 /*                  */
04451 /************************************************************************/
04452 
04453 void Parameters::traverse_improper_params(struct improper_params *list)
04454 
04455 {
04456   int i;
04457 
04458   while (list != NULL)
04459   {
04460     DebugM(3,"Improper  " << list->atom1name << "  " \
04461         << list->atom2name << "  " << list->atom3name \
04462         << "  " << list->atom4name << " index="  \
04463         << list->index  \
04464         << " multiplicity=" << list->multiplicity << "\n");
04465 
04466     for (i=0; i<list->multiplicity; i++)
04467     {
04468        DebugM(3,"k=" << list->values[i].k \
04469            << " n=" << list->values[i].n \
04470            << " delta=" << list->values[i].delta);
04471     }
04472 
04473     list=list->next;
04474   }
04475 }
04476 /*      END OF FUNCTION traverse_improper_params  */
04477 
04478 
04479 /************************************************************************/
04480 /*                  */
04481 /*      FUNCTION traverse_vdw_params      */
04482 /*                  */
04483 /*   INPUTS:                */
04484 /*  tree - the vw binary tree to traverse        */
04485 /*                  */
04486 /*  This is a recursive call used for debugging purposes that  */
04487 /*   prints out all the vdw paramters in the vdw parameter binary  */
04488 /*   search tree. It is only called by print_vdw_params      */
04489 /*                  */
04490 /************************************************************************/
04491 
04492 void Parameters::traverse_vdw_params(struct vdw_params *tree)
04493 
04494 {
04495   if (tree==NULL)
04496     return;
04497 
04498   if (tree->left != NULL)
04499   {
04500     traverse_vdw_params(tree->left);
04501   }
04502 
04503   DebugM(3,"vdW " << tree->atomname << " index=" << tree->index \
04504       << " sigma=" << tree->sigma << " epsilon=" << \
04505       tree->epsilon << " sigma 1-4=" << tree->sigma14 \
04506       << " epsilon 1-4=" << tree->epsilon14 << std::endl);
04507 
04508   if (tree->right != NULL)
04509   {
04510     traverse_vdw_params(tree->right);
04511   }
04512 }
04513 /*      END OF FUNCTION traverse_vdw_params    */
04514 
04515 
04516 /************************************************************************/
04517 /*                  */
04518 /*      FUNCTION traverse_vdw_pair_params    */
04519 /*                  */
04520 /*   INPUTS:                */
04521 /*  list - the vdw_pair list to traverse        */
04522 /*                  */
04523 /*  This call simply prints out the vdw_pair list      */
04524 /*                  */
04525 /************************************************************************/
04526 
04527 void Parameters::traverse_vdw_pair_params(struct vdw_pair_params *list)
04528 
04529 {
04530   if (list==NULL)
04531     return;
04532 
04533   DebugM(3,"vdW PAIR  " << list->atom1name << "  "  \
04534       << list->atom2name << " A=" << list->A \
04535       << " B=" << list->B << " A 1-4=" \
04536       << list->A14 << " B 1-4=" << list->B14 \
04537       );
04538 
04539   traverse_vdw_pair_params(list->next);
04540 }
04541 /*      END OF FUNCTION traverse_vdw_pair_params  */
04542 
04543 /************************************************************************/
04544 /*                  */
04545 /*      FUNCTION print_bond_params      */
04546 /*                  */
04547 /*  This is a debugging routine used to print out all the bond  */
04548 /*  parameters                */
04549 /*                  */
04550 /************************************************************************/
04551 
04552 void Parameters::print_bond_params()
04553 {
04554   DebugM(3,NumBondParams << " BOND PARAMETERS\n" \
04555       << "*****************************************"  \
04556       );
04557 
04558   traverse_bond_params(bondp);
04559 }
04560 
04561 /************************************************************************/
04562 /*                  */
04563 /*      FUNCTION print_angle_params      */
04564 /*                  */
04565 /*  This is a debugging routine used to print out all the angle  */
04566 /*  parameters                */
04567 /*                  */
04568 /************************************************************************/
04569 
04570 void Parameters::print_angle_params()
04571 {
04572   DebugM(3,NumAngleParams << " ANGLE PARAMETERS\n"
04573       << "*****************************************" );
04574   traverse_angle_params(anglep);
04575 }
04576 
04577 /************************************************************************/
04578 /*                  */
04579 /*      FUNCTION print_dihedral_params      */
04580 /*                  */
04581 /*  This is a debugging routine used to print out all the dihedral  */
04582 /*  parameters                */
04583 /*                  */
04584 /************************************************************************/
04585 
04586 void Parameters::print_dihedral_params()
04587 {
04588   DebugM(3,NumDihedralParams << " DIHEDRAL PARAMETERS\n" \
04589       << "*****************************************" );
04590 
04591   traverse_dihedral_params(dihedralp);
04592 }
04593 
04594 /************************************************************************/
04595 /*                  */
04596 /*      FUNCTION print_improper_params      */
04597 /*                  */
04598 /*  This is a debugging routine used to print out all the improper  */
04599 /*  parameters                */
04600 /*                  */
04601 /************************************************************************/
04602 
04603 void Parameters::print_improper_params()
04604 {
04605   DebugM(3,NumImproperParams << " IMPROPER PARAMETERS\n" \
04606       << "*****************************************" );
04607 
04608   traverse_improper_params(improperp);
04609 }
04610 
04611 /************************************************************************/
04612 /*                  */
04613 /*      FUNCTION print_vdw_params      */
04614 /*                  */
04615 /*  This is a debugging routine used to print out all the vdw  */
04616 /*  parameters                */
04617 /*                  */
04618 /************************************************************************/
04619 
04620 void Parameters::print_vdw_params()
04621 {
04622   DebugM(3,NumVdwParams << " vdW PARAMETERS\n" \
04623       << "*****************************************" );
04624 
04625   traverse_vdw_params(vdwp);
04626 }
04627 
04628 /************************************************************************/
04629 /*                  */
04630 /*      FUNCTION print_vdw_pair_params      */
04631 /*                  */
04632 /*  This is a debugging routine used to print out all the vdw_pair  */
04633 /*  parameters                */
04634 /*                  */
04635 /************************************************************************/
04636 
04637 void Parameters::print_vdw_pair_params()
04638 {
04639   DebugM(3,NumVdwPairParams << " vdW PAIR PARAMETERS\n" \
04640       << "*****************************************" );
04641 
04642   traverse_vdw_pair_params(vdw_pairp);
04643 }
04644 
04645 /************************************************************************/
04646 /*                  */
04647 /*      FUNCTION print_param_summary      */
04648 /*                  */
04649 /*  This function just prints out a brief summary of the paramters  */
04650 /*  that have been read in.  It is intended for debugging purposes  */
04651 /*                  */
04652 /************************************************************************/
04653 
04654 void Parameters::print_param_summary()
04655 {
04656   iout << iINFO << "SUMMARY OF PARAMETERS:\n" 
04657        << iINFO << NumBondParams << " BONDS\n" 
04658        << iINFO << NumAngleParams << " ANGLES\n" << endi;
04659        if (cosAngles) {
04660          iout << iINFO << "  " << (NumAngleParams - NumCosAngles) << " HARMONIC\n"
04661               << iINFO << "  " << NumCosAngles << " COSINE-BASED\n" << endi;
04662        }
04663   iout << iINFO << NumDihedralParams << " DIHEDRAL\n"
04664        << iINFO << NumImproperParams << " IMPROPER\n"
04665        << iINFO << NumCrosstermParams << " CROSSTERM\n"
04666        << iINFO << NumVdwParams << " VDW\n"
04667        << iINFO << NumVdwPairParams << " VDW_PAIRS\n" << endi;
04668 }
04669 
04670 
04671 /************************************************************************/
04672 /*                  */
04673 /*      FUNCTION done_reading_structure      */
04674 /*                  */
04675 /*  This function is used to tell the Parameters object that the    */
04676 /*  structure has been read in.  This is so that the Parameters object  */
04677 /*  can now release the binary trees and linked lists that it was using */
04678 /*  to search for parameters based on the atom type.  From this point   */
04679 /*  on, only the arrays of parameter data will be used.  If this object */
04680 /*  resides on any node BUT the master node, it will never even have    */
04681 /*  these trees and lists.  For the master node, this just frees up     */
04682 /*  some memory for better uses.          */
04683 /*                  */
04684 /************************************************************************/
04685 
04686 void Parameters::done_reading_structure()
04687 
04688 {
04689   if (bondp != NULL)
04690     free_bond_tree(bondp);
04691 
04692   if (anglep != NULL)
04693     free_angle_tree(anglep);
04694 
04695   if (dihedralp != NULL)
04696     free_dihedral_list(dihedralp);
04697 
04698   if (improperp != NULL)
04699     free_improper_list(improperp);
04700 
04701   if (crosstermp != NULL)
04702     free_crossterm_list(crosstermp);
04703 
04704   if (vdwp != NULL)
04705     free_vdw_tree(vdwp);
04706 
04707   //  Free the arrays used to track multiplicity for dihedrals
04708   //  and impropers
04709   if (maxDihedralMults != NULL)
04710     delete [] maxDihedralMults;
04711 
04712   if (maxImproperMults != NULL)
04713     delete [] maxImproperMults;
04714 
04715   bondp=NULL;
04716   anglep=NULL;
04717   dihedralp=NULL;
04718   improperp=NULL;
04719   crosstermp=NULL;
04720   vdwp=NULL;
04721   maxImproperMults=NULL;
04722   maxDihedralMults=NULL;
04723 }
04724 /*      END OF FUNCTION done_reading_structure    */
04725 
04726 /************************************************************************/
04727 /*                  */
04728 /*      FUNCTION send_Parameters      */
04729 /*                  */
04730 /*  This function is used by the master node to broadcast the       */
04731 /*   structure Parameters to all the other nodes.        */
04732 /*                  */
04733 /************************************************************************/
04734 
04735 void Parameters::send_Parameters(MOStream *msg)
04736 {
04737   Real *a1, *a2, *a3, *a4;        //  Temporary arrays for sending messages
04738   int *i1, *i2, *i3;      //  Temporary int array
04739   int i, j;      //  Loop counters
04740   Real **kvals;      //  Force constant values for dihedrals and impropers
04741   int **nvals;      //  Periodicity values for  dihedrals and impropers
04742   Real **deltavals;    //  Phase shift values for  dihedrals and impropers
04743   /*MOStream *msg=comm_obj->newOutputStream(ALLBUTME, STATICPARAMSTAG, BUFSIZE);
04744   if ( msg == NULL )
04745   {
04746     NAMD_die("memory allocation failed in Parameters::send_Parameters");
04747   }*/
04748 
04749   //  Send the bond parameters
04750   msg->put(NumBondParams);
04751 
04752   if (NumBondParams)
04753   {
04754     a1 = new Real[NumBondParams];
04755     a2 = new Real[NumBondParams];
04756 
04757     if ( (a1 == NULL) || (a2 == NULL) )
04758     {
04759       NAMD_die("memory allocation failed in Parameters::send_Parameters");
04760     }
04761 
04762     for (i=0; i<NumBondParams; i++)
04763     {
04764       a1[i] = bond_array[i].k;
04765       a2[i] = bond_array[i].x0;
04766     }
04767 
04768     msg->put(NumBondParams, a1)->put(NumBondParams, a2);
04769 
04770     delete [] a1;
04771     delete [] a2;
04772   }
04773 
04774   //  Send the angle parameters
04775   msg->put(NumAngleParams);
04776 
04777   if (NumAngleParams)
04778   {
04779     a1 = new Real[NumAngleParams];
04780     a2 = new Real[NumAngleParams];
04781     a3 = new Real[NumAngleParams];
04782     a4 = new Real[NumAngleParams];
04783     i1 = new int[NumAngleParams];
04784 
04785     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) ||
04786          (a4 == NULL) || (i1 == NULL))
04787     {
04788       NAMD_die("memory allocation failed in Parameters::send_Parameters");
04789     }
04790 
04791     for (i=0; i<NumAngleParams; i++)
04792     {
04793       a1[i] = angle_array[i].k;
04794       a2[i] = angle_array[i].theta0;
04795       a3[i] = angle_array[i].k_ub;
04796       a4[i] = angle_array[i].r_ub;
04797       i1[i] = angle_array[i].normal;
04798     }
04799 
04800     msg->put(NumAngleParams, a1)->put(NumAngleParams, a2);
04801     msg->put(NumAngleParams, a3)->put(NumAngleParams, a4);
04802     msg->put(NumAngleParams, i1);
04803 
04804     delete [] a1;
04805     delete [] a2;
04806     delete [] a3;
04807     delete [] a4;
04808     delete [] i1;
04809   }
04810 
04811   //  Send the dihedral parameters
04812   msg->put(NumDihedralParams);
04813 
04814   if (NumDihedralParams)
04815   {
04816     i1 = new int[NumDihedralParams];
04817     kvals = new Real *[MAX_MULTIPLICITY];
04818     nvals = new int *[MAX_MULTIPLICITY];
04819     deltavals = new Real *[MAX_MULTIPLICITY];
04820 
04821     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
04822          (deltavals == NULL) )
04823     {
04824       NAMD_die("memory allocation failed in Parameters::send_Parameters");
04825     }
04826 
04827     for (i=0; i<MAX_MULTIPLICITY; i++)
04828     {
04829       kvals[i] = new Real[NumDihedralParams];
04830       nvals[i] = new int[NumDihedralParams];
04831       deltavals[i] = new Real[NumDihedralParams];
04832 
04833       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
04834       {
04835         NAMD_die("memory allocation failed in Parameters::send_Parameters");
04836       }
04837     }
04838 
04839     for (i=0; i<NumDihedralParams; i++)
04840     {
04841       i1[i] = dihedral_array[i].multiplicity;
04842 
04843       for (j=0; j<MAX_MULTIPLICITY; j++)
04844       {
04845         kvals[j][i] = dihedral_array[i].values[j].k;
04846         nvals[j][i] = dihedral_array[i].values[j].n;
04847         deltavals[j][i] = dihedral_array[i].values[j].delta;
04848       }
04849     }
04850 
04851     msg->put(NumDihedralParams, i1);
04852 
04853     for (i=0; i<MAX_MULTIPLICITY; i++)
04854     {
04855       msg->put(NumDihedralParams, kvals[i]);
04856       msg->put(NumDihedralParams, nvals[i]);
04857       msg->put(NumDihedralParams, deltavals[i]);
04858 
04859       delete [] kvals[i];
04860       delete [] nvals[i];
04861       delete [] deltavals[i];
04862     }
04863 
04864     delete [] i1;
04865     delete [] kvals;
04866     delete [] nvals;
04867     delete [] deltavals;
04868   }
04869 
04870   //  Send the improper parameters
04871   msg->put(NumImproperParams);
04872 
04873   if (NumImproperParams)
04874   {
04875     i1 = new int[NumImproperParams];
04876     kvals = new Real *[MAX_MULTIPLICITY];
04877     nvals = new int *[MAX_MULTIPLICITY];
04878     deltavals = new Real *[MAX_MULTIPLICITY];
04879 
04880     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
04881          (deltavals == NULL) )
04882     {
04883       NAMD_die("memory allocation failed in Parameters::send_Parameters");
04884     }
04885 
04886     for (i=0; i<MAX_MULTIPLICITY; i++)
04887     {
04888       kvals[i] = new Real[NumImproperParams];
04889       nvals[i] = new int[NumImproperParams];
04890       deltavals[i] = new Real[NumImproperParams];
04891 
04892       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
04893       {
04894         NAMD_die("memory allocation failed in Parameters::send_Parameters");
04895       }
04896     }
04897 
04898     for (i=0; i<NumImproperParams; i++)
04899     {
04900       i1[i] = improper_array[i].multiplicity;
04901 
04902       for (j=0; j<MAX_MULTIPLICITY; j++)
04903       {
04904         kvals[j][i] = improper_array[i].values[j].k;
04905         nvals[j][i] = improper_array[i].values[j].n;
04906         deltavals[j][i] = improper_array[i].values[j].delta;
04907       }
04908     }
04909 
04910     msg->put(NumImproperParams, i1);
04911 
04912     for (i=0; i<MAX_MULTIPLICITY; i++)
04913     {
04914       msg->put(NumImproperParams, kvals[i]);
04915       msg->put(NumImproperParams, nvals[i]);
04916       msg->put(NumImproperParams, deltavals[i]);
04917 
04918       delete [] kvals[i];
04919       delete [] nvals[i];
04920       delete [] deltavals[i];
04921     }
04922 
04923     delete [] i1;
04924     delete [] kvals;
04925     delete [] nvals;
04926     delete [] deltavals;
04927   }
04928 
04929   //  Send the crossterm parameters
04930   msg->put(NumCrosstermParams);
04931 
04932   if (NumCrosstermParams)
04933   {
04934     for (i=0; i<NumCrosstermParams; ++i) {
04935       int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
04936       msg->put(nvals,&crossterm_array[i].c[0][0].d00);
04937     }
04938   }
04939   //
04940   //Send the energy table parameters
04941   msg->put(numenerentries);
04942 
04943   if (numenerentries) {
04944           /*
04945     b1 = new Real[numenerentries];
04946     if (b1 == NULL) 
04947     {
04948       NAMD_die("memory allocation failed in Parameters::send_Parameters");
04949     }
04950     */
04951     
04952     msg->put(numenerentries, table_ener);
04953   }
04954 
04955   //  Send the vdw parameters
04956   msg->put(NumVdwParams);
04957   msg->put(NumVdwParamsAssigned);
04958 
04959   if (NumVdwParams)
04960   {
04961     a1 = new Real[NumVdwParams];
04962     a2 = new Real[NumVdwParams];
04963     a3 = new Real[NumVdwParams];
04964     a4 = new Real[NumVdwParams];
04965 
04966     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) )
04967     {
04968       NAMD_die("memory allocation failed in Parameters::send_Parameters");
04969     }
04970 
04971     for (i=0; i<NumVdwParams; i++)
04972     {
04973       a1[i] = vdw_array[i].sigma;
04974       a2[i] = vdw_array[i].epsilon;
04975       a3[i] = vdw_array[i].sigma14;
04976       a4[i] = vdw_array[i].epsilon14;
04977     }
04978 
04979     msg->put(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
04980     msg->put(NumVdwParams, a1);
04981     msg->put(NumVdwParams, a2);
04982     msg->put(NumVdwParams, a3);
04983     msg->put(NumVdwParams, a4);
04984 
04985     delete [] a1;
04986     delete [] a2;
04987     delete [] a3;
04988     delete [] a4;
04989   }
04990   
04991   //  Send the vdw pair parameters
04992   msg->put(NumVdwPairParams);
04993   
04994   if (NumVdwPairParams)
04995   {
04996     a1 = new Real[NumVdwPairParams];
04997     a2 = new Real[NumVdwPairParams];
04998     a3 = new Real[NumVdwPairParams];
04999     a4 = new Real[NumVdwPairParams];
05000     i1 = new int[NumVdwPairParams];
05001     i2 = new int[NumVdwPairParams];    
05002 
05003     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
05004          (i1 == NULL) || (i2 == NULL) )
05005     {
05006       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05007     }
05008     
05009     vdw_pair_to_arrays(i1, i2, a1, a2, a3, a4, 0, vdw_pair_tree);
05010     
05011     msg->put(NumVdwPairParams, i1)->put(NumVdwPairParams, i2);
05012     msg->put(NumVdwPairParams, a1);
05013     msg->put(NumVdwPairParams, a2)->put(NumVdwPairParams, a3);
05014     msg->put(NumVdwPairParams, a4);
05015   }
05016   
05017   //  Send the table pair parameters
05018   //printf("Pairs: %i\n", NumTablePairParams);
05019   msg->put(NumTablePairParams);
05020   
05021   if (NumTablePairParams)
05022   {
05023     i1 = new int[NumTablePairParams];
05024     i2 = new int[NumTablePairParams];    
05025     i3 = new int[NumTablePairParams];
05026 
05027     if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05028     {
05029       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05030     }
05031     
05032     table_pair_to_arrays(i1, i2, i3, 0, tab_pair_tree);
05033     
05034     msg->put(NumTablePairParams, i1)->put(NumTablePairParams, i2);
05035     msg->put(NumTablePairParams, i3);
05036   }
05037 
05038   // send the hydrogen bond parameters
05039   // hbondParams.create_message(msg);
05040   msg->end();
05041   delete msg;
05042 }
05043 
05044 /************************************************************************/
05045 /*                  */
05046 /*      FUNCTION receive_Parameters      */
05047 /*                  */
05048 /*  This function is used by all the client processes to receive    */
05049 /*   the structure parameters from the master node.      */
05050 /*                  */
05051 /************************************************************************/
05052 
05053 void Parameters::receive_Parameters(MIStream *msg)
05054 
05055 {
05056   int i, j;      //  Loop counters
05057   Real *a1, *a2, *a3, *a4;  //  Temporary arrays to get data from message in
05058   int *i1, *i2, *i3;      //  Temporary int array to get data from message in
05059   IndexedVdwPair *new_node;  //  New node for vdw pair param tree
05060   IndexedTablePair *new_tab_node;
05061   Real **kvals;      //  Force constant values for dihedrals and impropers
05062   int **nvals;      //  Periodicity values for dihedrals and impropers
05063   Real **deltavals;    //  Phase shift values for dihedrals and impropers
05064 
05065   //  Get the bonded parameters
05066   msg->get(NumBondParams);
05067 
05068   if (NumBondParams)
05069   {
05070     bond_array = new BondValue[NumBondParams];
05071     a1 = new Real[NumBondParams];
05072     a2 = new Real[NumBondParams];
05073 
05074     if ( (bond_array == NULL) || (a1 == NULL) || (a2 == NULL) )
05075     {
05076       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05077     }
05078 
05079     msg->get(NumBondParams, a1);
05080     msg->get(NumBondParams, a2);
05081 
05082     for (i=0; i<NumBondParams; i++)
05083     {
05084       bond_array[i].k = a1[i];
05085       bond_array[i].x0 = a2[i];
05086     }
05087 
05088     delete [] a1;
05089     delete [] a2;
05090   }
05091 
05092   //  Get the angle parameters
05093   msg->get(NumAngleParams);
05094 
05095   if (NumAngleParams)
05096   {
05097     angle_array = new AngleValue[NumAngleParams];
05098     a1 = new Real[NumAngleParams];
05099     a2 = new Real[NumAngleParams];
05100     a3 = new Real[NumAngleParams];
05101     a4 = new Real[NumAngleParams];
05102     i1 = new int[NumAngleParams];
05103 
05104     if ( (angle_array == NULL) || (a1 == NULL) || (a2 == NULL) ||
05105          (a3 == NULL) || (a4 == NULL) || (i1 == NULL))
05106     {
05107       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05108     }
05109 
05110     msg->get(NumAngleParams, a1);
05111     msg->get(NumAngleParams, a2);
05112     msg->get(NumAngleParams, a3);
05113     msg->get(NumAngleParams, a4);
05114     msg->get(NumAngleParams, i1);
05115 
05116     for (i=0; i<NumAngleParams; i++)
05117     {
05118       angle_array[i].k = a1[i];
05119       angle_array[i].theta0 = a2[i];
05120       angle_array[i].k_ub = a3[i];
05121       angle_array[i].r_ub = a4[i];
05122       angle_array[i].normal = i1[i];
05123     }
05124 
05125     delete [] a1;
05126     delete [] a2;
05127     delete [] a3;
05128     delete [] a4;
05129     delete [] i1;
05130   }
05131 
05132   //  Get the dihedral parameters
05133   msg->get(NumDihedralParams);
05134 
05135   if (NumDihedralParams)
05136   {
05137     dihedral_array = new DihedralValue[NumDihedralParams];
05138 
05139     i1 = new int[NumDihedralParams];
05140     kvals = new Real *[MAX_MULTIPLICITY];
05141     nvals = new int *[MAX_MULTIPLICITY];
05142     deltavals = new Real *[MAX_MULTIPLICITY];
05143 
05144     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
05145          (deltavals == NULL) || (dihedral_array == NULL) )
05146     {
05147       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05148     }
05149 
05150     for (i=0; i<MAX_MULTIPLICITY; i++)
05151     {
05152       kvals[i] = new Real[NumDihedralParams];
05153       nvals[i] = new int[NumDihedralParams];
05154       deltavals[i] = new Real[NumDihedralParams];
05155 
05156       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05157       {
05158         NAMD_die("memory allocation failed in Parameters::send_Parameters");
05159       }
05160     }
05161 
05162     msg->get(NumDihedralParams, i1);
05163 
05164     for (i=0; i<MAX_MULTIPLICITY; i++)
05165     {
05166       msg->get(NumDihedralParams, kvals[i]);
05167       msg->get(NumDihedralParams, nvals[i]);
05168       msg->get(NumDihedralParams, deltavals[i]);
05169     }
05170 
05171     for (i=0; i<NumDihedralParams; i++)
05172     {
05173       dihedral_array[i].multiplicity = i1[i];
05174 
05175       for (j=0; j<MAX_MULTIPLICITY; j++)
05176       {
05177         dihedral_array[i].values[j].k = kvals[j][i];
05178         dihedral_array[i].values[j].n = nvals[j][i];
05179         dihedral_array[i].values[j].delta = deltavals[j][i];
05180       }
05181     }
05182 
05183     for (i=0; i<MAX_MULTIPLICITY; i++)
05184     {
05185       delete [] kvals[i];
05186       delete [] nvals[i];
05187       delete [] deltavals[i];
05188     }
05189 
05190     delete [] i1;
05191     delete [] kvals;
05192     delete [] nvals;
05193     delete [] deltavals;
05194   }
05195 
05196   //  Get the improper parameters
05197   msg->get(NumImproperParams);
05198 
05199   if (NumImproperParams)
05200   {
05201     improper_array = new ImproperValue[NumImproperParams];
05202     i1 = new int[NumImproperParams];
05203     kvals = new Real *[MAX_MULTIPLICITY];
05204     nvals = new int *[MAX_MULTIPLICITY];
05205     deltavals = new Real *[MAX_MULTIPLICITY];
05206 
05207     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
05208          (deltavals == NULL) || (improper_array==NULL) )
05209     {
05210       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05211     }
05212 
05213     for (i=0; i<MAX_MULTIPLICITY; i++)
05214     {
05215       kvals[i] = new Real[NumImproperParams];
05216       nvals[i] = new int[NumImproperParams];
05217       deltavals[i] = new Real[NumImproperParams];
05218 
05219       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05220       {
05221         NAMD_die("memory allocation failed in Parameters::send_Parameters");
05222       }
05223     }
05224 
05225     msg->get(NumImproperParams,i1);
05226 
05227     for (i=0; i<MAX_MULTIPLICITY; i++)
05228     {
05229       msg->get(NumImproperParams,kvals[i]);
05230       msg->get(NumImproperParams,nvals[i]);
05231       msg->get(NumImproperParams,deltavals[i]);
05232     }
05233 
05234     for (i=0; i<NumImproperParams; i++)
05235     {
05236       improper_array[i].multiplicity = i1[i];
05237 
05238       for (j=0; j<MAX_MULTIPLICITY; j++)
05239       {
05240         improper_array[i].values[j].k = kvals[j][i];
05241         improper_array[i].values[j].n = nvals[j][i];
05242         improper_array[i].values[j].delta = deltavals[j][i];
05243       }
05244     }
05245 
05246     for (i=0; i<MAX_MULTIPLICITY; i++)
05247     {
05248       delete [] kvals[i];
05249       delete [] nvals[i];
05250       delete [] deltavals[i];
05251     }
05252 
05253     delete [] i1;
05254     delete [] kvals;
05255     delete [] nvals;
05256     delete [] deltavals;
05257   }
05258 
05259   //  Get the crossterm parameters
05260   msg->get(NumCrosstermParams);
05261 
05262   if (NumCrosstermParams)
05263   {
05264     crossterm_array = new CrosstermValue[NumCrosstermParams];
05265 
05266     for (i=0; i<NumCrosstermParams; ++i) {
05267       int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
05268       msg->get(nvals,&crossterm_array[i].c[0][0].d00);
05269     }
05270   }
05271   
05272   //Get the energy table
05273   msg->get(numenerentries);
05274   if (numenerentries > 0) {
05275     //fprintf(stdout, "Getting tables\n");
05276     //fprintf(infofile, "Trying to allocate table\n");
05277     table_ener = new BigReal[numenerentries];
05278     //fprintf(infofile, "Finished table allocation\n");
05279     if (table_ener==NULL)
05280     {
05281       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05282     }
05283 
05284     msg->get(numenerentries, table_ener);
05285   }
05286 
05287   //  Get the vdw parameters
05288   msg->get(NumVdwParams);
05289   msg->get(NumVdwParamsAssigned);
05290 
05291   if (NumVdwParams)
05292   {
05293           atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
05294     vdw_array = new VdwValue[NumVdwParams];
05295     a1 = new Real[NumVdwParams];
05296     a2 = new Real[NumVdwParams];
05297     a3 = new Real[NumVdwParams];
05298     a4 = new Real[NumVdwParams];
05299 
05300     if ( (vdw_array==NULL) || (a1==NULL) || (a2==NULL) || (a3==NULL)
05301              || (a4==NULL) || (atomTypeNames==NULL))
05302     {
05303       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05304     }
05305 
05306     msg->get(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
05307     msg->get(NumVdwParams, a1);
05308     msg->get(NumVdwParams, a2);
05309     msg->get(NumVdwParams, a3);
05310     msg->get(NumVdwParams, a4);
05311 
05312     for (i=0; i<NumVdwParams; i++)
05313     {
05314       vdw_array[i].sigma = a1[i];
05315       vdw_array[i].epsilon = a2[i];
05316       vdw_array[i].sigma14 = a3[i];
05317       vdw_array[i].epsilon14 = a4[i];
05318     }
05319 
05320     delete [] a1;
05321     delete [] a2;
05322     delete [] a3;
05323     delete [] a4;
05324   }
05325   
05326   //  Get the vdw_pair_parameters
05327   msg->get(NumVdwPairParams);
05328   
05329   if (NumVdwPairParams)
05330   {
05331     a1 = new Real[NumVdwPairParams];
05332     a2 = new Real[NumVdwPairParams];
05333     a3 = new Real[NumVdwPairParams];
05334     a4 = new Real[NumVdwPairParams];
05335     i1 = new int[NumVdwPairParams];
05336     i2 = new int[NumVdwPairParams];
05337 
05338     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
05339          (i1 == NULL) || (i2 == NULL) )
05340     {
05341       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05342     }
05343     
05344     msg->get(NumVdwPairParams, i1);
05345     msg->get(NumVdwPairParams, i2);
05346     msg->get(NumVdwPairParams, a1);
05347     msg->get(NumVdwPairParams, a2);
05348     msg->get(NumVdwPairParams, a3);
05349     msg->get(NumVdwPairParams, a4);
05350     
05351     for (i=0; i<NumVdwPairParams; i++)
05352     {
05353       new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
05354       
05355       if (new_node == NULL)
05356       {
05357          NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05358       }
05359       
05360       new_node->ind1 = i1[i];
05361       new_node->ind2 = i2[i];
05362       new_node->A = a1[i];
05363       new_node->A14 = a2[i];
05364       new_node->B = a3[i];
05365       new_node->B14 = a4[i];
05366       new_node->left = NULL;
05367       new_node->right = NULL;
05368       
05369       vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
05370     }
05371     
05372     delete [] i1;
05373     delete [] i2;
05374     delete [] a1;
05375     delete [] a2;
05376     delete [] a3;
05377     delete [] a4;
05378   }
05379  
05380   //  Get the table_pair_parameters
05381   msg->get(NumTablePairParams);
05382   
05383   if (NumTablePairParams)
05384   {
05385     i1 = new int[NumTablePairParams];
05386     i2 = new int[NumTablePairParams];
05387     i3 = new int[NumTablePairParams];
05388 
05389     if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05390     {
05391       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05392     }
05393     
05394     msg->get(NumTablePairParams, i1);
05395     msg->get(NumTablePairParams, i2);
05396     msg->get(NumTablePairParams, i3);
05397     
05398     for (i=0; i<NumTablePairParams; i++)
05399     {
05400       new_tab_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
05401       
05402       if (new_tab_node == NULL)
05403       {
05404          NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05405       }
05406       
05407 //      printf("Adding new table pair with ind1 %i ind2 %i k %i\n", i1[i], i2[i],i3[i]);
05408       new_tab_node->ind1 = i1[i];
05409       new_tab_node->ind2 = i2[i];
05410       new_tab_node->K = i3[i];
05411       new_tab_node->left = NULL;
05412       new_tab_node->right = NULL;
05413       
05414       tab_pair_tree = add_to_indexed_table_pairs(new_tab_node, tab_pair_tree);
05415     }
05416     
05417     delete [] i1;
05418     delete [] i2;
05419     delete [] i3;
05420   }
05421   
05422   // receive the hydrogen bond parameters
05423   // hbondParams.receive_message(msg);
05424 
05425   AllFilesRead = TRUE;
05426 
05427   delete msg;
05428 }
05429 /*      END OF FUNCTION receive_Parameters    */
05430 
05431 /************************************************************************/
05432 /*                  */
05433 /*      FUNCTION convert_vdw_pairs      */
05434 /*                  */
05435 /*  This function converts the linked list of vdw_pairs indexed by  */
05436 /*  atom name into a binary search tree of parameters stored by vdw     */
05437 /*  type index.  This tree is what will be used for real when searching */
05438 /*  for parameters during the simulation.        */
05439 /*                  */
05440 /************************************************************************/
05441 
05442 void Parameters::convert_vdw_pairs()
05443    
05444 {
05445    #ifdef MEM_OPT_VERSION
05446    AtomCstInfo atom_struct;
05447    #else
05448    Atom atom_struct;    //  Dummy structure for getting indexes
05449    #endif
05450    Index index1, index2;  //  Indexes for the two atoms
05451    IndexedVdwPair *new_node;  //  New node for tree
05452    struct vdw_pair_params *ptr, *next;  //  Pointers for traversing list
05453    
05454    ptr = vdw_pairp;
05455    
05456    //  Go down then entire list and insert each node into the 
05457    //  binary search tree
05458    while (ptr != NULL)
05459    {
05460       new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
05461       
05462       if (new_node == NULL)
05463       {
05464    NAMD_die("memory allocation failed in Parameters::convert_vdw_pairs");
05465       }
05466       
05467       //  Get the vdw indexes for the two atoms.  This is kind of a hack
05468       //  using the goofy Atom structure, but hey, it works
05469       assign_vdw_index(ptr->atom1name, &atom_struct);
05470       index1 = atom_struct.vdw_type;
05471       assign_vdw_index(ptr->atom2name, &atom_struct);
05472       index2 = atom_struct.vdw_type;
05473       
05474       if (index1 > index2)
05475       {
05476    new_node->ind1 = index2;
05477    new_node->ind2 = index1;
05478       }
05479       else
05480       {
05481    new_node->ind1 = index1;
05482    new_node->ind2 = index2;
05483       }
05484            
05485       new_node->A = ptr->A;
05486       new_node->B = ptr->B;
05487       new_node->A14 = ptr->A14;
05488       new_node->B14 = ptr->B14;
05489       
05490       new_node->left = NULL;
05491       new_node->right = NULL;
05492       
05493       //  Add it to the tree
05494       vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
05495       
05496       //  Free the current node and move to the next
05497       next = ptr->next;
05498       
05499       delete ptr;
05500       
05501       ptr = next;
05502    }
05503    
05504    vdw_pairp = NULL;
05505 
05506 }
05507 /*      END OF FUNCTION convert_vdw_pairs    */
05508 
05509 /************************************************************************/
05510 /*                  */
05511 /*      FUNCTION convert_table_pairs      */
05512 /*                  */
05513 /*  This function converts the linked list of table_pairs indexed by  */
05514 /*  atom name into a binary search tree of parameters stored by table     */
05515 /*  type index.  This tree is what will be used for real when searching */
05516 /*  for parameters during the simulation.        */
05517 /*                  */
05518 /************************************************************************/
05519 
05520 void Parameters::convert_table_pairs()
05521    
05522 {
05523    #ifdef MEM_OPT_VERSION
05524    AtomCstInfo atom_struct;
05525    #else
05526    Atom atom_struct;    //  Dummy structure for getting indexes
05527    #endif
05528    Index index1, index2;  //  Indexes for the two atoms
05529    IndexedTablePair *new_node;  //  New node for tree
05530    struct table_pair_params *ptr, *next;  //  Pointers for traversing list
05531    
05532    ptr = table_pairp;
05533    
05534    //  Go down then entire list and insert each node into the 
05535    //  binary search tree
05536    while (ptr != NULL)
05537    {
05538       new_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
05539       
05540       if (new_node == NULL)
05541       {
05542    NAMD_die("memory allocation failed in Parameters::convert_table_pairs");
05543       }
05544       
05545       //  Get the vdw indexes for the two atoms.  This is kind of a hack
05546       //  using the goofy Atom structure, but hey, it works
05547       assign_vdw_index(ptr->atom1name, &atom_struct);
05548       index1 = atom_struct.vdw_type;
05549       assign_vdw_index(ptr->atom2name, &atom_struct);
05550       index2 = atom_struct.vdw_type;
05551 
05552 //      printf("Adding new table pair with index1 %i, index2 %i, k %i\n", index1, index2, ptr->K);
05553       
05554       if (index1 > index2)
05555       {
05556    new_node->ind1 = index2;
05557    new_node->ind2 = index1;
05558       }
05559       else
05560       {
05561    new_node->ind1 = index1;
05562    new_node->ind2 = index2;
05563       }
05564            
05565       new_node->K = ptr->K;
05566 
05567       new_node->left = NULL;
05568       new_node->right = NULL;
05569       
05570       //  Add it to the tree
05571       tab_pair_tree = add_to_indexed_table_pairs(new_node, tab_pair_tree);
05572       
05573       //  Free the current node and move to the next
05574       next = ptr->next;
05575       
05576       delete ptr;
05577       
05578       ptr = next;
05579    }
05580    
05581    table_pairp = NULL;
05582 
05583 }
05584 /*      END OF FUNCTION convert_table_pairs    */
05585 
05586 /************************************************************************/
05587 /*                  */
05588 /*      FUNCTION add_to_indexed_table_pairs    */
05589 /*                  */
05590 /*   INPUTS:                */
05591 /*  new_node - new node to be added to the tree      */
05592 /*  tree - tree to add the node to          */
05593 /*                  */
05594 /*  This is a recursive function that adds a node to the    */
05595 /*   binary search tree of table_pair parameters        */
05596 /*                  */
05597 /************************************************************************/
05598 
05599 IndexedTablePair *Parameters::add_to_indexed_table_pairs(IndexedTablePair *new_node,
05600                  IndexedTablePair *tree)
05601    
05602 {
05603    if (tree == NULL)
05604       return(new_node);
05605    
05606    if ( (new_node->ind1 < tree->ind1) || 
05607         ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
05608    {
05609       tree->left = add_to_indexed_table_pairs(new_node, tree->left);
05610    }
05611    else
05612    {
05613       tree->right = add_to_indexed_table_pairs(new_node, tree->right);
05614    }
05615    
05616    return(tree);
05617 }
05618 /*      END OF FUNCTION add_to_indexed_table_pairs  */
05619 
05620 /************************************************************************/
05621 /*                  */
05622 /*      FUNCTION add_to_indexed_vdw_pairs    */
05623 /*                  */
05624 /*   INPUTS:                */
05625 /*  new_node - new node to be added to the tree      */
05626 /*  tree - tree to add the node to          */
05627 /*                  */
05628 /*  This is a recursive function that adds a node to the    */
05629 /*   binary search tree of vdw_pair parameters        */
05630 /*                  */
05631 /************************************************************************/
05632 
05633 IndexedVdwPair *Parameters::add_to_indexed_vdw_pairs(IndexedVdwPair *new_node,
05634                  IndexedVdwPair *tree)
05635    
05636 {
05637    if (tree == NULL)
05638       return(new_node);
05639    
05640    if ( (new_node->ind1 < tree->ind1) || 
05641         ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
05642    {
05643       tree->left = add_to_indexed_vdw_pairs(new_node, tree->left);
05644    }
05645    else
05646    {
05647       tree->right = add_to_indexed_vdw_pairs(new_node, tree->right);
05648    }
05649    
05650    return(tree);
05651 }
05652 /*      END OF FUNCTION add_to_indexed_vdw_pairs  */
05653 
05654 /************************************************************************/
05655 /*                  */
05656 /*      FUNCTION vdw_pair_to_arrays      */
05657 /*                  */
05658 /*   INPUTS:                */
05659 /*  ind1_array - Array of index 1 values        */
05660 /*  ind2_array - Array of index 2 values        */
05661 /*  A - Array of A values            */
05662 /*  A14 - Array of A 1-4 values          */
05663 /*  B - Array of B values            */
05664 /*  B14 - Array of B 1-4 values          */
05665 /*  arr_index - current position in arrays        */
05666 /*  tree - tree to traverse            */
05667 /*                  */
05668 /*  This is a recursive function that places all the entries of     */
05669 /*   the tree passed in into arrays of values.  This is done so that    */
05670 /*   the parameters can be sent from the master node to the other       */
05671 /*   nodes.                */
05672 /*                  */
05673 /************************************************************************/
05674 
05675 int Parameters::vdw_pair_to_arrays(int *ind1_array, int *ind2_array,
05676           Real *A, Real *A14,
05677           Real *B, Real *B14,
05678           int arr_index, IndexedVdwPair *tree)
05679       
05680 {
05681    if (tree == NULL)
05682       return(arr_index);
05683    
05684    ind1_array[arr_index] = tree->ind1;
05685    ind2_array[arr_index] = tree->ind2;
05686    A[arr_index] = tree->A;
05687    A14[arr_index] = tree->A14;
05688    B[arr_index] = tree->B;
05689    B14[arr_index] = tree->B14;
05690    
05691    arr_index++;
05692    
05693    arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
05694           arr_index, tree->left);
05695    arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
05696           arr_index, tree->right);
05697    
05698    return(arr_index);
05699 }
05700 /*      END OF FUNCTION vdw_pair_to_arrays    */
05701 
05702 /************************************************************************/
05703 /*                  */
05704 /*      FUNCTION table_pair_to_arrays      */
05705 /*                  */
05706 /*   INPUTS:                */
05707 /*  ind1_array - Array of index 1 values        */
05708 /*  ind2_array - Array of index 2 values        */
05709 /*  K - Array of K values            */
05710 /*                  */
05711 /*  This is a recursive function that places all the entries of     */
05712 /*   the tree passed in into arrays of values.  This is done so that    */
05713 /*   the parameters can be sent from the master node to the other       */
05714 /*   nodes.                */
05715 /*                  */
05716 /************************************************************************/
05717 
05718 int Parameters::table_pair_to_arrays(int *ind1_array, int *ind2_array,
05719           int *K,
05720           int arr_index, IndexedTablePair *tree)
05721       
05722 {
05723    if (tree == NULL)
05724       return(arr_index);
05725    
05726    ind1_array[arr_index] = tree->ind1;
05727    ind2_array[arr_index] = tree->ind2;
05728    K[arr_index] = tree->K;
05729 
05730    arr_index++;
05731    
05732    arr_index = table_pair_to_arrays(ind1_array, ind2_array, K,
05733           arr_index, tree->left);
05734    arr_index = table_pair_to_arrays(ind1_array, ind2_array, K,
05735           arr_index, tree->right);
05736    
05737    return(arr_index);
05738 }
05739 /*      END OF FUNCTION table_pair_to_arrays    */
05740 
05741 /************************************************************************/
05742 /*                  */
05743 /*      FUNCTION Parameters        */
05744 /*                  */
05745 /*  This is the constructor for reading AMBER parameter data    */
05746 /*                  */
05747 /************************************************************************/
05748 
05749 Parameters::Parameters(Ambertoppar *amber_data, BigReal vdw14)
05750 {
05751   initialize();
05752 
05753   // Read in parm parameters
05754   read_parm(amber_data,vdw14);
05755 }
05756 /*      END OF FUNCTION Parameters      */
05757 
05758 
05759 /************************************************************************/
05760 /*                  */
05761 /*      FUNCTION read_parm   */
05762 /*                  */
05763 /*   INPUTS:                */
05764 /*  amber_data - AMBER data structure    */
05765 /*                  */
05766 /*  This function copys AMBER parameter data to the corresponding data  */
05767 /*   structures      */
05768 /*                  */
05769 /************************************************************************/
05770 
05771 void Parameters::read_parm(Ambertoppar *amber_data, BigReal vdw14)
05772 { 
05773   int i,j,ico,numtype,mul;
05774   IndexedVdwPair *new_node;
05775 
05776   if (!amber_data->data_read)
05777     NAMD_die("No data read from parm file yet!");
05778 
05779   // Copy bond parameters
05780   NumBondParams = amber_data->Numbnd;
05781   if (NumBondParams)
05782   { bond_array = new BondValue[NumBondParams];
05783     if (bond_array == NULL)
05784       NAMD_die("memory allocation of bond_array failed!");
05785   }
05786   for (i=0; i<NumBondParams; ++i)
05787   { bond_array[i].k = amber_data->Rk[i];
05788     bond_array[i].x0 = amber_data->Req[i];
05789   }
05790 
05791   // Copy angle parameters
05792   NumAngleParams = amber_data->Numang;
05793   if (NumAngleParams)
05794   { angle_array = new AngleValue[NumAngleParams];
05795     if (angle_array == NULL)
05796       NAMD_die("memory allocation of angle_array failed!");
05797   }
05798   for (i=0; i<NumAngleParams; ++i)
05799   { angle_array[i].k = amber_data->Tk[i];
05800     angle_array[i].theta0 = amber_data->Teq[i];
05801     // Amber has no k_ub and r_ub for angle parameters, so they're set to 0
05802     angle_array[i].k_ub = angle_array[i].r_ub = 0;
05803     // All angles are harmonic
05804     angle_array[i].normal = 1;
05805   }
05806 
05807   // Copy dihedral parameters
05808   // Note: If the periodicity is negative, it means the following
05809   //  entry is another term in a multitermed dihedral, and the
05810   //  absolute value is the true periodicity; in this case the
05811   //  following entry in "dihedral_array" should be left empty,
05812   //  NOT be skipped, in order to keep the next dihedral's index
05813   //  unchanged.
05814   NumDihedralParams = amber_data->Nptra;
05815   if (NumDihedralParams)
05816   { dihedral_array = new DihedralValue[amber_data->Nptra];
05817     if (dihedral_array == NULL)
05818       NAMD_die("memory allocation of dihedral_array failed!");
05819   }
05820   mul = 0;
05821   for (i=0; i<NumDihedralParams; ++i)
05822   { dihedral_array[i-mul].values[mul].k = amber_data->Pk[i];
05823     dihedral_array[i-mul].values[mul].n = int(fabs(amber_data->Pn[i])+0.5);
05824     if (dihedral_array[i-mul].values[mul].n == 0)
05825     { char err_msg[128];
05826       sprintf(err_msg, "The periodicity of dihedral # %d is zero!", i+1);
05827       NAMD_die(err_msg);
05828     }
05829     dihedral_array[i-mul].values[mul].delta = amber_data->Phase[i];
05830     // If the periodicity is positive, it means the following
05831     // entry is a new dihedral term.
05832     if (amber_data->Pn[i] > 0)
05833     { dihedral_array[i-mul].multiplicity = mul+1;
05834       mul = 0;
05835     }
05836     else if (++mul >= MAX_MULTIPLICITY)
05837     { char err_msg[181];
05838       sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
05839          mul+1, MAX_MULTIPLICITY);
05840       NAMD_die(err_msg);
05841     }
05842   }
05843   if (mul > 0)
05844     NAMD_die("Negative periodicity in the last dihedral entry!");
05845 
05846   // Copy VDW parameters: AMBER explicitly gives VDW parameters between every
05847   // 2 atom types, so the data are copied to vdw_pair_tree
05848   // In AMBER, all 1-4 VDW interactions are divided by factor vdw14
05849   NumVdwParamsAssigned = numtype = amber_data->Ntypes; // Number of atom types
05850   NumVdwPairParams = numtype * (numtype+1) / 2;
05851   for (i=0; i<numtype; ++i)
05852     for (j=i; j<numtype; ++j)
05853     { new_node = new IndexedVdwPair;
05854       if (new_node == NULL)
05855         NAMD_die("memory allocation of vdw_pair_tree failed!");
05856       new_node->ind1 = i;
05857       new_node->ind2 = j;
05858       new_node->left = new_node->right = NULL;
05859       // ico is the index of interaction between atom type i and j into
05860       // the parameter arrays. If it's positive, the interaction is
05861       // 6-12 VDW; otherwise it's 10-12 H-bond interaction. NAMD doesn't
05862       // have 10-12 term, so if ico is negative, then the 10-12
05863       // coefficients must be 0, otherwise die.
05864       ico = amber_data->Cno[numtype*i+j];
05865       if (ico>0)
05866       { new_node->A = amber_data->Cn1[ico-1];
05867         new_node->A14 = new_node->A / vdw14;
05868         new_node->B = amber_data->Cn2[ico-1];
05869         new_node->B14 = new_node->B / vdw14;
05870       }
05871       else if (amber_data->HB12[abs(ico)-1]==0.0 && amber_data->HB6[abs(ico)-1]==0.0)
05872       { new_node->A = new_node->A14 = new_node->B = new_node->B14 = 0.0;
05873         iout << iWARN << "Encounter 10-12 H-bond term\n";
05874       }
05875       else
05876         NAMD_die("Encounter non-zero 10-12 H-bond term!");
05877       // Add the node to the binary tree
05878       vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
05879     }
05880 }
05881 /*      END OF FUNCTION read_parm    */
05882 
05883 /************************************************************************/
05884 /*                                                                      */
05885 /*      FUNCTION Parameters                                             */
05886 /*                                                                      */
05887 /*  This is the constructor for reading GROMACS parameter data          */
05888 /*                                                                      */
05889 /************************************************************************/
05890 
05891 Parameters::Parameters(const GromacsTopFile *gf, Bool min)
05892 {
05893   initialize();
05894 
05895   // Read in parm parameters
05896   read_parm(gf,min);
05897 }
05898 /*      END OF FUNCTION Parameters      */
05899 
05900 
05901 /************************************************************************/
05902 /*                                                                      */
05903 /*      FUNCTION read_parm                                              */
05904 /*                                                                      */
05905 /*   INPUTS:                                                            */
05906 /*  gf - GROMACS topology file                                          */
05907 /*                                                                      */
05908 /* This function copys GROMACS parameter data to the corresponding data */
05909 /*   structures                                                         */
05910 /*                                                                      */
05911 /************************************************************************/
05912 
05913 void Parameters::read_parm(const GromacsTopFile *gf, Bool min)
05914 { 
05915   int numtype;
05916   IndexedVdwPair *new_node;
05917   int i,j,funct;
05918   Real test1,test2;
05919 
05920   // Allocate space for all of the arrays first
05921   NumBondParams = gf->getNumBondParams();
05922   NumAngleParams = gf->getNumAngleParams();
05923   NumDihedralParams = gf->getNumDihedralParams();
05924   if (NumBondParams) {
05925     bond_array = new BondValue[NumBondParams];
05926     if (bond_array == NULL)
05927       NAMD_die("memory allocation of bond_array failed!");
05928   }
05929   if (NumDihedralParams) {
05930     dihedral_array = new DihedralValue[NumDihedralParams];
05931     if (dihedral_array == NULL)
05932       NAMD_die("memory allocation of dihedral_array failed!");
05933   }
05934   if (NumAngleParams) {
05935     angle_array = new AngleValue[NumAngleParams];
05936     if (angle_array == NULL)
05937       NAMD_die("memory allocation of angle_array failed!");
05938   }
05939 
05940   // Copy bond parameters
05941   // XXX Warning: we are discarding the GROMACS function type - since
05942   // NAMD does not let us choose between different spring models.
05943   for (i=0;i<NumBondParams;i++) {
05944     Real x0,k;
05945     gf->getBondParams(i,
05946                       &x0, // the bond length
05947                       &k,  // the spring constant
05948                       &funct);           // discarded
05949     bond_array[i].x0 = x0;
05950     bond_array[i].k = k;
05951   }
05952 
05953   // Copy angle parameters
05954   // XXX Warning: we are discarding the GROMACS function type here
05955   // too.
05956   for (i=0;i<NumAngleParams;i++) {
05957     Real theta0,k;
05958     gf->getAngleParams(i,
05959                        &theta0, // the angle size
05960                        &k,      // the spring constant
05961                        &funct);                // discarded
05962     angle_array[i].theta0 = theta0*PI/180;
05963     angle_array[i].k = k;
05964     // Gromacs has no Urey-Bradley angle parameters, so they're set to 0
05965     angle_array[i].k_ub = angle_array[i].r_ub = 0;
05966     angle_array[i].normal = 1;
05967   }
05968 
05969   // Copy dihedral parameters
05970   // Here we use the function type (carefully)
05971   for (i=0; i<NumDihedralParams; ++i) { // iterate over all dihedral angles
05972     Real c[6];
05973     int num_periods; // number of periods in one full rotation
05974     int funct;
05975 
05976     gf->getDihedralParams(i,c,&num_periods,&funct); // get the parameters
05977     
05978     switch(funct) {
05979     case 1: 
05980       dihedral_array[i].values[0].delta = c[0]*PI/180; // the phase shift
05981       dihedral_array[i].values[0].k = c[1]; // the spring constant
05982       dihedral_array[i].values[0].n = num_periods; // the periodicity
05983       dihedral_array[i].multiplicity = 1;
05984       break;
05985     case 2: 
05986       dihedral_array[i].values[0].delta = c[0]*PI/180; // the phase shift
05987       dihedral_array[i].values[0].k = c[1]; // the spring constant
05988       dihedral_array[i].values[0].n = 0; // 'periodicity'=0 for impropers
05989       dihedral_array[i].multiplicity = 1;
05990       break;
05991     case 3: 
05992 
05993       // first a quick check to make sure this is legal
05994       if(MAX_MULTIPLICITY < 5)
05995         NAMD_die("I can't do RB dihedrals with MAX_MULTIPLICITY < 5");
05996       dihedral_array[i].multiplicity = 5;
05997 
05998       // Next we negate every other term, since GROMACS does this
05999       // silly thing with psi = 180 - phi:
06000       for(j=0;j<6;j++) {
06001         if(j%2 == 1) c[j] = -c[j];
06002       }
06003 
06004       // Now fill up all the terms.  Each is k(1 + cos(n*x - delta))
06005       // so we first let delta = 0 and let n range from 1-6:
06006       for(j=0;j<5;j++) dihedral_array[i].values[j].delta = 0;
06007       for(j=0;j<5;j++) dihedral_array[i].values[j].n = j+1;
06008 
06009       // and now we have a sum of kn(1+cos(nx))
06010       // Gromacs RB gives you a sum of cn(cos(x))^n, so we have to
06011       // use trigonometry to compute the kn:
06012       dihedral_array[i].values[0].k =    1*c[1] + 3/4.*c[3] + 10/16.*c[5];
06013       dihedral_array[i].values[1].k =       1/2.*c[2] + 4/8.*c[4]        ;
06014       dihedral_array[i].values[2].k =             1/4.*c[3] +  5/16.*c[5];
06015       dihedral_array[i].values[3].k =                   1/8.*c[4]        ;
06016       dihedral_array[i].values[4].k =                          1/16.*c[5];
06017 
06018       // That was a lot of math, so we had better check it:
06019       // The constant term (which is missing) is c0 + 1/2 c2 + 3/8 c4
06020       test1 = 0;
06021       for(j=5;j>=0;j--) { // sum (..(c5 cos x + c4) cos x + c3)..) + c1
06022         test1 *= cos(0.5);
06023         test1 += c[j];
06024       }
06025 
06026       test2 = c[0]+1/2.*c[2]+3/8.*c[4];
06027       for(j=0;j<5;j++) { // sum k1 cos x + k2 cos 2x + ... 
06028         test2 += dihedral_array[i].values[j].k * cos((j+1)*0.5);
06029       }
06030 
06031       if(fabs(test1-test2) > 0.0001)
06032         NAMD_die("Internal error: failed to handle RB dihedrals");
06033       
06034       // Turn this on to have a look at the values if you *still*
06035       // don't believe that they are right!
06036 
06037       /*      iout << iINFO << "k: ";
06038               for(j=0;j<5;j++)
06039               iout  << dihedral_array[i].values[j].k << " ";
06040               iout << "\n" << endi;
06041               
06042               iout << iINFO << "c: ";
06043               for(j=0;j<6;j++)
06044               iout  << c[j] << " ";
06045               iout << "\n" << endi;*/
06046       
06047       break;
06048     default:
06049       NAMD_die("unknown dihedral type found");
06050     }
06051   }
06052 
06053   // Copy VDW parameters.
06054 
06055   Bool warned=false; // warned the user about extra LJ term yet?
06056 
06057   NumVdwParamsAssigned = numtype = gf->getNumAtomParams(); // # of atom types
06058   NumVdwPairParams = numtype * (numtype+1) / 2;
06059   for (i=0; i<numtype; i++) {
06060     for (j=i; j<numtype; j++) {
06061 
06062       // set up a node to put one set of VDW parameters in
06063       new_node = new IndexedVdwPair;
06064       if (new_node == NULL)
06065         NAMD_die("memory allocation of vdw_pair_tree failed!");
06066       new_node->ind1 = i;
06067       new_node->ind2 = j;
06068       new_node->left = new_node->right = NULL;
06069 
06070       gf->getVDWParams(i,j, &(new_node->B), &(new_node->A),
06071                        &(new_node->B14), &(new_node->A14));
06072 
06073       /* If we have any VDW radii equal to zero, atoms can just sit on
06074          each other during minimization.  So, we'll give a minimum of
06075          1.0 kcal*A^12 to the LJ-repulsion when we are minimizing.
06076          But a warning should be displayed to the user... */
06077       if(min && ( fabs(new_node->A) < 1.0 )) {
06078         new_node->A = 1.0;
06079         if(!warned) {
06080           iout << iWARN <<
06081             "Adding small LJ repulsion term to some atoms.\n" << endi;
06082           warned=true;
06083         }
06084       }
06085 
06086       vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
06087     }
06088   }
06089 }
06090 /*      END OF FUNCTION read_parm    */
06091 
06092 /************************************************************************/
06093 /*                                                                      */
06094 /*      FUNCTION read_ener_table                                          */
06095 /*                                                                      */
06096 /*   INPUTS:                                                            */
06097 /*  simParams -- Simulation Parameters   */
06098 /*                                                                      */
06099 /*   This function reads energy tables from a file and places them into                                                       */
06100 /*   memory.                                                                      */
06101 /************************************************************************/
06102 
06103 void Parameters::read_ener_table(SimParameters *simParams) {
06104         char* table_file = simParams->tabulatedEnergiesFile;
06105   char* interp_type = simParams->tableInterpType;
06106         FILE* enertable;
06107         enertable = fopen(table_file, "r");
06108 
06109         if (enertable == NULL) {
06110                 NAMD_die("ERROR: Could not open energy table file!\n");
06111         }
06112 
06113         char tableline[121];
06114   char* newtypename;
06115   int numtypes;
06116         int atom1;
06117         int atom2;
06118         int distbin;
06119   int readcount;
06120         Real dist;
06121         Real energy;
06122         Real force;
06123         Real table_spacing;
06124         Real maxdist;
06125 
06126 /* First find the header line and set the variables we need */
06127         iout << "Beginning table read\n" << endi;
06128         while(fgets(tableline,120,enertable)) {
06129                 if (strncmp(tableline,"#",1)==0) {
06130                         continue;
06131                 }
06132     readcount = sscanf(tableline, "%i %f %f", &numtypes, &table_spacing, &maxdist);
06133     if (readcount != 3) {
06134       NAMD_die("ERROR: Couldn't parse table header line\n");
06135     }
06136     break;
06137   }
06138 
06139   if (maxdist < simParams->cutoff) {
06140     NAMD_die("Tabulated energies must at least extend to the cutoff distance\n");
06141   }
06142 
06143         if (maxdist > simParams->cutoff) {
06144                 iout << "Info: Reducing max table size to match nonbond cutoff\n";
06145                 maxdist = ceil(simParams->cutoff);
06146         }
06147 
06148 /* Now allocate memory for the table; we know what we should be getting */
06149         numenerentries = 2 * numtypes * int (mynearbyint(maxdist/table_spacing) + 1);
06150         //Set up a full energy lookup table from a file
06151         //Allocate the table; layout is atom1 x atom2 x distance energy force
06152         fprintf(stdout, "Table has %i entries\n",numenerentries);
06153         //iout << "Allocating original energy table\n" << endi;
06154         if ( table_ener ) delete [] table_ener;
06155         table_ener = new BigReal[numenerentries];
06156   if ( table_types ) delete [] table_types;
06157   table_types = new char*[numtypes];
06158 
06159 /* Finally, start reading the table */
06160   int numtypesread = 0; //Number of types read so far
06161 
06162         while(fgets(tableline,120,enertable)) {
06163                 if (strncmp(tableline,"#",1)==0) {
06164                         continue;
06165                 }
06166     if (strncmp(tableline,"TYPE",4)==0) {
06167       // Read a new type
06168       newtypename = new char[5];
06169       int readcount = sscanf(tableline, "%*s %s", newtypename);
06170       if (readcount != 1) {
06171         NAMD_die("Couldn't read interaction type from TYPE line\n");
06172       }
06173 //      printf("Setting table_types[%i] to %s\n", numtypesread, newtypename);
06174       table_types[numtypesread] = newtypename;
06175 
06176       if (numtypesread == numtypes) {
06177         NAMD_die("Error: Number of types in table doesn't match table header\n");
06178       }
06179 
06180       // Read the current energy type with the proper interpolation
06181       if (!strncasecmp(interp_type, "linear", 6)) {
06182         if (read_energy_type(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
06183           char err_msg[512];
06184           sprintf(err_msg, "Failed to read table energy (linear) type %s\n", newtypename);
06185           NAMD_die(err_msg);
06186         }
06187       } else if (!strncasecmp(interp_type, "cubic", 5)) {
06188         if (read_energy_type_bothcubspline(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
06189           char err_msg[512];
06190           sprintf(err_msg, "Failed to read table energy (cubic) type %s\n", newtypename);
06191           NAMD_die(err_msg);
06192         }
06193       } else {
06194         NAMD_die("Table interpolation type must be linear or cubic\n");
06195       }
06196 
06197       numtypesread++;
06198       continue;
06199     }
06200     //char err_msg[512];
06201     //sprintf(err_msg, "Unrecognized line in energy table file: %s\n", tableline);
06202     //NAMD_die(err_msg);
06203   }
06204 
06205   // See if we got what we expected
06206   if (numtypesread != numtypes) {
06207     char err_msg[512];
06208     sprintf(err_msg, "ERROR: Expected %i tabulated energy types but got %i\n", numtypes, numtypesread);
06209     NAMD_die(err_msg);
06210   }
06211 
06212   // Move the necessary information into simParams
06213   simParams->tableNumTypes = numtypes;
06214   simParams->tableSpacing = table_spacing;
06215   simParams->tableMaxDist = maxdist;
06216   tablenumtypes = numtypes;
06217 
06218   /*
06219 xxxxxx
06220         int numtypes = simParams->tableNumTypes;
06221 
06222         //This parameter controls the table spacing
06223         BigReal table_spacing = 0.01;
06224         BigReal maxdist = 20.0;
06225         if (maxdist > simParams->cutoff) {
06226                 iout << "Info: Reducing max table size to match nonbond cutoff\n";
06227                 maxdist = ceil(simParams->cutoff);
06228         }
06229 
06230         numenerentries = (numtypes + 1) * numtypes * int (ceil(maxdist/table_spacing));
06231 //      fprintf(stdout, "Table arithmetic: maxdist %f, table_spacing %f, numtypes %i, numentries %i\n", maxdist, table_spacing, numtypes, numenerentries);
06232         columnsize = 2 * mynearbyint(maxdist/table_spacing);
06233         rowsize = numtypes * columnsize;
06234         //Set up a full energy lookup table from a file
06235         //Allocate the table; layout is atom1 x atom2 x distance energy force
06236         fprintf(stdout, "Table has %i entries\n",numenerentries);
06237         //iout << "Allocating original energy table\n" << endi;
06238         if ( table_ener ) delete [] table_ener;
06239         table_ener = new Real[numenerentries];
06240         //
06241         //Set sentinel values for uninitialized table energies
06242         for (int i=0 ; i<numenerentries ; i++) {
06243                 table_ener[i] = 1337.0;
06244         }
06245         Real compval = 1337.0;
06246 
06247         //    iout << "Entering table reading\n" << endi;
06248         //iout << "Finished allocating table\n" << endi;
06249 
06250         //Counter to make sure we read the right number of energy entries
06251         int numentries = 0;
06252 
06253         //Now, start reading from the file
06254         char* table_file = simParams->tabulatedEnergiesFile;
06255         FILE* enertable;
06256         enertable = fopen(table_file, "r");
06257 
06258         if (enertable == NULL) {
06259                 NAMD_die("ERROR: Could not open energy table file!\n");
06260         }
06261 
06262         char tableline[121];
06263         int atom1;
06264         int atom2;
06265         int distbin;
06266         Real dist;
06267         Real energy;
06268         Real force;
06269 
06270         iout << "Beginning table read\n" << endi;
06271         //Read the table entries
06272         while(fgets(tableline,120,enertable)) {
06273 //              IOut << "Processing line " << tableline << "\n" << endi;
06274                 if (strncmp(tableline,"#",1)==0) {
06275                         continue;
06276                 }
06277 
06278 
06279                 sscanf(tableline, "%i %i %f %f %f\n", &atom1, &atom2, &dist, &energy, &force);
06280                 distbin = int(mynearbyint(dist/table_spacing));
06281 //              fprintf(stdout, "Working on atoms %i and %i at distance %f\n",atom1,atom2,dist);
06282                 if ((atom2 > atom1) || (distbin > int(mynearbyint(maxdist/table_spacing)))) {
06283 //                      fprintf(stdout,"Rejected\n");
06284 //                      fprintf(stdout, "Error: Throwing out energy line beyond bounds\n");
06285         //              fprintf(stdout, "Bad line: Atom 1: %i Atom 2: %i Distance Bin: %i Max Distance Bin: %i \n", atom1, atom2, distbin, int(mynearbyint(maxdist/table_spacing)));
06286                 } else {
06287                         //The magic formula for the number of columns from previous rows
06288                         //in the triangular matrix is (2ni+i-i**2)/2
06289                         //Here n is the number of types, and i is atom2
06290 //                      fprintf(stdout, "Input: atom1 %f atom2 %f columnsize %f \n", float(atom1), float(atom2), float(columnsize));
06291 //                      fprintf(stdout, "Testing arithmetic: Part1: %i Part2: %i Part3: %i Total: %i\n", columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2), columnsize*(atom1-atom2), 2*distbin, columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2) + columnsize*(atom1-atom2) + 2*distbin - 2);
06292                         int eneraddress = columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2) + columnsize*(atom1-atom2) + 2*distbin - 2;
06293                         int forceaddress = eneraddress + 1;
06294 //                              fprintf(stdout, "Tableloc: %i Atom 1: %i Atom 2: %i Distance Bin: %i Energy: %f Force: %f\n", eneraddress, atom1, atom2, distbin, table_ener[eneraddress], table_ener[forceaddress]);
06295                 fflush(stdout);
06296 //                      fprintf(stdout, "Checking for dupes: Looking at: %f %f \n", table_ener[eneraddress], table_ener[forceaddress]);
06297                         if ((table_ener[eneraddress] == compval && table_ener[forceaddress] == compval)) {
06298                                 numentries++;
06299                                 table_ener[eneraddress] = energy;
06300                                 table_ener[forceaddress] = force;
06301 //                              fprintf(stdout, "Tableloc: %i Atom 1: %i Atom 2: %i Distance Bin: %i Energy: %f Force: %f\n", eneraddress, atom1, atom2, distbin, table_ener[eneraddress], table_ener[forceaddress]);
06302                                 //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin] = energy;
06303                                 //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin + 1] = force;
06304 //                              fflush(stdout);
06305                         } else {
06306 //                              fprintf(stdout,"Rejecting duplicate entry\n");
06307                         }
06308                 }
06309                 //      iout << "Done with line\n"<< endi;
06310         }
06311 
06312         //    int should = numtypes * numtypes * (maxdist/table_spacing);
06313         //    cout << "Entries: " << numentries << " should be " << should << "\n" << endi;
06314 //      int exptypes = ceil((numtypes+1) * numtypes * (maxdist/table_spacing));
06315 //fprintf(stdout,"numtypes: %i maxdist: %f table_spacing: %f exptypes: %i\n",numtypes,maxdist,table_spacing);
06316         if (numentries != int(numenerentries/2)) {
06317                 fprintf(stdout,"ERROR: Expected %i entries but got %i\n",numenerentries/2, numentries);
06318                 NAMD_die("Number of energy table entries did not match expected value\n");
06319         }
06320         //      iout << "Done with table\n"<< endi;
06321         fprintf(stdout, "Numenerentries: %i\n",numenerentries/2);
06322   */
06323 } /* END of function read_ener_table */ 
06324 
06325 /**************************************************************************
06326  * FUNCTION read_energy_type_bothcubspline
06327  *
06328  * Read a single type block from an energy table file, using cubic spline interpolation
06329  * Unlike _cubspline, the forces are interpolated separately
06330  *
06331  * Inputs:
06332  *  enertable - File stream positioned at the start of the type block
06333  *  typeindex  - integer index of current type
06334  *  table_ener - pointer to array to be filled with table entries
06335  *  table_spacing - Spacing between table points (A)
06336  *  maxdist - Longest distance needed in table
06337  *
06338  * Return values:
06339  *  0 on normal exit
06340  *  1 if not enough entries were present to fill out the table
06341  *
06342  *  Note: enertable should be left positioned immediately BEFORE the next
06343  *  TYPE block starts
06344  *  **********************************************************************/
06345 
06346 int Parameters::read_energy_type_bothcubspline(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
06347 
06348   char tableline[120];
06349   int i,j;
06350 
06351   /* Last values read from table */
06352   BigReal readdist;
06353   BigReal readener;
06354   BigReal readforce;
06355   BigReal spacing;
06356 //  BigReal readforce;
06357   BigReal lastdist;
06358 //  BigReal lastener;
06359 //  BigReal lastforce;
06360 //  readdist = -1.0;
06361 //  readener = 0.0;
06362 //  readforce = 0.0;
06363 
06364   /* Create arrays for holding the input values */
06365   std::vector<BigReal>  dists;
06366   std::vector<BigReal> enervalues;
06367   std::vector<BigReal> forcevalues;
06368   int numentries = 0;
06369 
06370 
06371   /* Keep track of where in the table we are */
06372   BigReal currdist;
06373   int distbin;
06374   currdist = 0.0;
06375   lastdist = -1.0;
06376   distbin = 0;
06377 
06378   // Read all of the values first -- we'll interpolate later
06379         while(fgets(tableline,120,enertable) && distbin <= (int) (mynearbyint(maxdist / table_spacing) + 1)) {
06380                 if (strncmp(tableline,"#",1)==0) {
06381                         continue;
06382                 }
06383     if (strncmp(tableline,"TYPE",4)==0) {
06384       fseek(enertable, -1 * strlen(tableline), SEEK_CUR); 
06385       break;
06386     }
06387 
06388     // Read an energy line from the table
06389     int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
06390 
06391     //printf("Read an energy line: %g %g %g\n", readdist, readener, readforce);
06392     if (readcount != 3) {
06393       char err_msg[512];
06394       sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
06395       NAMD_die(err_msg);
06396     }
06397 
06398     //Sanity check the current entry
06399     if (readdist < lastdist) {
06400       NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
06401     }
06402 
06403     lastdist = readdist;
06404     dists.push_back(readdist);
06405     enervalues.push_back(readener);
06406     forcevalues.push_back(readforce);
06407     numentries++;
06408   }
06409 
06410   // Check the spacing and make sure it is uniform
06411   if (dists[0] != 0.0) {
06412     NAMD_die("ERROR: First data point for energy table must be at r=0\n");
06413   }
06414   spacing = dists[1] - dists[0];
06415   for (i=2; i<(numentries - 1); i++) {
06416     BigReal myspacing;
06417     myspacing = dists[i] - dists[i-1];
06418     if (fabs(myspacing - spacing) > 0.00001) {
06419       printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
06420       NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
06421     }
06422   }
06423 
06424 // Perform cubic spline interpolation to get the energies and forces
06425 
06426   /* allocate spline coefficient matrix */
06427   // xe and xf / be and bf for energies and forces, respectively
06428   double* m = new double[numentries*numentries];
06429   double* xe = new double[numentries];
06430   double* xf = new double[numentries];
06431   double* be = new double[numentries];
06432   double* bf = new double[numentries];
06433 
06434   be[0] = 3 * (enervalues[1] - enervalues[0]);
06435   for (i=1; i < (numentries - 1); i++) {
06436 //    printf("Control point %i at %f\n", i, enervalues[i]);
06437     be[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
06438 //    printf("be is %f\n", be[i]);
06439   }
06440   be[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
06441 
06442   bf[0] = 3 * (forcevalues[1] - forcevalues[0]);
06443   for (i=1; i < (numentries - 1); i++) {
06444 //    printf("Control point %i at %f\n", i, forcevalues[i]);
06445     bf[i] = 3 * (forcevalues[i+1] - forcevalues[i-1]);
06446 //    printf("bf is %f\n", bf[i]);
06447   }
06448   bf[numentries - 1] = 3 * (forcevalues[numentries - 1] - forcevalues[numentries - 2]);
06449 
06450   memset(m,0,numentries*numentries*sizeof(double));
06451 
06452   /* initialize spline coefficient matrix */
06453   m[0] = 2;
06454   for (i = 1;  i < numentries;  i++) {
06455     m[INDEX(numentries,i-1,i)] = 1;
06456     m[INDEX(numentries,i,i-1)] = 1;
06457     m[INDEX(numentries,i,i)] = 4;
06458   }
06459   m[INDEX(numentries,numentries-1,numentries-1)] = 2;
06460 
06461   /* Now we need to solve the equation M X = b for X */
06462   // Do this simultaneously for energy and force -- ONLY because we use the same matrix
06463 
06464   //Put m in upper triangular form and apply corresponding operations to b
06465   for (i=0; i<numentries; i++) {
06466     // zero the ith column in all rows below i
06467     const BigReal baseval = m[INDEX(numentries,i,i)];
06468     for (j=i+1; j<numentries; j++) {
06469       const BigReal myval = m[INDEX(numentries,j,i)];
06470       const BigReal factor = myval / baseval;
06471 
06472       for (int k=i; k<numentries; k++) {
06473         const BigReal subval = m[INDEX(numentries,i,k)];
06474         m[INDEX(numentries,j,k)] -= (factor * subval);
06475       }
06476 
06477       be[j] -= (factor * be[i]);
06478       bf[j] -= (factor * bf[i]);
06479 
06480     }
06481   }
06482 
06483   //Now work our way diagonally up from the bottom right to assign values to X
06484   for (i=numentries-1; i>=0; i--) {
06485 
06486     //Subtract the effects of previous columns
06487     for (j=i+1; j<numentries; j++) {
06488       be[i] -= ( m[INDEX(numentries,i,j)] * xe[j] );
06489       bf[i] -= ( m[INDEX(numentries,i,j)] * xf[j] );
06490     }
06491 
06492     xe[i] = be[i] / m[INDEX(numentries,i,i)];
06493     xf[i] = bf[i] / m[INDEX(numentries,i,i)];
06494 
06495   }
06496 
06497   // We now have the coefficient information we need to make the table
06498   // Now interpolate on each interval we want
06499 
06500   distbin = 0;
06501   int entriesperseg = (int) ceil(spacing / table_spacing);
06502   int table_prefix = 2 * typeindex * (int) (mynearbyint(maxdist / table_spacing) + 1);
06503 
06504   for (i=0; i<numentries-1; i++) {
06505     BigReal Ae,Be,Ce,De;
06506     BigReal Af,Bf,Cf,Df;
06507     currdist = dists[i];
06508 
06509 //    printf("Interpolating on interval %i\n", i);
06510 
06511     // Set up the coefficients for this section
06512     Ae = enervalues[i];
06513     Be = xe[i];
06514     Ce = 3 * (enervalues[i+1] - enervalues[i]) - (2 * xe[i]) - xe[i+1];
06515     De = 2 * (enervalues[i] - enervalues[i+1]) + xe[i] + xe[i+1];
06516 
06517     Af = forcevalues[i];
06518     Bf = xf[i];
06519     Cf = 3 * (forcevalues[i+1] - forcevalues[i]) - (2 * xf[i]) - xf[i+1];
06520     Df = 2 * (forcevalues[i] - forcevalues[i+1]) + xf[i] + xf[i+1];
06521 
06522     // Go over the region of interest and fill in the table
06523     for (j=0; j<entriesperseg; j++) {
06524       const BigReal mydist = currdist + (j * table_spacing);
06525       const BigReal mydistfrac = (float) j / (entriesperseg - 1);
06526       distbin = (int) mynearbyint(mydist / table_spacing);
06527       if (distbin >= (int) mynearbyint(maxdist / table_spacing)) break;
06528       BigReal energy;
06529       BigReal force;
06530 
06531       energy = Ae + (Be * mydistfrac) + (Ce * mydistfrac * mydistfrac) + (De * mydistfrac * mydistfrac * mydistfrac);
06532       force = Af + (Bf * mydistfrac) + (Cf * mydistfrac * mydistfrac) + (Df * mydistfrac * mydistfrac * mydistfrac);
06533 
06534 //      printf("Adding energy/force entry %f / %f in bins %i / %i for distance %f (%f)\n", energy, force, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1), mydist, mydistfrac);
06535       table_ener[table_prefix + 2 * distbin] = energy;
06536       table_ener[table_prefix + 2 * distbin + 1] = force;
06537       distbin++;
06538     }
06539     if (currdist >= maxdist) break;
06540   }
06541 
06542   //The procedure above leaves out the last entry -- add it explicitly
06543   distbin = (int) mynearbyint(maxdist / table_spacing);
06544 //  printf("Adding energy/force entry %f / %f in bins %i / %i\n", enervalues[numentries - 1], 0.0, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1));
06545   table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
06546   table_ener[table_prefix + 2 * distbin + 1] = 0.0;
06547   distbin++;
06548 
06549 
06550   // Clean up and make sure everything worked ok
06551   delete m;
06552   delete xe;
06553   delete xf;
06554   delete be;
06555   delete bf;
06556   distbin--;
06557   printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (mynearbyint(maxdist / table_spacing)), maxdist, table_spacing);
06558   if (distbin != (int) (mynearbyint(maxdist / table_spacing))) return 1;
06559   return 0;
06560 } /* end read_energy_type_bothcubspline */
06561 
06562 /**************************************************************************
06563  * FUNCTION read_energy_type_cubspline
06564  *
06565  * Read a single type block from an energy table file, using cubic spline interpolation
06566  *
06567  * Inputs:
06568  *  enertable - File stream positioned at the start of the type block
06569  *  typeindex  - integer index of current type
06570  *  table_ener - pointer to array to be filled with table entries
06571  *  table_spacing - Spacing between table points (A)
06572  *  maxdist - Longest distance needed in table
06573  *
06574  * Return values:
06575  *  0 on normal exit
06576  *  1 if not enough entries were present to fill out the table
06577  *
06578  *  Note: enertable should be left positioned immediately BEFORE the next
06579  *  TYPE block starts
06580  *  **********************************************************************/
06581 
06582 int Parameters::read_energy_type_cubspline(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
06583 
06584   char tableline[120];
06585   int i,j;
06586 
06587   /* Last values read from table */
06588   BigReal readdist;
06589   BigReal readener;
06590   BigReal spacing;
06591 //  BigReal readforce;
06592   BigReal lastdist;
06593 //  BigReal lastener;
06594 //  BigReal lastforce;
06595 //  readdist = -1.0;
06596 //  readener = 0.0;
06597 //  readforce = 0.0;
06598 
06599   /* Create arrays for holding the input values */
06600   std::vector<BigReal>  dists;
06601   std::vector<BigReal> enervalues;
06602   int numentries = 0;
06603 
06604 
06605   /* Keep track of where in the table we are */
06606   BigReal currdist;
06607   int distbin;
06608   currdist = 0.0;
06609   lastdist = -1.0;
06610   distbin = 0;
06611 
06612   // Read all of the values first -- we'll interpolate later
06613         while(fgets(tableline,120,enertable) && distbin <= (int) (mynearbyint(maxdist / table_spacing) + 1)) {
06614                 if (strncmp(tableline,"#",1)==0) {
06615                         continue;
06616                 }
06617     if (strncmp(tableline,"TYPE",4)==0) {
06618       fseek(enertable, -1 * strlen(tableline), SEEK_CUR); 
06619       break;
06620     }
06621 
06622     // Read an energy line from the table
06623     int readcount = sscanf(tableline, "%lf %lf", &readdist, &readener);
06624 
06625    // printf("Read an energy line: %f %f %f\n", readdist, readener, readforce);
06626     if (readcount != 2) {
06627       char err_msg[512];
06628       sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
06629       NAMD_die(err_msg);
06630     }
06631 
06632     //Sanity check the current entry
06633     if (readdist < lastdist) {
06634       NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
06635     }
06636 
06637     lastdist = readdist;
06638     dists.push_back(readdist);
06639     enervalues.push_back(readener);
06640     numentries++;
06641   }
06642 
06643   // Check the spacing and make sure it is uniform
06644   if (dists[0] != 0.0) {
06645     NAMD_die("ERROR: First data point for energy table must be at r=0\n");
06646   }
06647   spacing = dists[1] - dists[0];
06648   for (i=2; i<(numentries - 1); i++) {
06649     BigReal myspacing;
06650     myspacing = dists[i] - dists[i-1];
06651     if (fabs(myspacing - spacing) > 0.00001) {
06652       printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
06653       NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
06654     }
06655   }
06656 
06657 // Perform cubic spline interpolation to get the energies and forces
06658 
06659   /* allocate spline coefficient matrix */
06660   double* m = new double[numentries*numentries];
06661   double* x = new double[numentries];
06662   double* b = new double[numentries];
06663 
06664   b[0] = 3 * (enervalues[1] - enervalues[0]);
06665   for (i=1; i < (numentries - 1); i++) {
06666     printf("Control point %i at %f\n", i, enervalues[i]);
06667     b[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
06668     printf("b is %f\n", b[i]);
06669   }
06670   b[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
06671 
06672   memset(m,0,numentries*numentries*sizeof(double));
06673 
06674   /* initialize spline coefficient matrix */
06675   m[0] = 2;
06676   for (i = 1;  i < numentries;  i++) {
06677     m[INDEX(numentries,i-1,i)] = 1;
06678     m[INDEX(numentries,i,i-1)] = 1;
06679     m[INDEX(numentries,i,i)] = 4;
06680   }
06681   m[INDEX(numentries,numentries-1,numentries-1)] = 2;
06682 
06683   /* Now we need to solve the equation M X = b for X */
06684 
06685   printf("Solving the matrix equation: \n");
06686 
06687   for (i=0; i<numentries; i++) {
06688     printf(" ( ");
06689     for (j=0; j<numentries; j++) {
06690       printf(" %6.3f,", m[INDEX(numentries, i, j)]);
06691     }
06692     printf(" )  ( D%-3i )  =  ( %6.3f )\n", i, b[i]);
06693   }
06694 
06695   //Put m in upper triangular form and apply corresponding operations to b
06696   for (i=0; i<numentries; i++) {
06697     // zero the ith column in all rows below i
06698     const BigReal baseval = m[INDEX(numentries,i,i)];
06699     for (j=i+1; j<numentries; j++) {
06700       const BigReal myval = m[INDEX(numentries,j,i)];
06701       const BigReal factor = myval / baseval;
06702 
06703       for (int k=i; k<numentries; k++) {
06704         const BigReal subval = m[INDEX(numentries,i,k)];
06705         m[INDEX(numentries,j,k)] -= (factor * subval);
06706       }
06707 
06708       b[j] -= (factor * b[i]);
06709 
06710     }
06711   }
06712 
06713   printf(" In upper diagonal form, equation is:\n");
06714   for (i=0; i<numentries; i++) {
06715     printf(" ( ");
06716     for (j=0; j<numentries; j++) {
06717       printf(" %6.3f,", m[INDEX(numentries, i, j)]);
06718     }
06719     printf(" )  ( D%-3i )  =  ( %6.3f )\n", i, b[i]);
06720   }
06721 
06722   //Now work our way diagonally up from the bottom right to assign values to X
06723   for (i=numentries-1; i>=0; i--) {
06724 
06725     //Subtract the effects of previous columns
06726     for (j=i+1; j<numentries; j++) {
06727       b[i] -= ( m[INDEX(numentries,i,j)] * x[j] );
06728     }
06729 
06730     x[i] = b[i] / m[INDEX(numentries,i,i)];
06731 
06732   }
06733 
06734   printf(" Solution vector is:\n\t(");
06735   for (i=0; i<numentries; i++) {
06736     printf(" %6.3f ", x[i]);
06737   }
06738   printf(" ) \n");
06739 
06740   // We now have the coefficient information we need to make the table
06741   // Now interpolate on each interval we want
06742 
06743   distbin = 0;
06744   int entriesperseg = (int) ceil(spacing / table_spacing);
06745   int table_prefix = 2 * typeindex * (int) (mynearbyint(maxdist / table_spacing) + 1);
06746 
06747   for (i=0; i<numentries-1; i++) {
06748     BigReal A,B,C,D;
06749     currdist = dists[i];
06750 
06751     printf("Interpolating on interval %i\n", i);
06752 
06753     // Set up the coefficients for this section
06754     A = enervalues[i];
06755     B = x[i];
06756     C = 3 * (enervalues[i+1] - enervalues[i]) - (2 * x[i]) - x[i+1];
06757     D = 2 * (enervalues[i] - enervalues[i+1]) + x[i] + x[i+1];
06758 
06759     printf("Coefficients for this interval: %f %f %f %f\n", A, B, C, D);
06760 
06761     // Go over the region of interest and fill in the table
06762     for (j=0; j<entriesperseg; j++) {
06763       const BigReal mydist = currdist + (j * table_spacing);
06764       const BigReal mydistfrac = (float) j / (entriesperseg - 1);
06765       distbin = (int) mynearbyint(mydist / table_spacing);
06766       if (distbin >= (int) mynearbyint(maxdist / table_spacing)) break;
06767       BigReal energy;
06768       BigReal force;
06769 
06770       energy = A + (B * mydistfrac) + (C * mydistfrac * mydistfrac) + (D * mydistfrac * mydistfrac * mydistfrac);
06771       force = B + (2 * C * mydistfrac) + (3 * D * mydistfrac * mydistfrac);
06772       // Multiply force by 1 / (interval length)
06773       force *= (1.0 / spacing);
06774 
06775       printf("Adding energy/force entry %f / %f in bins %i / %i for distance %f (%f)\n", energy, force, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1), mydist, mydistfrac);
06776       table_ener[table_prefix + 2 * distbin] = energy;
06777       table_ener[table_prefix + 2 * distbin + 1] = force;
06778       distbin++;
06779     }
06780     if (currdist >= maxdist) break;
06781   }
06782 
06783   //The procedure above leaves out the last entry -- add it explicitly
06784   distbin = (int) mynearbyint(maxdist / table_spacing);
06785   printf("Adding energy/force entry %f / %f in bins %i / %i\n", enervalues[numentries - 1], 0.0, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1));
06786   table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
06787   table_ener[table_prefix + 2 * distbin + 1] = 0.0;
06788   distbin++;
06789 
06790 
06791   // Clean up and make sure everything worked ok
06792   delete m;
06793   delete x;
06794   delete b;
06795   distbin--;
06796   printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (mynearbyint(maxdist / table_spacing)), maxdist, table_spacing);
06797   if (distbin != (int) (mynearbyint(maxdist / table_spacing))) return 1;
06798   return 0;
06799 } /* end read_energy_type_cubspline */
06800 
06801 /**************************************************************************
06802  * FUNCTION read_energy_type
06803  *
06804  * Read a single type block from an energy table file
06805  *
06806  * Inputs:
06807  *  enertable - File stream positioned at the start of the type block
06808  *  typeindex  - integer index of current type
06809  *  table_ener - pointer to array to be filled with table entries
06810  *  table_spacing - Spacing between table points (A)
06811  *  maxdist - Longest distance needed in table
06812  *
06813  * Return values:
06814  *  0 on normal exit
06815  *  1 if not enough entries were present to fill out the table
06816  *
06817  *  Note: enertable should be left positioned immediately BEFORE the next
06818  *  TYPE block starts
06819  *  **********************************************************************/
06820 
06821 int Parameters::read_energy_type(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
06822 
06823   char tableline[120];
06824 
06825   /* Last values read from table */
06826   BigReal readdist;
06827   BigReal readener;
06828   BigReal readforce;
06829   BigReal lastdist;
06830   BigReal lastener;
06831   BigReal lastforce;
06832   readdist = -1.0;
06833   readener = 0.0;
06834   readforce = 0.0;
06835 
06836   /* Keep track of where in the table we are */
06837   float currdist;
06838   int distbin;
06839   currdist = -1.0;
06840   distbin = -1;
06841 
06842         while(fgets(tableline,120,enertable) && distbin <= (int) (mynearbyint(maxdist / table_spacing) + 1)) {
06843     printf("At distance %f + %f vs. %f\n", currdist, table_spacing, maxdist);
06844                 if (strncmp(tableline,"#",1)==0) {
06845                         continue;
06846                 }
06847     if (strncmp(tableline,"TYPE",4)==0) {
06848       break;
06849     }
06850 
06851     // Read an energy line from the table
06852     lastdist = readdist;
06853     lastener = readener;
06854     lastforce = readforce;
06855