Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

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

Generated on Fri May 25 04:07:16 2012 for NAMD by  doxygen 1.3.9.1