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 #ifndef WIN32
00019 #include <strings.h>
00020 #endif
00021 #include "InfoStream.h"
00022 #include <charm++.h>
00023 #include "Parameters.h"
00024 #include "Communicate.h"
00025 #include "ConfigList.h"
00026 //****** BEGIN CHARMM/XPLOR type changes
00027 #include "SimParameters.h"
00028 //****** END CHARMM/XPLOR type changes
00029 
00030 #define MIN_DEBUG_LEVEL 3
00031 //#define DEBUGM
00032 #include "Debug.h"
00033 
00034 
00035 //  struct bond_params is used to form a binary tree of bond parameters.
00036 //  The two atom names are used to determine the order of the nodes in the
00037 //  tree.  atom1name should ALWAYS be lexically before atom2name
00038 
00039 struct bond_params
00040 {
00041   char atom1name[11];
00042   char atom2name[11];
00043   Real forceconstant;
00044   Real distance;
00045   Index index;
00046   struct bond_params *left;
00047   struct bond_params *right;
00048 };
00049 
00050 //  struct angle_params is used to form a binary tree of bond parameters.
00051 //  The three atom names are used to determine the order of the nodes in
00052 //  the tree.  atom1name should ALWAYS be lexically before atom3name
00053 
00054 struct angle_params
00055 {
00056   char atom1name[11];
00057   char atom2name[11];
00058   char atom3name[11];
00059   Real forceconstant;
00060   int normal;
00061   Real angle;
00062   Real k_ub;
00063   Real r_ub;
00064   Index index;
00065   struct angle_params *left;
00066   struct angle_params *right;
00067 };
00068 
00069 //  struct dihedral_params is used to form a linked list of the dihedral
00070 //  parameters.  The linked list is arranged in such a way that any
00071 //  bonds with wildcards are at the end of the list so that a linear
00072 //  search can be done but we will still find exact matches before
00073 //  wildcard matches
00074 
00075 struct dihedral_params
00076 {
00077   char atom1name[11];
00078   char atom2name[11];
00079   char atom3name[11];
00080   char atom4name[11];
00081   char atom1wild;
00082   char atom2wild;
00083   char atom3wild;
00084   char atom4wild;
00085   int multiplicity;
00086   FourBodyConsts values[MAX_MULTIPLICITY];
00087   Index index;
00088   struct dihedral_params *next;
00089   dihedral_params() { memset(this, 0, sizeof(dihedral_params)); }
00090 };
00091 
00092 //  struct improper_params is used to form a linked list of the improper
00093 //  parameters.  The linked list is arranged in such a way that any
00094 //  bonds with wildcards are at the end of the list so that a linear
00095 //  search can be done but we will still find exact matches before
00096 //  wildcard matches
00097 
00098 struct improper_params
00099 {
00100   char atom1name[11];
00101   char atom2name[11];
00102   char atom3name[11];
00103   char atom4name[11];
00104   int multiplicity;
00105   FourBodyConsts values[MAX_MULTIPLICITY];
00106   Index index;
00107   struct improper_params *next;
00108 };
00109 
00110 struct crossterm_params
00111 {
00112   crossterm_params(int dim) : dimension(dim) {
00113     values = new double[dimension*dimension];
00114   }
00115   ~crossterm_params() {
00116     delete [] values;
00117   }
00118   char atom1name[11];
00119   char atom2name[11];
00120   char atom3name[11];
00121   char atom4name[11];
00122   char atom5name[11];
00123   char atom6name[11];
00124   char atom7name[11];
00125   char atom8name[11];
00126   int dimension;  // usually 24
00127   double *values;  // dimension * dimension data
00128   Index index;
00129   struct crossterm_params *next;
00130 };
00131 
00132 //  struct vdw_params is used to form a binary serach tree of the
00133 //  vdw paramters for a single atom.
00134 
00135 struct vdw_params
00136 {
00137   char atomname[11];
00138   Real sigma;
00139   Real epsilon;
00140   Real sigma14;
00141   Real epsilon14;
00142   Index index;
00143   struct vdw_params *left;
00144   struct vdw_params *right;
00145 };
00146 
00147 //  struct vdw_pair_params is used to form a linked list of the
00148 //  vdw parameters for a pair of atoms
00149 
00150 struct vdw_pair_params
00151 {
00152   char atom1name[11];
00153   char atom2name[11];
00154   Real A;
00155   Real B;
00156   Real A14;
00157   Real B14;
00158   struct vdw_pair_params *next;
00159 };
00160 
00161 
00162 Parameters::Parameters() {
00163   initialize();
00164 }
00165 
00166 void Parameters::initialize() {
00167 
00168   paramType = -1;
00169 
00170   /*  Set all the pointers to NULL        */
00171   atomTypeNames=NULL;
00172   bondp=NULL;
00173   anglep=NULL;
00174   improperp=NULL;
00175   dihedralp=NULL;
00176   crosstermp=NULL;
00177   vdwp=NULL;
00178   vdw_pairp=NULL;
00179   bond_array=NULL;
00180   angle_array=NULL;
00181   dihedral_array=NULL;
00182   improper_array=NULL;
00183   crossterm_array=NULL;
00184   vdw_array=NULL;
00185   vdw_pair_tree=NULL;
00186   maxDihedralMults=NULL;
00187   maxImproperMults=NULL;
00188 
00189   /*  Set all the counts to 0          */
00190   NumBondParams=0;
00191   NumAngleParams=0;
00192   NumDihedralParams=0;
00193   NumImproperParams=0;
00194   NumCrosstermParams=0;
00195   NumVdwParams=0;
00196   NumVdwPairParams=0;
00197   NumCosAngles=0;
00198 }
00199 
00200 /************************************************************************/
00201 /*                  */
00202 /*      FUNCTION Parameters        */
00203 /*                  */
00204 /*  This is the constructor for the class.  It simlpy sets the      */
00205 /*  pointers to the list and trees to NULL and the count of all the     */
00206 /*  parameters to 0.              */
00207 /*  The type (format) of the input parameters (Xplor,Charmm) is set here. */
00208 /*                  */
00209 /************************************************************************/
00210 
00211 Parameters::Parameters(SimParameters *simParams, StringList *f)
00212 {
00213   initialize();
00214 
00215   //****** BEGIN CHARMM/XPLOR type changes
00217   if (simParams->paraTypeXplorOn)
00218   {
00219     paramType = paraXplor;
00220   }
00221   else if (simParams->paraTypeCharmmOn)
00222   {
00223     paramType = paraCharmm;
00224   }
00225   //****** END CHARMM/XPLOR type changes
00226   //Test for cos-based angles
00227   if (simParams->cosAngles) {
00228     cosAngles = true;
00229   } else {
00230     cosAngles = false;
00231   }
00232 
00233   /* Set up AllFilesRead flag to FALSE.  Once all of the files    */
00234   /* have been read in, then this will be set to true and the     */
00235   /* arrays of parameters will be set up        */
00236   AllFilesRead = FALSE;
00237 
00238   if (NULL != f) 
00239   {
00240     do
00241     {
00242       //****** BEGIN CHARMM/XPLOR type changes
00243       if (paramType == paraXplor)
00244       {
00245         read_parameter_file(f->data);
00246       }
00247       else if (paramType == paraCharmm)
00248       {
00249         read_charmm_parameter_file(f->data);
00250       }
00251       //****** END CHARMM/XPLOR type changes
00252       f = f->next;
00253     } while ( f != NULL );
00254 
00255     done_reading_files();
00256   }
00257 
00258 }
00259 /*      END OF FUNCTION Parameters      */
00260 
00261 /************************************************************************/
00262 /*                  */
00263 /*      FUNCTION ~Parameters        */
00264 /*                  */
00265 /*  This is the destructor for this class.  It basically just       */
00266 /*  frees all of the memory allocated for the parameters.    */
00267 /*                  */
00268 /************************************************************************/
00269 
00270 Parameters::~Parameters()
00271 
00272 {
00273         if (atomTypeNames)
00274           delete [] atomTypeNames;
00275 
00276   if (bondp != NULL)
00277     free_bond_tree(bondp);
00278 
00279   if (anglep != NULL)
00280     free_angle_tree(anglep);
00281 
00282   if (dihedralp != NULL)
00283     free_dihedral_list(dihedralp);
00284 
00285   if (improperp != NULL)
00286     free_improper_list(improperp);
00287 
00288   if (crosstermp != NULL)
00289     free_crossterm_list(crosstermp);
00290 
00291   if (vdwp != NULL)
00292     free_vdw_tree(vdwp);
00293 
00294   if (vdw_pairp != NULL)
00295     free_vdw_pair_list();
00296 
00297   if (bond_array != NULL)
00298     delete [] bond_array;
00299 
00300   if (angle_array != NULL)
00301     delete [] angle_array;
00302 
00303   if (dihedral_array != NULL)
00304     delete [] dihedral_array;
00305 
00306   if (improper_array != NULL)
00307     delete [] improper_array;
00308 
00309   if (crossterm_array != NULL)
00310     delete [] crossterm_array;
00311 
00312   if (vdw_array != NULL)
00313     delete [] vdw_array;
00314   
00315   if (vdw_pair_tree != NULL)
00316     free_vdw_pair_tree(vdw_pair_tree);
00317 
00318   if (maxDihedralMults != NULL)
00319     delete [] maxDihedralMults;
00320 
00321   if (maxImproperMults != NULL)
00322     delete [] maxImproperMults;
00323 
00324   for( int i = 0; i < error_msgs.size(); ++i ) {
00325     delete [] error_msgs[i];
00326   }
00327   error_msgs.resize(0);
00328 }
00329 /*      END OF FUNCTION ~Parameters      */
00330 
00331 /************************************************************************/
00332 /*                  */
00333 /*      FUNCTION read_paramter_file      */
00334 /*                  */
00335 /*   INPUTS:                */
00336 /*  fname - name of the parameter file to read      */
00337 /*                  */
00338 /*  This function reads in a parameter file and adds the parameters */
00339 /*   from this file to the current group of parameters.  The basic      */
00340 /*   logic of the routine is to read in a line from the file, looks at  */
00341 /*   the first word of the line to determine what kind of parameter we  */
00342 /*   have, and then call the appropriate routine to add the parameter   */
00343 /*   to the parameters that we have.          */
00344 /*                  */
00345 /************************************************************************/
00346 
00347 void Parameters::read_parameter_file(char *fname)
00348 
00349 {
00350   char buffer[512];  //  Buffer to store each line of the file
00351   char first_word[512];  //  First word of the current line
00352   FILE *pfile;    //  File descriptor for the parameter file
00353 
00354   /*  Check to make sure that we haven't previously been told     */
00355   /*  that all the files were read        */
00356   if (AllFilesRead)
00357   {
00358     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00359   }
00360 
00361   /*  Try and open the file          */
00362   if ( (pfile = Fopen(fname, "r")) == NULL)
00363   {
00364     char err_msg[256];
00365 
00366     sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
00367     NAMD_die(err_msg);
00368   }
00369 
00370   /*  Keep reading in lines until we hit the EOF      */
00371   while (NAMD_read_line(pfile, buffer) != -1)
00372   {
00373     /*  Get the first word of the line      */
00374     NAMD_find_first_word(buffer, first_word);
00375 
00376     /*  First, screen out things that we ignore, such as    */
00377     /*  blank lines, lines that start with '!', lines that  */
00378     /*  start with "REMARK", lines that start with set",    */
00379     /*  and most of the HBOND parameters which include      */
00380     /*  AEXP, REXP, HAEX, AAEX, but not the HBOND statement */
00381     /*  which is parsed.                                    */
00382     if ((buffer[0] != '!') && 
00383         !NAMD_blank_string(buffer) &&
00384         (strncasecmp(first_word, "REMARK", 6) != 0) &&
00385         (strcasecmp(first_word, "set")!=0) &&
00386         (strncasecmp(first_word, "AEXP", 4) != 0) &&
00387         (strncasecmp(first_word, "REXP", 4) != 0) &&
00388         (strncasecmp(first_word, "HAEX", 4) != 0) &&
00389         (strncasecmp(first_word, "AAEX", 4) != 0) &&
00390         (strncasecmp(first_word, "NBOND", 5) != 0) &&
00391         (strncasecmp(first_word, "CUTNB", 5) != 0) &&
00392         (strncasecmp(first_word, "END", 3) != 0) &&
00393         (strncasecmp(first_word, "CTONN", 5) != 0) &&
00394         (strncasecmp(first_word, "EPS", 3) != 0) &&
00395         (strncasecmp(first_word, "VSWI", 4) != 0) &&
00396         (strncasecmp(first_word, "NBXM", 4) != 0) &&
00397         (strncasecmp(first_word, "INHI", 4) != 0) )
00398     {
00399       /*  Now, call the appropriate function based    */
00400       /*  on the type of parameter we have    */
00401       if (strncasecmp(first_word, "bond", 4)==0)
00402       {
00403         add_bond_param(buffer);
00404         NumBondParams++;
00405       }
00406       else if (strncasecmp(first_word, "angl", 4)==0)
00407       {
00408         add_angle_param(buffer);
00409         NumAngleParams++;
00410       }
00411       else if (strncasecmp(first_word, "dihe", 4)==0)
00412       {
00413         add_dihedral_param(buffer, pfile);
00414         NumDihedralParams++;
00415       }
00416       else if (strncasecmp(first_word, "impr", 4)==0)
00417       {
00418         add_improper_param(buffer, pfile);
00419         NumImproperParams++;
00420       }
00421       else if (strncasecmp(first_word, "nonb", 4)==0)
00422       {
00423         add_vdw_param(buffer);
00424         NumVdwParams++; 
00425       }
00426       else if (strncasecmp(first_word, "nbfi", 4)==0)
00427       {
00428         add_vdw_pair_param(buffer);
00429         NumVdwPairParams++; 
00430       }
00431       else if (strncasecmp(first_word, "hbon", 4)==0)
00432       {
00433         add_hb_pair_param(buffer);
00434       }
00435       else
00436       {
00437         /*  This is an unknown paramter.        */
00438         /*  This is BAD        */
00439         char err_msg[512];
00440 
00441         sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
00442            fname, buffer);
00443         NAMD_die(err_msg);
00444       }
00445     }
00446   }
00447 
00448   /*  Close the file            */
00449   Fclose(pfile);
00450 
00451   return;
00452 }
00453 /*      END OF FUNCTION read_paramter_file    */
00454 
00455 //****** BEGIN CHARMM/XPLOR type changes
00456 /************************************************************************/
00457 /*                                                                        */
00458 /*                        FUNCTION read_charmm_paramter_file                */
00459 /*                                                                        */
00460 /*   INPUTS:                                                                */
00461 /*        fname - name of the parameter file to read                        */
00462 /*                                                                        */
00463 /*        This function reads in a CAHRMM parameter file and adds the     */ 
00464 /*   parameters from this file to the current group of parameters.      */
00465 /*   The basic logic of the routine is to first find out what type of   */
00466 /*   parameter we have in the file. Then look at each line in turn      */
00467 /*   and call the appropriate routine to add the parameters until we hit*/
00468 /*   a new type of parameter or EOF.                                    */
00469 /*                                                                        */
00470 /************************************************************************/
00471 
00472 void Parameters::read_charmm_parameter_file(char *fname)
00473 
00474 {
00475   int  par_type=0;         //  What type of parameter are we currently
00476                            //  dealing with? (vide infra)
00477   int  skipline;           //  skip this line?
00478   int  skipall = 0;        //  skip rest of file;
00479   char buffer[512];           //  Buffer to store each line of the file
00480   char first_word[512];           //  First word of the current line
00481   FILE *pfile;                   //  File descriptor for the parameter file
00482 
00483   /*  Check to make sure that we haven't previously been told     */
00484   /*  that all the files were read                                */
00485   if (AllFilesRead)
00486   {
00487     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00488   }
00489 
00490   /*  Try and open the file                                        */
00491   if ( (pfile = fopen(fname, "r")) == NULL)
00492   {
00493     char err_msg[256];
00494 
00495     sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
00496     NAMD_die(err_msg);
00497   }
00498 
00499   /*  Keep reading in lines until we hit the EOF                        */
00500   while (NAMD_read_line(pfile, buffer) != -1)
00501   {
00502     /*  Get the first word of the line                        */
00503     NAMD_find_first_word(buffer, first_word);
00504     skipline=0;
00505 
00506     /*  First, screen out things that we ignore.                   */   
00507     /*  blank lines, lines that start with '!' or '*', lines that  */
00508     /*  start with "END".                                          */
00509     if (!NAMD_blank_string(buffer) &&
00510         (strncmp(first_word, "!", 1) != 0) &&
00511          (strncmp(first_word, "*", 1) != 0) &&
00512         (strncasecmp(first_word, "END", 3) != 0))
00513     {
00514       if ( skipall ) {
00515         iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
00516         break;
00517       }
00518       /*  Now, determine the apropriate parameter type.   */
00519       if (strncasecmp(first_word, "bond", 4)==0)
00520       {
00521         par_type=1; skipline=1;
00522       }
00523       else if (strncasecmp(first_word, "angl", 4)==0)
00524       {
00525         par_type=2; skipline=1;
00526       }
00527       else if (strncasecmp(first_word, "thet", 4)==0)
00528       {
00529         par_type=2; skipline=1;
00530       }
00531       else if (strncasecmp(first_word, "dihe", 4)==0)
00532       {
00533         par_type=3; skipline=1;
00534       }
00535       else if (strncasecmp(first_word, "phi", 3)==0)
00536       {
00537         par_type=3; skipline=1;
00538       }
00539       else if (strncasecmp(first_word, "impr", 4)==0)
00540       {
00541         par_type=4; skipline=1;
00542       }
00543       else if (strncasecmp(first_word, "imph", 4)==0)
00544       {
00545         par_type=4; skipline=1;
00546       }
00547       else if (strncasecmp(first_word, "nbon", 4)==0)
00548       {
00549         par_type=5; skipline=1;
00550       }
00551       else if (strncasecmp(first_word, "nonb", 4)==0)
00552       {
00553         par_type=5; skipline=1;
00554       }
00555       else if (strncasecmp(first_word, "nbfi", 4)==0)
00556       {
00557         par_type=6; skipline=1;
00558       }
00559       else if (strncasecmp(first_word, "hbon", 4)==0)
00560       {
00561         par_type=7; skipline=1;
00562       }
00563       else if (strncasecmp(first_word, "cmap", 4)==0)
00564       {
00565         par_type=8; skipline=1;
00566       }
00567       else if (strncasecmp(first_word, "read", 4)==0)
00568       {
00569         skip_stream_read(buffer, pfile);  skipline=1;
00570       }
00571       else if (strncasecmp(first_word, "return", 4)==0)
00572       {
00573         skipall=1;  skipline=1;
00574       }
00575       else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
00576                (strncasecmp(first_word, "grou", 4) == 0) ||
00577                (strncasecmp(first_word, "cdie", 4) == 0) ||
00578                (strncasecmp(first_word, "shif", 4) == 0) ||
00579                (strncasecmp(first_word, "vgro", 4) == 0) ||
00580                (strncasecmp(first_word, "vdis", 4) == 0) ||
00581                (strncasecmp(first_word, "vswi", 4) == 0) ||
00582                (strncasecmp(first_word, "cutn", 4) == 0) ||
00583                (strncasecmp(first_word, "ctof", 4) == 0) ||
00584                (strncasecmp(first_word, "cton", 4) == 0) ||
00585                (strncasecmp(first_word, "eps", 3) == 0) ||
00586                (strncasecmp(first_word, "e14f", 4) == 0) ||
00587                (strncasecmp(first_word, "wmin", 4) == 0) ||
00588                (strncasecmp(first_word, "aexp", 4) == 0) ||
00589                (strncasecmp(first_word, "rexp", 4) == 0) ||
00590                (strncasecmp(first_word, "haex", 4) == 0) ||
00591                (strncasecmp(first_word, "aaex", 4) == 0) ||
00592                (strncasecmp(first_word, "noac", 4) == 0) ||
00593                (strncasecmp(first_word, "hbno", 4) == 0) ||
00594                (strncasecmp(first_word, "cuth", 4) == 0) ||
00595                (strncasecmp(first_word, "ctof", 4) == 0) ||
00596                (strncasecmp(first_word, "cton", 4) == 0) ||
00597                (strncasecmp(first_word, "cuth", 4) == 0) ||
00598                (strncasecmp(first_word, "ctof", 4) == 0) ||
00599                (strncasecmp(first_word, "cton", 4) == 0) ) 
00600       {
00601         if ((par_type != 5) && (par_type != 6) && (par_type != 7))
00602         {
00603           char err_msg[512];
00604 
00605           sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00606           NAMD_die(err_msg);
00607         }
00608         else 
00609         {
00610           skipline = 1;
00611         }
00612       }        
00613       else if (par_type == 0)
00614       {
00615         /*  This is an unknown paramter.        */
00616         /*  This is BAD                                */
00617         char err_msg[512];
00618 
00619         sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00620         NAMD_die(err_msg);
00621       }
00622     }
00623     else
00624     {
00625       skipline=1;
00626     }
00627 
00628     if ( (par_type != 0) && (!skipline) )
00629     {
00630       /*  Now, call the appropriate function based    */
00631       /*  on the type of parameter we have                */
00632       /*  I know, this should really be a switch ...  */
00633       if (par_type == 1)
00634       {
00635         add_bond_param(buffer);
00636         NumBondParams++;
00637       }
00638       else if (par_type == 2)
00639       {
00640         add_angle_param(buffer);
00641         NumAngleParams++;
00642       }
00643       else if (par_type == 3)
00644       {
00645         add_dihedral_param(buffer, pfile);
00646         NumDihedralParams++;
00647       }
00648       else if (par_type == 4)
00649       {
00650         add_improper_param(buffer, pfile);
00651         NumImproperParams++;
00652       }
00653       else if (par_type == 5)
00654       {
00655         add_vdw_param(buffer);
00656         NumVdwParams++;
00657       }
00658       else if (par_type == 6)
00659       {
00660         add_vdw_pair_param(buffer);
00661         NumVdwPairParams++; 
00662       }
00663       else if (par_type == 7)
00664       {
00665         add_hb_pair_param(buffer);                  
00666       }
00667       else if (par_type == 8)
00668       {
00669         add_crossterm_param(buffer, pfile);                  
00670         NumCrosstermParams++;
00671       }
00672       else
00673       {
00674         /*  This really should not occour!      */
00675         /*  This is an internal error.          */
00676         /*  This is VERY BAD                        */
00677         char err_msg[512];
00678 
00679         sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00680         NAMD_die(err_msg);
00681       }
00682     }
00683   }
00684 
00685   /*  Close the file                                                */
00686   fclose(pfile);
00687 
00688   return;
00689 }
00690 /*                        END OF FUNCTION read_charmm_paramter_file                */
00691 //****** END CHARMM/XPLOR type changes
00692 
00693 
00694 void Parameters::skip_stream_read(char *buf, FILE *fd) {
00695 
00696   char buffer[513];
00697   char first_word[513];
00698   char s1[128];
00699   char s2[128];
00700   int rval = sscanf(buf, "%s %s", s1, s2);
00701   if (rval != 2) {
00702         char err_msg[512];
00703         sprintf(err_msg, "BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
00704         NAMD_die(err_msg);
00705   }
00706   if ( ! strncasecmp(s2,"PARA",4) ) return;  // read parameters
00707   
00708   iout << iINFO << "SKIPPING " << s2 << " SECTION IN STREAM FILE\n" << endi;
00709 
00710   while (NAMD_read_line(fd, buffer) != -1)
00711   {
00712     // read until we find "END"
00713     NAMD_find_first_word(buffer, first_word);
00714     if (!NAMD_blank_string(buffer) &&
00715         (strncmp(first_word, "!", 1) != 0) &&
00716          (strncmp(first_word, "*", 1) != 0) &&
00717          (strncasecmp(first_word, "END", 3) == 0) ) return;
00718   }
00719 
00720 }
00721 
00722 
00723 /************************************************************************/
00724 /*                  */
00725 /*      FUNCTION add_bond_param        */
00726 /*                  */
00727 /*   INPUTS:                */
00728 /*  buf - Line from parameter file containing bond parameters  */
00729 /*                  */
00730 /*  This function adds a new bond paramter to the binary tree of    */
00731 /*   angle paramters that we have.  If a duplicate is found, a warning  */
00732 /*   message is printed and the new parameters are used.    */
00733 /*                  */
00734 /************************************************************************/
00735 
00736 void Parameters::add_bond_param(char *buf)
00737 
00738 {
00739   char atom1name[11];    //  Atom type for atom 1
00740   char atom2name[11];    //  Atom type for atom 2
00741   Real forceconstant;    //  Force constant for bond
00742   Real distance;      //  Rest distance for bond
00743   int read_count;      //  Count from sscanf
00744   struct bond_params *new_node;  //  New node in tree
00745 
00746   //****** BEGIN CHARMM/XPLOR type changes
00747   /*  Use sscanf to parse up the input line      */
00748   if (paramType == paraXplor)
00749   {
00750     /* read XPLOR format */
00751     read_count=sscanf(buf, "%*s %s %s %f %f\n", atom1name, atom2name, 
00752        &forceconstant, &distance);
00753   }
00754   else if (paramType == paraCharmm)
00755   {
00756     /* read CHARMM format */
00757     read_count=sscanf(buf, "%s %s %f %f\n", atom1name, atom2name, 
00758        &forceconstant, &distance);
00759   }
00760   //****** END CHARMM/XPLOR type changes
00761 
00762   /*  Check to make sure we found everything we expeceted    */
00763   if (read_count != 4)
00764   {
00765     char err_msg[512];
00766 
00767     if (paramType == paraXplor)
00768     {
00769       sprintf(err_msg, "BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
00770     }
00771     else if (paramType == paraCharmm)
00772     {
00773       sprintf(err_msg, "BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
00774     }
00775     NAMD_die(err_msg);
00776   }
00777 
00778   /*  Allocate a new node            */
00779   new_node = new bond_params;
00780 
00781   if (new_node == NULL)
00782   {
00783     NAMD_die("memory allocation failed in Parameters::add_bond_param\n");
00784   }
00785 
00786   /*  Order the atoms so that the atom that comes alphabetically  */
00787   /*  first is atom 1.  Since the bond is symmetric, it doesn't   */
00788   /*  matter physically which atom is first.  And this allows the */
00789   /*  search of the binary tree to be done in a logical manner    */
00790   if (strcasecmp(atom1name, atom2name) < 0)
00791   {
00792     strcpy(new_node->atom1name, atom1name);
00793     strcpy(new_node->atom2name, atom2name);
00794   }
00795   else
00796   {
00797     strcpy(new_node->atom2name, atom1name);
00798     strcpy(new_node->atom1name, atom2name);
00799   }
00800 
00801   /*  Assign force constant and distance        */
00802   new_node->forceconstant = forceconstant;
00803   new_node->distance = distance;
00804 
00805   /*  Set pointers to null          */
00806   new_node->left = NULL;
00807   new_node->right = NULL;
00808 
00809   /*  Make call to recursive call to actually add the node to the */
00810   /*  tree              */
00811   bondp=add_to_bond_tree(new_node, bondp);
00812 
00813   return;
00814 }
00815 /*      END OF FUNCTION add_bond_param      */
00816 
00817 /************************************************************************/
00818 /*                  */
00819 /*      FUNCTION add_to_bond_tree      */
00820 /*                  */
00821 /*   INPUTS:                */
00822 /*  new_node - Node to add to the tree        */
00823 /*  tree - tree to add the node to          */
00824 /*                  */
00825 /*   OUTPUTS:                */
00826 /*  ths function returns a pointer to the new tree with the node    */
00827 /*   added to it.  Most of the time it will be the same pointer as was  */
00828 /*   passed in, but not if the current tree is empty.      */
00829 /*                  */
00830 /*  this is a receursive function that adds a node to the binary    */
00831 /*   tree used to store bond parameters.        */
00832 /*                  */
00833 /************************************************************************/
00834 
00835 struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
00836              struct bond_params *tree)
00837 
00838 {
00839   int compare_code;  //  Results from strcasecmp
00840 
00841   /*  If the tree is currently empty, then the new tree consists  */
00842   /*  only of the new node          */
00843   if (tree == NULL)
00844     return(new_node);
00845 
00846   /*  Compare the atom1 name from the new node and the head of    */
00847   /*  the tree              */
00848   compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
00849 
00850   /*  Check to see if they are the same        */
00851   if (compare_code == 0)
00852   {
00853     /*  The atom 1 names are the same, compare atom 2  */
00854     compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
00855 
00856     /*  If atom 1 AND atom 2 are the same, we have a duplicate */
00857     if (compare_code == 0)
00858     {
00859       /*  We have a duplicate.  So print out a warning*/
00860       /*  message.  Then assign the new values to the */
00861       /*  tree and free the new_node      */
00862       //****** BEGIN CHARMM/XPLOR type changes
00863       /* we do not care about identical replacement */
00864       if ((tree->forceconstant != new_node->forceconstant) || 
00865           (tree->distance != new_node->distance))
00866       {
00867         iout << "\n" << iWARN << "DUPLICATE BOND ENTRY FOR "
00868           << new_node->atom1name << "-"
00869           << new_node->atom2name
00870           << "\nPREVIOUS VALUES  k=" << tree->forceconstant
00871           << "  x0=" << tree->distance
00872           << "\n   USING VALUES  k=" << new_node->forceconstant
00873           << "  x0=" << new_node->distance
00874           << "\n" << endi;
00875 
00876         tree->forceconstant=new_node->forceconstant;
00877         tree->distance=new_node->distance;
00878       }
00879       //****** END CHARMM/XPLOR type changes
00880 
00881       delete new_node;
00882 
00883       return(tree);
00884     }
00885   }
00886 
00887   /*  We don't have a duplicate, so if the new value is less      */
00888   /*  than the head of the tree, add it to the left child,   */
00889   /*  otherwise add it to the right child        */
00890   if (compare_code < 0)
00891   {
00892     tree->left = add_to_bond_tree(new_node, tree->left);
00893   }
00894   else
00895   {
00896     tree->right = add_to_bond_tree(new_node, tree->right);
00897   }
00898 
00899   return(tree);
00900 }
00901 /*    END OF FUNCTION add_to_bond_tree      */
00902 
00903 /************************************************************************/
00904 /*                  */
00905 /*      FUNCTION add_angle_param      */
00906 /*                  */
00907 /*   INPUTS:                */
00908 /*  buf - line from paramter file with angle parameters    */
00909 /*                  */
00910 /*  this function adds an angle parameter.  It parses up the input  */
00911 /*   line and then adds it to the binary tree used to store the angle   */
00912 /*   parameters.              */
00913 /*                  */
00914 /************************************************************************/
00915 
00916 void Parameters::add_angle_param(char *buf)
00917 
00918 {
00919   char atom1name[11];    // Type for atom 1
00920   char atom2name[11];    // Type for atom 2
00921   char atom3name[11];    // Type for atom 3
00922   char norm[4]="xxx";
00923   Real forceconstant;    // Force constant
00924   Real angle;      // Theta 0
00925   Real k_ub;      // Urey-Bradley force constant
00926   Real r_ub;      // Urey-Bradley distance
00927   int read_count;      // count from sscanf
00928   struct angle_params *new_node;  // new node in tree
00929 
00930   //****** BEGIN CHARMM/XPLOR type changes
00931   /*  parse up the input line with sscanf        */
00932   if (paramType == paraXplor)
00933   {
00934     /* read XPLOR format */
00935     read_count=sscanf(buf, "%*s %s %s %s %f %f UB %f %f\n", 
00936        atom1name, atom2name, atom3name, &forceconstant, &angle,
00937        &k_ub, &r_ub);
00938   }
00939   else if ((paramType == paraCharmm) && cosAngles) {
00940     read_count=sscanf(buf, "%s %s %s %f %f %3s %f %f\n", 
00941        atom1name, atom2name, atom3name, &forceconstant, &angle, norm,
00942        &k_ub, &r_ub);
00943 //    printf("%s\n", buf);
00944 //    printf("Data: %s %s %s %f %f %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle, norm, k_ub, r_ub);
00945   }  
00946   else if (paramType == paraCharmm)
00947   {
00948     /* read CHARMM format */
00949     read_count=sscanf(buf, "%s %s %s %f %f %f %f\n", 
00950        atom1name, atom2name, atom3name, &forceconstant, &angle,
00951        &k_ub, &r_ub);
00952 //    printf("%s\n", buf);
00953 //    printf("Data: %s %s %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle);
00954   }
00955   //****** END CHARMM/XPLOR type changes
00956 
00957   /*  Check to make sure we got what we expected      */
00958   if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
00959   {
00960     char err_msg[512];
00961 
00962     if (paramType == paraXplor)
00963     {
00964       sprintf(err_msg, "BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
00965     }
00966     else if (paramType == paraCharmm)
00967     {
00968       sprintf(err_msg, "BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
00969     }
00970     NAMD_die(err_msg);
00971   }
00972 
00973   /*  Allocate the new node          */
00974   new_node = new angle_params;
00975 
00976   if (new_node == NULL)
00977   {
00978     NAMD_die("memory allocation failed in Parameters::add_angle_param");
00979   }
00980 
00981   /*  As with the bond, we want the atom type is comes first  */
00982   /*  alphbetically first between atom 1 and atom 3 to be in      */
00983   /*  atom 1 so that we can search the tree reliably.    */
00984   if (strcasecmp(atom1name, atom3name) < 0)
00985   {
00986     strcpy(new_node->atom1name, atom1name);
00987     strcpy(new_node->atom2name, atom2name);
00988     strcpy(new_node->atom3name, atom3name);
00989   }
00990   else
00991   {
00992     strcpy(new_node->atom3name, atom1name);
00993     strcpy(new_node->atom2name, atom2name);
00994     strcpy(new_node->atom1name, atom3name);
00995   }
00996 
00997   /*  Assign the constants and pointer values      */
00998   new_node->forceconstant = forceconstant;
00999   new_node->angle = angle;
01000 
01001   if (cosAngles) {
01002     if (strcmp("cos",norm)==0) {
01003 //      iout << "Info: Using cos mode for angle " << buf << endl;
01004       NumCosAngles++;
01005       new_node->normal = 0;
01006     } else {
01007 //      iout << "Info: Using x^2 mode for angle " << buf << endl;
01008       new_node->normal = 1;
01009     }
01010   } else {
01011     new_node->normal = 1;
01012   }
01013 
01014   if (read_count == 7)
01015   {
01016     //  Urey-Bradley constants
01017     if (new_node->normal == 0) {
01018       NAMD_die("ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
01019     }
01020     new_node->k_ub = k_ub;
01021     new_node->r_ub = r_ub;
01022   }
01023   else
01024   {
01025     new_node->k_ub = 0.0;
01026     new_node->r_ub = 0.0;
01027   }
01028 
01029   new_node->left = NULL;
01030   new_node->right = NULL;
01031 
01032   /*  Insert it into the tree          */
01033   anglep = add_to_angle_tree(new_node, anglep);
01034 
01035   return;
01036 }
01037 /*      END OF FUNCTION add_angle_param      */
01038 
01039 /************************************************************************/
01040 /*                  */
01041 /*      FUNCTION add_to_angle_tree      */
01042 /*                  */
01043 /*   INPUTS:                */
01044 /*  new_node - new node to add to the angle tree      */
01045 /*  tree - tree to add the node to          */
01046 /*                  */
01047 /*   OUTPUTS:                */
01048 /*  the function returns a pointer to the new tree with the node    */
01049 /*   added.  Most of the time, this will be the same as the value passed*/
01050 /*   in, but not in the case where the tree is empty.      */
01051 /*                  */
01052 /*  this is a recursive function that adds an angle parameter  */
01053 /*   to the binary tree storing the angle parameters.  If a duplicate   */
01054 /*   is found, a warning message is printed, the current values in the  */
01055 /*   tree are replaced with the new values, and the new node is free'd  */
01056 /*                  */
01057 /************************************************************************/
01058 
01059 struct angle_params *Parameters::add_to_angle_tree(struct angle_params *new_node,
01060              struct angle_params *tree)
01061 
01062 {
01063   int compare_code;  //  Return code from strcasecmp
01064 
01065   /*  If the tree is empty, then the new_node is the tree    */
01066   if (tree == NULL)
01067     return(new_node);
01068 
01069   /*  Compare atom 1 from the new node and the head of the tree   */
01070   compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
01071 
01072   if (compare_code == 0)
01073   {
01074     /*  Atom 1 is the same, compare atom 2      */
01075     compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
01076 
01077     if (compare_code == 0)
01078     {
01079       /*  Atoms 1 & 2 are the same, compare atom 3  */
01080       compare_code = strcasecmp(new_node->atom3name, 
01081             tree->atom3name);
01082 
01083       if (compare_code == 0)
01084       {
01085         /*  All three atoms were the same, this */
01086         /*  is a duplicate.  Print a warning    */
01087         /*  message, replace the current values,*/
01088         /*  and free the new node    */
01089         //****** BEGIN CHARMM/XPLOR type changes
01090         /* we do not care about identical replacement */
01091         if ((tree->forceconstant != new_node->forceconstant) ||
01092             (tree->angle != new_node->angle) ||
01093             (tree->k_ub != new_node->k_ub) ||
01094             (tree->r_ub != new_node->r_ub))
01095         {
01096           iout << "\n" << iWARN << "DUPLICATE ANGLE ENTRY FOR "
01097             << new_node->atom1name << "-"
01098             << new_node->atom2name << "-"
01099             << new_node->atom3name
01100             << "\nPREVIOUS VALUES  k="
01101             << tree->forceconstant << "  theta0="
01102             << tree->angle << " k_ub="
01103             << tree->k_ub << " r_ub="
01104             << tree->r_ub
01105             << "\n   USING VALUES  k="
01106             << new_node->forceconstant << "  theta0="
01107             << new_node->angle << " k_ub="
01108             << new_node->k_ub << " r_ub=" << new_node->r_ub 
01109             << "\n" << endi;
01110 
01111           tree->forceconstant=new_node->forceconstant;
01112           tree->angle=new_node->angle;
01113           tree->k_ub=new_node->k_ub;
01114           tree->r_ub=new_node->r_ub;
01115         }
01116         //****** END CHARMM/XPLOR type changes
01117 
01118         delete new_node;
01119 
01120         return(tree);
01121       }
01122     }
01123   }
01124 
01125   /*  Didn't find a duplicate, so if the new_node is smaller  */
01126   /*  than the current head, add it to the left child.  Otherwise */
01127   /*  add it to the right child.          */
01128   if (compare_code < 0)
01129   {
01130     tree->left = add_to_angle_tree(new_node, tree->left);
01131   }
01132   else
01133   {
01134     tree->right = add_to_angle_tree(new_node, tree->right);
01135   }
01136 
01137   return(tree);
01138 }
01139 /*      END OF FUNCTION add_to_angle_tree    */
01140 
01141 /************************************************************************/
01142 /*                  */
01143 /*      FUNCTION add_dihedral_param      */
01144 /*                  */
01145 /*   INPUTS:                */
01146 /*  buf - line from paramter file with dihedral parameters    */
01147 /*                  */
01148 /*  this function adds an dihedral parameter.  It parses up the     */
01149 /*   input line and then adds it to the binary tree used to store the   */
01150 /*   dihedral parameters.            */
01151 /*                  */
01152 /************************************************************************/
01153 
01154 void Parameters::add_dihedral_param(char *buf, FILE *fd)
01155 
01156 {
01157   char atom1name[11];       //  Type of atom 1
01158   char atom2name[11];       //  Type of atom 2
01159   char atom3name[11];       //  Type of atom 3
01160   char atom4name[11];       //  Type of atom 4
01161   Real forceconstant;       //  Force constant
01162   int periodicity;       //  Periodicity
01163   Real phase_shift;       //  Phase shift
01164   int read_count;         //  Count from sscanf
01165   struct dihedral_params *new_node;  //  New node
01166   int multiplicity;       //  Multiplicity for bonds
01167   int i;           //  Loop counter
01168   char buffer[513];       //  Buffer for new line
01169   int ret_code;         //  Return code
01170 
01171   //****** BEGIN CHARMM/XPLOR type changes
01172   /*  Parse up the input line using sscanf      */
01173   if (paramType == paraXplor)
01174   {
01175     /* read XPLOR format */
01176     read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n", 
01177        atom1name, atom2name, atom3name, atom4name, &multiplicity,
01178        &forceconstant, &periodicity, &phase_shift);
01179   }
01180   else if (paramType == paraCharmm)
01181   {
01182     /* read CHARMM format */
01183     read_count=sscanf(buf, "%s %s %s %s %f %d %f\n", 
01184        atom1name, atom2name, atom3name, atom4name,
01185        &forceconstant, &periodicity, &phase_shift);
01186     multiplicity=1; 
01187   }
01188 
01189   if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
01190   {
01191     char err_msg[512];
01192 
01193     sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01194     NAMD_die(err_msg);
01195   }
01196   else if ( (read_count != 7) && (paramType == paraCharmm) )
01197   {
01198     char err_msg[512];
01199 
01200     sprintf(err_msg, "BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
01201     NAMD_die(err_msg);
01202   }
01203 
01204   if ( (read_count == 4) && (paramType == paraXplor) )
01205   //****** END CHARMM/XPLOR type changes
01206   {
01207     read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n", 
01208           &forceconstant, &periodicity, &phase_shift);
01209 
01210     /*  Check to make sure we got what we expected    */
01211     if (read_count != 3)
01212     {
01213       char err_msg[512];
01214 
01215       sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01216       NAMD_die(err_msg);
01217     }
01218 
01219     multiplicity = 1;
01220   }
01221 
01222   if (multiplicity > MAX_MULTIPLICITY)
01223   {
01224     char err_msg[181];
01225 
01226     sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
01227        multiplicity, MAX_MULTIPLICITY);
01228     NAMD_die(err_msg);
01229   }
01230 
01231   /*  Allocate new node            */
01232   new_node = new dihedral_params;
01233 
01234   if (new_node == NULL)
01235   {
01236     NAMD_die("memory allocation failed in Parameters::add_dihedral_param\n");
01237   }
01238 
01239   /*  Assign all of the values for this node.  Notice that since  */
01240   /*  the dihedrals and impropers are implemented with a linked   */
01241   /*  list rather than a binary tree, we don't really care about  */
01242   /*  the order of the atoms any more        */
01243   strcpy(new_node->atom1name, atom1name);
01244   strcpy(new_node->atom2name, atom2name);
01245   strcpy(new_node->atom3name, atom3name);
01246   strcpy(new_node->atom4name, atom4name);
01247   new_node->atom1wild = ! strcasecmp(atom1name, "X");
01248   new_node->atom2wild = ! strcasecmp(atom2name, "X");
01249   new_node->atom3wild = ! strcasecmp(atom3name, "X");
01250   new_node->atom4wild = ! strcasecmp(atom4name, "X");
01251   new_node->multiplicity = multiplicity;
01252   if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01253   new_node->values[0].k = forceconstant;
01254   new_node->values[0].n = periodicity;
01255   new_node->values[0].delta = phase_shift;
01256 
01257   new_node->next = NULL;
01258 
01259   //  If the multiplicity is greater than 1, then read in other parameters
01260   if (multiplicity > 1)
01261   {
01262     for (i=1; i<multiplicity; i++)
01263     {
01264       ret_code = NAMD_read_line(fd, buffer);
01265 
01266       //  Get rid of comments at the end of a line
01267       if (ret_code == 0)
01268       {
01269         NAMD_remove_comment(buffer);
01270       }
01271 
01272       //  Keep reading lines until we get one that isn't blank
01273       while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
01274       {
01275         ret_code = NAMD_read_line(fd, buffer);
01276       }
01277 
01278       if (ret_code != 0)
01279       {
01280         NAMD_die("EOF encoutner in middle of multiple dihedral");
01281       }
01282 
01283       read_count=sscanf(buffer, "%f %d %f\n", 
01284             &forceconstant, &periodicity, &phase_shift);
01285 
01286       if (read_count != 3)
01287       {
01288         char err_msg[512];
01289 
01290         sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
01291         NAMD_die(err_msg);
01292       }
01293 
01294       if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01295       new_node->values[i].k = forceconstant;
01296       new_node->values[i].n = periodicity;
01297       new_node->values[i].delta = phase_shift;
01298     }
01299   }
01300 
01301   //****** BEGIN CHARMM/XPLOR type changes
01302   /*  Add this node to the list          */
01303   if (paramType == paraXplor)
01304   {
01305     add_to_dihedral_list(new_node); // XPLOR
01306   }
01307   else if (paramType == paraCharmm)
01308   {
01309     add_to_charmm_dihedral_list(new_node); // CHARMM
01310   }
01311  //****** END CHARMM/XPLOR type changes
01312 
01313   return;
01314 }
01315 /*      END OF FUNCTION add_dihedral_param    */
01316 
01317 /************************************************************************/
01318 /*                  */
01319 /*      FUNCTION add_to_dihedral_list      */
01320 /*                  */
01321 /*   INPUTS:                */
01322 /*  new_node - node that is to be added to dihedral_list    */
01323 /*                  */
01324 /*  this function adds a new dihedral parameter to the linked list  */
01325 /*   of dihedral parameters.  First, it checks for duplicates.  If a    */
01326 /*   duplicate is found, a warning message is printed, the old values   */
01327 /*   are replaced with the new values, and the new node is freed.  If   */
01328 /*   Otherwise, the node is added to the list.  This list is arranged   */
01329 /*   so that bods with wildcards are placed at the tail of the list.    */
01330 /*   This will guarantee that if we just do a linear search, we will    */
01331 /*   always find an exact match before a wildcard match.    */
01332 /*                  */
01333 /************************************************************************/
01334 
01335 void Parameters::add_to_dihedral_list(
01336         struct dihedral_params *new_node)
01337 
01338 {
01339   static struct dihedral_params *ptr;   //  position within list
01340   static struct dihedral_params *tail;  //  Pointer to the end of 
01341                 //  the list so we can add
01342                 //  entries to the end of the
01343                 //  list in constant time
01344   int i;              //  Loop counter
01345 
01346   /*  If the list is currently empty, then the new node is the list*/
01347   if (dihedralp == NULL)
01348   {
01349     dihedralp=new_node;
01350     tail=new_node;
01351 
01352     return;
01353   }
01354 
01355   /*  The list isn't empty, so check for a duplicate    */
01356   ptr=dihedralp;
01357 
01358   while (ptr != NULL)
01359   {
01360     if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01361            (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01362            (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01363            (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01364          ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01365            (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01366            (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01367            (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
01368     {
01369       /*  Found a duplicate        */
01370       //****** BEGIN CHARMM/XPLOR type changes
01371       /* we do not care about identical replacement */
01372       int echoWarn=0;  // echo warning messages ?
01373 
01374       if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
01375       
01376       if (!echoWarn)
01377       {
01378         for (i=0; i<ptr->multiplicity; i++)
01379         {
01380           if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
01381           if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
01382           if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
01383         }
01384       }
01385 
01386       if (echoWarn)
01387       {
01388         iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
01389           << ptr->atom1name << "-"
01390           << ptr->atom2name << "-"
01391           << ptr->atom3name << "-"
01392           << ptr->atom4name
01393           << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
01394         
01395         for (i=0; i<ptr->multiplicity; i++)
01396         {
01397           iout     << "  k=" << ptr->values[i].k
01398                    << "  n=" << ptr->values[i].n
01399                    << "  delta=" << ptr->values[i].delta;
01400         }
01401 
01402         iout << "\nUSING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
01403 
01404         for (i=0; i<new_node->multiplicity; i++)
01405         {
01406           iout <<     "  k=" << new_node->values[i].k
01407                    << "  n=" << new_node->values[i].n
01408                    << "  delta=" << new_node->values[i].delta;
01409         }
01410 
01411         iout << endi;
01412 
01413         ptr->multiplicity = new_node->multiplicity;
01414 
01415         for (i=0; i<new_node->multiplicity; i++)
01416         {
01417           ptr->values[i].k = new_node->values[i].k;
01418           ptr->values[i].n = new_node->values[i].n;
01419           ptr->values[i].delta = new_node->values[i].delta;
01420         }
01421 
01422       }
01423       //****** END CHARMM/XPLOR type changes
01424 
01425       delete new_node;
01426 
01427       return;
01428     }
01429 
01430     ptr=ptr->next;
01431   }
01432 
01433   /*  Check to see if we have any wildcards.  Since specific  */
01434   /*  entries are to take precedence, we'll put anything without  */
01435   /*  wildcards at the begining of the list and anything with     */
01436   /*  wildcards at the end of the list.  Then, we can just do a   */
01437   /*  linear search for a bond and be guaranteed to have specific */
01438   /*  entries take precendence over over wildcards          */
01439   if ( new_node->atom1wild ||
01440        new_node->atom2wild ||
01441        new_node->atom3wild ||
01442        new_node->atom4wild )
01443   {
01444     /*  add to the end of the list        */
01445     tail->next=new_node;
01446     tail=new_node;
01447 
01448     return;
01449   }
01450   else
01451   {
01452     /*  add to the head of the list        */
01453     new_node->next=dihedralp;
01454     dihedralp=new_node;
01455 
01456     return;
01457   }
01458 
01459 }
01460 /*    END OF FUNCTION add_to_dihedral_list      */
01461 
01462 //****** BEGIN CHARMM/XPLOR type changes
01463 /************************************************************************/
01464 /*                                                                        */
01465 /*                        FUNCTION add_to_charmm_dihedral_list                */
01466 /*                                                                        */
01467 /*   INPUTS:                                                                */
01468 /*        new_node - node that is to be added to dihedral_list                */
01469 /*                                                                        */
01470 /*        this function adds a new dihedral parameter to the linked list  */
01471 /*   of dihedral parameters in CHARMM format.                           */
01472 /*   First, it checks for duplicates.  If a duplicate is found, a       */
01473 /*   warning message is printed. If the periodicity is the same as of   */
01474 /*   a previous dihedral the old values are replaced with the new       */
01475 /*   values, otherwise, the dihedral is added and the multiplicity is   */
01476 /*   increased.                                                         */
01477 /*   Otherwise, the node is added to the list.  This list is arranged   */
01478 /*   so that bonds with wildcards are placed at the tail of the list.   */
01479 /*   This will guarantee that if we just do a linear search, we will    */
01480 /*   always find an exact match before a wildcard match.                */
01481 /*                                                                        */
01482 /************************************************************************/
01483 
01484 void Parameters::add_to_charmm_dihedral_list(
01485                                 struct dihedral_params *new_node)
01486 
01487 {
01488         static struct dihedral_params *ptr;   //  position within list
01489         static struct dihedral_params *tail;  //  Pointer to the end of 
01490                                               //  the list so we can add
01491                                               //  entries to the end of the
01492                                               //  list in constant time
01493         int i;                                      //  Loop counter
01494         int replace;                          //  replace values?
01495 
01496         // keep track of the last dihedral param read to avoid spurious
01497         // error messages.
01498         static struct dihedral_params last_dihedral; 
01499 
01500         /*  If the list is currently empty, then the new node is the list*/
01501         if (dihedralp == NULL)
01502         {
01503                 dihedralp=new_node;
01504                 tail=new_node;
01505                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01506 
01507                 return;
01508         }
01509 
01510         /*  The list isn't empty, so check for a duplicate                */
01511         ptr=dihedralp;
01512 
01513         while (ptr != NULL)
01514         {
01515                 int same_as_last = 0;
01516                 if (  ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01517                        (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01518                        (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01519                        (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01520                      ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01521                        (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01522                        (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01523                        (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) )
01524                        )
01525                 {
01526                         /*  Found a duplicate                                */
01527                         
01528                         // check for same_as_last.  Note: don't believe the echoWarn crap; it controls
01529                         // not just whether we print warning messages, but whether we actually change
01530                         // values or not!  
01531 
01532                         if ( ( !strcmp(ptr->atom1name, last_dihedral.atom1name) && 
01533                                !strcmp(ptr->atom2name, last_dihedral.atom2name) &&
01534                                !strcmp(ptr->atom3name, last_dihedral.atom3name) &&
01535                                !strcmp(ptr->atom4name, last_dihedral.atom4name)))
01536                           same_as_last = 1;
01537 
01538                         //****** BEGIN CHARMM/XPLOR type changes
01539                         /* we do not care about identical replacement */
01540                         int echoWarn=1;  // echo warning messages ?
01541 
01542                         // ptr->multiplicity will always be >= new_node->multiplicity
01543                         for (i=0; i<ptr->multiplicity; i++)
01544                         {
01545                           if ((ptr->values[i].k == new_node->values[0].k) && 
01546                               (ptr->values[i].n == new_node->values[0].n) &&
01547                               (ptr->values[i].delta == new_node->values[0].delta)) 
01548                           {
01549                             // found an identical replacement
01550                             echoWarn=0; 
01551                             break;
01552                           }
01553 
01554                         }
01555                   
01556                         if (echoWarn)
01557                         {
01558                           if (!same_as_last) {
01559                             iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
01560                                  << ptr->atom1name << "-"
01561                                  << ptr->atom2name << "-"
01562                                  << ptr->atom3name << "-"
01563                                  << ptr->atom4name
01564                                  << "\nPREVIOUS VALUES MULTIPLICITY: " << ptr->multiplicity << "\n";
01565                           }
01566                           replace=0;
01567                           
01568                           for (i=0; i<ptr->multiplicity; i++)
01569                           {
01570                             if (!same_as_last) {
01571                               iout << "  k=" << ptr->values[i].k
01572                                    << "  n=" << ptr->values[i].n
01573                                    << "  delta=" << ptr->values[i].delta << "\n";
01574                             }
01575                             if (ptr->values[i].n == new_node->values[0].n)
01576                             {
01577                               iout << iWARN << "IDENTICAL PERIODICITY! REPLACING OLD VALUES BY: \n";
01578                               ptr->values[i].k = new_node->values[0].k;
01579                               ptr->values[i].delta = new_node->values[0].delta;
01580                               iout << "  k=" << ptr->values[i].k
01581                                    << "  n=" << ptr->values[i].n
01582                                    << "  delta=" << ptr->values[i].delta<< "\n";
01583                               replace=1;
01584                               break;
01585                             }
01586                           }
01587 
01588                           if (!replace)
01589                           {
01590                             ptr->multiplicity += 1;
01591 
01592                             if (ptr->multiplicity > MAX_MULTIPLICITY)
01593                             {
01594                               char err_msg[181];
01595 
01596                               sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
01597                                       ptr->multiplicity, MAX_MULTIPLICITY);
01598                               NAMD_die(err_msg);
01599                             }
01600                             if (!same_as_last) 
01601                               iout << "INCREASING MULTIPLICITY TO: " << ptr->multiplicity << "\n";
01602 
01603                             i= ptr->multiplicity - 1; 
01604                             ptr->values[i].k = new_node->values[0].k;
01605                             ptr->values[i].n = new_node->values[0].n;
01606                             ptr->values[i].delta = new_node->values[0].delta;
01607 
01608                             if (!same_as_last) 
01609                               iout << "  k=" << ptr->values[i].k
01610                                    << "  n=" << ptr->values[i].n
01611                                    << "  delta=" << ptr->values[i].delta<< "\n";
01612                           }
01613                         
01614                           iout << endi;
01615                         } 
01616                         //****** END CHARMM/XPLOR type changes
01617 
01618                         memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01619                         delete new_node;
01620 
01621                         return;
01622                 }
01623 
01624                 ptr=ptr->next;
01625         }
01626 
01627         /*  CHARMM and XPLOR wildcards for dihedrals are luckily the same */
01628         /*  Check to see if we have any wildcards.  Since specific        */
01629         /*  entries are to take precedence, we'll put anything without  */
01630         /*  wildcards at the begining of the list and anything with     */
01631         /*  wildcards at the end of the list.  Then, we can just do a   */
01632         /*  linear search for a bond and be guaranteed to have specific */
01633         /*  entries take precendence over over wildcards                */
01634         if ( new_node->atom1wild ||
01635              new_node->atom2wild ||
01636              new_node->atom3wild ||
01637              new_node->atom4wild )
01638         {
01639                 /*  add to the end of the list                                */
01640                 tail->next=new_node;
01641                 tail=new_node;
01642 
01643                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01644                 return;
01645         }
01646         else
01647         {
01648                 /*  add to the head of the list                                */
01649                 new_node->next=dihedralp;
01650                 dihedralp=new_node;
01651 
01652                 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01653