Parameters.C

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

Generated on Tue Sep 19 01:17:13 2017 for NAMD by  doxygen 1.4.7