Parameters.C

Go to the documentation of this file.

/*
   The class Parameters is used to hold all of the parameters read
   in from the parameter files.  The class provides a routine to read in
   parameter files (as many parameter files as desired can be read in) and
   a series of routines that allow the parameters that have been read in
   to be queried.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifndef WIN32
#include <strings.h>
#endif
#include "InfoStream.h"
#include <charm++.h>
#include "Parameters.h"
#include "Communicate.h"
#include "ConfigList.h"
//****** BEGIN CHARMM/XPLOR type changes
#include "SimParameters.h"
//****** END CHARMM/XPLOR type changes

#define MIN_DEBUG_LEVEL 3
//#define DEBUGM
#include "Debug.h"


//  struct bond_params is used to form a binary tree of bond parameters.
//  The two atom names are used to determine the order of the nodes in the
//  tree.  atom1name should ALWAYS be lexically before atom2name

struct bond_params
{
  char atom1name[11];
  char atom2name[11];
  Real forceconstant;
  Real distance;
  Index index;
  struct bond_params *left;
  struct bond_params *right;
};

//  struct angle_params is used to form a binary tree of bond parameters.
//  The three atom names are used to determine the order of the nodes in
//  the tree.  atom1name should ALWAYS be lexically before atom3name

struct angle_params
{
  char atom1name[11];
  char atom2name[11];
  char atom3name[11];
  Real forceconstant;
  int normal;
  Real angle;
  Real k_ub;
  Real r_ub;
  Index index;
  struct angle_params *left;
  struct angle_params *right;
};

//  struct dihedral_params is used to form a linked list of the dihedral
//  parameters.  The linked list is arranged in such a way that any
//  bonds with wildcards are at the end of the list so that a linear
//  search can be done but we will still find exact matches before
//  wildcard matches

struct dihedral_params
{
  char atom1name[11];
  char atom2name[11];
  char atom3name[11];
  char atom4name[11];
  char atom1wild;
  char atom2wild;
  char atom3wild;
  char atom4wild;
  int multiplicity;
  FourBodyConsts values[MAX_MULTIPLICITY];
  Index index;
  struct dihedral_params *next;
  dihedral_params() { memset(this, 0, sizeof(dihedral_params)); }
};

//  struct improper_params is used to form a linked list of the improper
//  parameters.  The linked list is arranged in such a way that any
//  bonds with wildcards are at the end of the list so that a linear
//  search can be done but we will still find exact matches before
//  wildcard matches

struct improper_params
{
  char atom1name[11];
  char atom2name[11];
  char atom3name[11];
  char atom4name[11];
  int multiplicity;
  FourBodyConsts values[MAX_MULTIPLICITY];
  Index index;
  struct improper_params *next;
};

struct crossterm_params
{
  crossterm_params(int dim) : dimension(dim) {
    values = new double[dimension*dimension];
  }
  ~crossterm_params() {
    delete [] values;
  }
  char atom1name[11];
  char atom2name[11];
  char atom3name[11];
  char atom4name[11];
  char atom5name[11];
  char atom6name[11];
  char atom7name[11];
  char atom8name[11];
  int dimension;  // usually 24
  double *values;  // dimension * dimension data
  Index index;
  struct crossterm_params *next;
};

//  struct vdw_params is used to form a binary serach tree of the
//  vdw paramters for a single atom.

struct vdw_params
{
  char atomname[11];
  Real sigma;
  Real epsilon;
  Real sigma14;
  Real epsilon14;
  Index index;
  struct vdw_params *left;
  struct vdw_params *right;
};

//  struct vdw_pair_params is used to form a linked list of the
//  vdw parameters for a pair of atoms

struct vdw_pair_params
{
  char atom1name[11];
  char atom2name[11];
  Real A;
  Real B;
  Real A14;
  Real B14;
  struct vdw_pair_params *next;
};


Parameters::Parameters() {
  initialize();
}

void Parameters::initialize() {

  paramType = -1;

  /*  Set all the pointers to NULL        */
  atomTypeNames=NULL;
  bondp=NULL;
  anglep=NULL;
  improperp=NULL;
  dihedralp=NULL;
  crosstermp=NULL;
  vdwp=NULL;
  vdw_pairp=NULL;
  bond_array=NULL;
  angle_array=NULL;
  dihedral_array=NULL;
  improper_array=NULL;
  crossterm_array=NULL;
  vdw_array=NULL;
  vdw_pair_tree=NULL;
  maxDihedralMults=NULL;
  maxImproperMults=NULL;

  /*  Set all the counts to 0          */
  NumBondParams=0;
  NumAngleParams=0;
  NumDihedralParams=0;
  NumImproperParams=0;
  NumCrosstermParams=0;
  NumVdwParams=0;
  NumVdwPairParams=0;
  NumCosAngles=0;
}

/************************************************************************/
/*                  */
/*      FUNCTION Parameters        */
/*                  */
/*  This is the constructor for the class.  It simlpy sets the      */
/*  pointers to the list and trees to NULL and the count of all the     */
/*  parameters to 0.              */
/*  The type (format) of the input parameters (Xplor,Charmm) is set here. */
/*                  */
/************************************************************************/

Parameters::Parameters(SimParameters *simParams, StringList *f)
{
  initialize();

  //****** BEGIN CHARMM/XPLOR type changes
  if (simParams->paraTypeXplorOn)
  {
    paramType = paraXplor;
  }
  else if (simParams->paraTypeCharmmOn)
  {
    paramType = paraCharmm;
  }
  //****** END CHARMM/XPLOR type changes
  //Test for cos-based angles
  if (simParams->cosAngles) {
    cosAngles = true;
  } else {
    cosAngles = false;
  }

  /* Set up AllFilesRead flag to FALSE.  Once all of the files    */
  /* have been read in, then this will be set to true and the     */
  /* arrays of parameters will be set up        */
  AllFilesRead = FALSE;

  if (NULL != f) 
  {
    do
    {
      //****** BEGIN CHARMM/XPLOR type changes
      if (paramType == paraXplor)
      {
        read_parameter_file(f->data);
      }
      else if (paramType == paraCharmm)
      {
        read_charmm_parameter_file(f->data);
      }
      //****** END CHARMM/XPLOR type changes
      f = f->next;
    } while ( f != NULL );

    done_reading_files();
  }

}
/*      END OF FUNCTION Parameters      */

/************************************************************************/
/*                  */
/*      FUNCTION ~Parameters        */
/*                  */
/*  This is the destructor for this class.  It basically just       */
/*  frees all of the memory allocated for the parameters.    */
/*                  */
/************************************************************************/

Parameters::~Parameters()

{
        if (atomTypeNames)
          delete [] atomTypeNames;

  if (bondp != NULL)
    free_bond_tree(bondp);

  if (anglep != NULL)
    free_angle_tree(anglep);

  if (dihedralp != NULL)
    free_dihedral_list(dihedralp);

  if (improperp != NULL)
    free_improper_list(improperp);

  if (crosstermp != NULL)
    free_crossterm_list(crosstermp);

  if (vdwp != NULL)
    free_vdw_tree(vdwp);

  if (vdw_pairp != NULL)
    free_vdw_pair_list();

  if (bond_array != NULL)
    delete [] bond_array;

  if (angle_array != NULL)
    delete [] angle_array;

  if (dihedral_array != NULL)
    delete [] dihedral_array;

  if (improper_array != NULL)
    delete [] improper_array;

  if (crossterm_array != NULL)
    delete [] crossterm_array;

  if (vdw_array != NULL)
    delete [] vdw_array;
  
  if (vdw_pair_tree != NULL)
    free_vdw_pair_tree(vdw_pair_tree);

  if (maxDihedralMults != NULL)
    delete [] maxDihedralMults;

  if (maxImproperMults != NULL)
    delete [] maxImproperMults;

  for( int i = 0; i < error_msgs.size(); ++i ) {
    delete [] error_msgs[i];
  }
  error_msgs.resize(0);
}
/*      END OF FUNCTION ~Parameters      */

/************************************************************************/
/*                  */
/*      FUNCTION read_paramter_file      */
/*                  */
/*   INPUTS:                */
/*  fname - name of the parameter file to read      */
/*                  */
/*  This function reads in a parameter file and adds the parameters */
/*   from this file to the current group of parameters.  The basic      */
/*   logic of the routine is to read in a line from the file, looks at  */
/*   the first word of the line to determine what kind of parameter we  */
/*   have, and then call the appropriate routine to add the parameter   */
/*   to the parameters that we have.          */
/*                  */
/************************************************************************/

void Parameters::read_parameter_file(char *fname)

{
  char buffer[512];  //  Buffer to store each line of the file
  char first_word[512];  //  First word of the current line
  FILE *pfile;    //  File descriptor for the parameter file

  /*  Check to make sure that we haven't previously been told     */
  /*  that all the files were read        */
  if (AllFilesRead)
  {
    NAMD_die("Tried to read another parameter file after being told that all files were read!");
  }

  /*  Try and open the file          */
  if ( (pfile = Fopen(fname, "r")) == NULL)
  {
    char err_msg[256];

    sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
    NAMD_die(err_msg);
  }

  /*  Keep reading in lines until we hit the EOF      */
  while (NAMD_read_line(pfile, buffer) != -1)
  {
    /*  Get the first word of the line      */
    NAMD_find_first_word(buffer, first_word);

    /*  First, screen out things that we ignore, such as    */
    /*  blank lines, lines that start with '!', lines that  */
    /*  start with "REMARK", lines that start with set",    */
    /*  and most of the HBOND parameters which include      */
    /*  AEXP, REXP, HAEX, AAEX, but not the HBOND statement */
    /*  which is parsed.                                    */
    if ((buffer[0] != '!') && 
        !NAMD_blank_string(buffer) &&
        (strncasecmp(first_word, "REMARK", 6) != 0) &&
        (strcasecmp(first_word, "set")!=0) &&
        (strncasecmp(first_word, "AEXP", 4) != 0) &&
        (strncasecmp(first_word, "REXP", 4) != 0) &&
        (strncasecmp(first_word, "HAEX", 4) != 0) &&
        (strncasecmp(first_word, "AAEX", 4) != 0) &&
        (strncasecmp(first_word, "NBOND", 5) != 0) &&
        (strncasecmp(first_word, "CUTNB", 5) != 0) &&
        (strncasecmp(first_word, "END", 3) != 0) &&
        (strncasecmp(first_word, "CTONN", 5) != 0) &&
        (strncasecmp(first_word, "EPS", 3) != 0) &&
        (strncasecmp(first_word, "VSWI", 4) != 0) &&
        (strncasecmp(first_word, "NBXM", 4) != 0) &&
        (strncasecmp(first_word, "INHI", 4) != 0) )
    {
      /*  Now, call the appropriate function based    */
      /*  on the type of parameter we have    */
      if (strncasecmp(first_word, "bond", 4)==0)
      {
        add_bond_param(buffer);
        NumBondParams++;
      }
      else if (strncasecmp(first_word, "angl", 4)==0)
      {
        add_angle_param(buffer);
        NumAngleParams++;
      }
      else if (strncasecmp(first_word, "dihe", 4)==0)
      {
        add_dihedral_param(buffer, pfile);
        NumDihedralParams++;
      }
      else if (strncasecmp(first_word, "impr", 4)==0)
      {
        add_improper_param(buffer, pfile);
        NumImproperParams++;
      }
      else if (strncasecmp(first_word, "nonb", 4)==0)
      {
        add_vdw_param(buffer);
        NumVdwParams++; 
      }
      else if (strncasecmp(first_word, "nbfi", 4)==0)
      {
        add_vdw_pair_param(buffer);
        NumVdwPairParams++; 
      }
      else if (strncasecmp(first_word, "hbon", 4)==0)
      {
        add_hb_pair_param(buffer);
      }
      else
      {
        /*  This is an unknown paramter.        */
        /*  This is BAD        */
        char err_msg[512];

        sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
           fname, buffer);
        NAMD_die(err_msg);
      }
    }
  }

  /*  Close the file            */
  Fclose(pfile);

  return;
}
/*      END OF FUNCTION read_paramter_file    */

//****** BEGIN CHARMM/XPLOR type changes
/************************************************************************/
/*                                                                        */
/*                        FUNCTION read_charmm_paramter_file                */
/*                                                                        */
/*   INPUTS:                                                                */
/*        fname - name of the parameter file to read                        */
/*                                                                        */
/*        This function reads in a CAHRMM parameter file and adds the     */ 
/*   parameters from this file to the current group of parameters.      */
/*   The basic logic of the routine is to first find out what type of   */
/*   parameter we have in the file. Then look at each line in turn      */
/*   and call the appropriate routine to add the parameters until we hit*/
/*   a new type of parameter or EOF.                                    */
/*                                                                        */
/************************************************************************/

void Parameters::read_charmm_parameter_file(char *fname)

{
  int  par_type=0;         //  What type of parameter are we currently
                           //  dealing with? (vide infra)
  int  skipline;           //  skip this line?
  int  skipall = 0;        //  skip rest of file;
  char buffer[512];           //  Buffer to store each line of the file
  char first_word[512];           //  First word of the current line
  FILE *pfile;                   //  File descriptor for the parameter file

  /*  Check to make sure that we haven't previously been told     */
  /*  that all the files were read                                */
  if (AllFilesRead)
  {
    NAMD_die("Tried to read another parameter file after being told that all files were read!");
  }

  /*  Try and open the file                                        */
  if ( (pfile = fopen(fname, "r")) == NULL)
  {
    char err_msg[256];

    sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
    NAMD_die(err_msg);
  }

  /*  Keep reading in lines until we hit the EOF                        */
  while (NAMD_read_line(pfile, buffer) != -1)
  {
    /*  Get the first word of the line                        */
    NAMD_find_first_word(buffer, first_word);
    skipline=0;

    /*  First, screen out things that we ignore.                   */   
    /*  blank lines, lines that start with '!' or '*', lines that  */
    /*  start with "END".                                          */
    if (!NAMD_blank_string(buffer) &&
        (strncmp(first_word, "!", 1) != 0) &&
         (strncmp(first_word, "*", 1) != 0) &&
        (strncasecmp(first_word, "END", 3) != 0))
    {
      if ( skipall ) {
        iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
        break;
      }
      /*  Now, determine the apropriate parameter type.   */
      if (strncasecmp(first_word, "bond", 4)==0)
      {
        par_type=1; skipline=1;
      }
      else if (strncasecmp(first_word, "angl", 4)==0)
      {
        par_type=2; skipline=1;
      }
      else if (strncasecmp(first_word, "thet", 4)==0)
      {
        par_type=2; skipline=1;
      }
      else if (strncasecmp(first_word, "dihe", 4)==0)
      {
        par_type=3; skipline=1;
      }
      else if (strncasecmp(first_word, "phi", 3)==0)
      {
        par_type=3; skipline=1;
      }
      else if (strncasecmp(first_word, "impr", 4)==0)
      {
        par_type=4; skipline=1;
      }
      else if (strncasecmp(first_word, "imph", 4)==0)
      {
        par_type=4; skipline=1;
      }
      else if (strncasecmp(first_word, "nbon", 4)==0)
      {
        par_type=5; skipline=1;
      }
      else if (strncasecmp(first_word, "nonb", 4)==0)
      {
        par_type=5; skipline=1;
      }
      else if (strncasecmp(first_word, "nbfi", 4)==0)
      {
        par_type=6; skipline=1;
      }
      else if (strncasecmp(first_word, "hbon", 4)==0)
      {
        par_type=7; skipline=1;
      }
      else if (strncasecmp(first_word, "cmap", 4)==0)
      {
        par_type=8; skipline=1;
      }
      else if (strncasecmp(first_word, "read", 4)==0)
      {
        skip_stream_read(buffer, pfile);  skipline=1;
      }
      else if (strncasecmp(first_word, "return", 4)==0)
      {
        skipall=1;  skipline=1;
      }
      else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
               (strncasecmp(first_word, "grou", 4) == 0) ||
               (strncasecmp(first_word, "cdie", 4) == 0) ||
               (strncasecmp(first_word, "shif", 4) == 0) ||
               (strncasecmp(first_word, "vgro", 4) == 0) ||
               (strncasecmp(first_word, "vdis", 4) == 0) ||
               (strncasecmp(first_word, "vswi", 4) == 0) ||
               (strncasecmp(first_word, "cutn", 4) == 0) ||
               (strncasecmp(first_word, "ctof", 4) == 0) ||
               (strncasecmp(first_word, "cton", 4) == 0) ||
               (strncasecmp(first_word, "eps", 3) == 0) ||
               (strncasecmp(first_word, "e14f", 4) == 0) ||
               (strncasecmp(first_word, "wmin", 4) == 0) ||
               (strncasecmp(first_word, "aexp", 4) == 0) ||
               (strncasecmp(first_word, "rexp", 4) == 0) ||
               (strncasecmp(first_word, "haex", 4) == 0) ||
               (strncasecmp(first_word, "aaex", 4) == 0) ||
               (strncasecmp(first_word, "noac", 4) == 0) ||
               (strncasecmp(first_word, "hbno", 4) == 0) ||
               (strncasecmp(first_word, "cuth", 4) == 0) ||
               (strncasecmp(first_word, "ctof", 4) == 0) ||
               (strncasecmp(first_word, "cton", 4) == 0) ||
               (strncasecmp(first_word, "cuth", 4) == 0) ||
               (strncasecmp(first_word, "ctof", 4) == 0) ||
               (strncasecmp(first_word, "cton", 4) == 0) ) 
      {
        if ((par_type != 5) && (par_type != 6) && (par_type != 7))
        {
          char err_msg[512];

          sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
          NAMD_die(err_msg);
        }
        else 
        {
          skipline = 1;
        }
      }        
      else if (par_type == 0)
      {
        /*  This is an unknown paramter.        */
        /*  This is BAD                                */
        char err_msg[512];

        sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
        NAMD_die(err_msg);
      }
    }
    else
    {
      skipline=1;
    }

    if ( (par_type != 0) && (!skipline) )
    {
      /*  Now, call the appropriate function based    */
      /*  on the type of parameter we have                */
      /*  I know, this should really be a switch ...  */
      if (par_type == 1)
      {
        add_bond_param(buffer);
        NumBondParams++;
      }
      else if (par_type == 2)
      {
        add_angle_param(buffer);
        NumAngleParams++;
      }
      else if (par_type == 3)
      {
        add_dihedral_param(buffer, pfile);
        NumDihedralParams++;
      }
      else if (par_type == 4)
      {
        add_improper_param(buffer, pfile);
        NumImproperParams++;
      }
      else if (par_type == 5)
      {
        add_vdw_param(buffer);
        NumVdwParams++;
      }
      else if (par_type == 6)
      {
        add_vdw_pair_param(buffer);
        NumVdwPairParams++; 
      }
      else if (par_type == 7)
      {
        add_hb_pair_param(buffer);                  
      }
      else if (par_type == 8)
      {
        add_crossterm_param(buffer, pfile);                  
        NumCrosstermParams++;
      }
      else
      {
        /*  This really should not occour!      */
        /*  This is an internal error.          */
        /*  This is VERY BAD                        */
        char err_msg[512];

        sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
        NAMD_die(err_msg);
      }
    }
  }

  /*  Close the file                                                */
  fclose(pfile);

  return;
}
/*                        END OF FUNCTION read_charmm_paramter_file                */
//****** END CHARMM/XPLOR type changes


void Parameters::skip_stream_read(char *buf, FILE *fd) {

  char buffer[513];
  char first_word[513];
  char s1[128];
  char s2[128];
  int rval = sscanf(buf, "%s %s", s1, s2);
  if (rval != 2) {
        char err_msg[512];
        sprintf(err_msg, "BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
        NAMD_die(err_msg);
  }
  if ( ! strncasecmp(s2,"PARA",4) ) return;  // read parameters
  
  iout << iINFO << "SKIPPING " << s2 << " SECTION IN STREAM FILE\n" << endi;

  while (NAMD_read_line(fd, buffer) != -1)
  {
    // read until we find "END"
    NAMD_find_first_word(buffer, first_word);
    if (!NAMD_blank_string(buffer) &&
        (strncmp(first_word, "!", 1) != 0) &&
         (strncmp(first_word, "*", 1) != 0) &&
         (strncasecmp(first_word, "END", 3) == 0) ) return;
  }

}


/************************************************************************/
/*                  */
/*      FUNCTION add_bond_param        */
/*                  */
/*   INPUTS:                */
/*  buf - Line from parameter file containing bond parameters  */
/*                  */
/*  This function adds a new bond paramter to the binary tree of    */
/*   angle paramters that we have.  If a duplicate is found, a warning  */
/*   message is printed and the new parameters are used.    */
/*                  */
/************************************************************************/

void Parameters::add_bond_param(char *buf)

{
  char atom1name[11];    //  Atom type for atom 1
  char atom2name[11];    //  Atom type for atom 2
  Real forceconstant;    //  Force constant for bond
  Real distance;      //  Rest distance for bond
  int read_count;      //  Count from sscanf
  struct bond_params *new_node;  //  New node in tree

  //****** BEGIN CHARMM/XPLOR type changes
  /*  Use sscanf to parse up the input line      */
  if (paramType == paraXplor)
  {
    /* read XPLOR format */
    read_count=sscanf(buf, "%*s %s %s %f %f\n", atom1name, atom2name, 
       &forceconstant, &distance);
  }
  else if (paramType == paraCharmm)
  {
    /* read CHARMM format */
    read_count=sscanf(buf, "%s %s %f %f\n", atom1name, atom2name, 
       &forceconstant, &distance);
  }
  //****** END CHARMM/XPLOR type changes

  /*  Check to make sure we found everything we expeceted    */
  if (read_count != 4)
  {
    char err_msg[512];

    if (paramType == paraXplor)
    {
      sprintf(err_msg, "BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
    }
    else if (paramType == paraCharmm)
    {
      sprintf(err_msg, "BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
    }
    NAMD_die(err_msg);
  }

  /*  Allocate a new node            */
  new_node = new bond_params;

  if (new_node == NULL)
  {
    NAMD_die("memory allocation failed in Parameters::add_bond_param\n");
  }

  /*  Order the atoms so that the atom that comes alphabetically  */
  /*  first is atom 1.  Since the bond is symmetric, it doesn't   */
  /*  matter physically which atom is first.  And this allows the */
  /*  search of the binary tree to be done in a logical manner    */
  if (strcasecmp(atom1name, atom2name) < 0)
  {
    strcpy(new_node->atom1name, atom1name);
    strcpy(new_node->atom2name, atom2name);
  }
  else
  {
    strcpy(new_node->atom2name, atom1name);
    strcpy(new_node->atom1name, atom2name);
  }

  /*  Assign force constant and distance        */
  new_node->forceconstant = forceconstant;
  new_node->distance = distance;

  /*  Set pointers to null          */
  new_node->left = NULL;
  new_node->right = NULL;

  /*  Make call to recursive call to actually add the node to the */
  /*  tree              */
  bondp=add_to_bond_tree(new_node, bondp);

  return;
}
/*      END OF FUNCTION add_bond_param      */

/************************************************************************/
/*                  */
/*      FUNCTION add_to_bond_tree      */
/*                  */
/*   INPUTS:                */
/*  new_node - Node to add to the tree        */
/*  tree - tree to add the node to          */
/*                  */
/*   OUTPUTS:                */
/*  ths function returns a pointer to the new tree with the node    */
/*   added to it.  Most of the time it will be the same pointer as was  */
/*   passed in, but not if the current tree is empty.      */
/*                  */
/*  this is a receursive function that adds a node to the binary    */
/*   tree used to store bond parameters.        */
/*                  */
/************************************************************************/

struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
             struct bond_params *tree)

{
  int compare_code;  //  Results from strcasecmp

  /*  If the tree is currently empty, then the new tree consists  */
  /*  only of the new node          */
  if (tree == NULL)
    return(new_node);

  /*  Compare the atom1 name from the new node and the head of    */
  /*  the tree              */
  compare_code = strcasecmp(new_node->atom1name, tree->atom1name);

  /*  Check to see if they are the same        */
  if (compare_code == 0)
  {
    /*  The atom 1 names are the same, compare atom 2  */
    compare_code = strcasecmp(new_node->atom2name, tree->atom2name);

    /*  If atom 1 AND atom 2 are the same, we have a duplicate */
    if (compare_code == 0)
    {
      /*  We have a duplicate.  So print out a warning*/
      /*  message.  Then assign the new values to the */
      /*  tree and free the new_node      */
      //****** BEGIN CHARMM/XPLOR type changes
      /* we do not care about identical replacement */
      if ((tree->forceconstant != new_node->forceconstant) || 
          (tree->distance != new_node->distance))
      {
        iout << "\n" << iWARN << "DUPLICATE BOND ENTRY FOR "
          << new_node->atom1name << "-"
          << new_node->atom2name
          << "\nPREVIOUS VALUES  k=" << tree->forceconstant
          << "  x0=" << tree->distance
          << "\n   USING VALUES  k=" << new_node->forceconstant
          << "  x0=" << new_node->distance
          << "\n" << endi;

        tree->forceconstant=new_node->forceconstant;
        tree->distance=new_node->distance;
      }
      //****** END CHARMM/XPLOR type changes

      delete new_node;

      return(tree);
    }
  }

  /*  We don't have a duplicate, so if the new value is less      */
  /*  than the head of the tree, add it to the left child,   */
  /*  otherwise add it to the right child        */
  if (compare_code < 0)
  {
    tree->left = add_to_bond_tree(new_node, tree->left);
  }
  else
  {
    tree->right = add_to_bond_tree(new_node, tree->right);
  }

  return(tree);
}
/*    END OF FUNCTION add_to_bond_tree      */

/************************************************************************/
/*                  */
/*      FUNCTION add_angle_param      */
/*                  */
/*   INPUTS:                */
/*  buf - line from paramter file with angle parameters    */
/*                  */
/*  this function adds an angle parameter.  It parses up the input  */
/*   line and then adds it to the binary tree used to store the angle   */
/*   parameters.              */
/*                  */
/************************************************************************/

void Parameters::add_angle_param(char *buf)

{
  char atom1name[11];    // Type for atom 1
  char atom2name[11];    // Type for atom 2
  char atom3name[11];    // Type for atom 3
  char norm[4]="xxx";
  Real forceconstant;    // Force constant
  Real angle;      // Theta 0
  Real k_ub;      // Urey-Bradley force constant
  Real r_ub;      // Urey-Bradley distance
  int read_count;      // count from sscanf
  struct angle_params *new_node;  // new node in tree

  //****** BEGIN CHARMM/XPLOR type changes
  /*  parse up the input line with sscanf        */
  if (paramType == paraXplor)
  {
    /* read XPLOR format */
    read_count=sscanf(buf, "%*s %s %s %s %f %f UB %f %f\n", 
       atom1name, atom2name, atom3name, &forceconstant, &angle,
       &k_ub, &r_ub);
  }
  else if ((paramType == paraCharmm) && cosAngles) {
    read_count=sscanf(buf, "%s %s %s %f %f %3s %f %f\n", 
       atom1name, atom2name, atom3name, &forceconstant, &angle, norm,
       &k_ub, &r_ub);
//    printf("%s\n", buf);
//    printf("Data: %s %s %s %f %f %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle, norm, k_ub, r_ub);
  }  
  else if (paramType == paraCharmm)
  {
    /* read CHARMM format */
    read_count=sscanf(buf, "%s %s %s %f %f %f %f\n", 
       atom1name, atom2name, atom3name, &forceconstant, &angle,
       &k_ub, &r_ub);
//    printf("%s\n", buf);
//    printf("Data: %s %s %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle);
  }
  //****** END CHARMM/XPLOR type changes

  /*  Check to make sure we got what we expected      */
  if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
  {
    char err_msg[512];

    if (paramType == paraXplor)
    {
      sprintf(err_msg, "BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
    }
    else if (paramType == paraCharmm)
    {
      sprintf(err_msg, "BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
    }
    NAMD_die(err_msg);
  }

  /*  Allocate the new node          */
  new_node = new angle_params;

  if (new_node == NULL)
  {
    NAMD_die("memory allocation failed in Parameters::add_angle_param");
  }

  /*  As with the bond, we want the atom type is comes first  */
  /*  alphbetically first between atom 1 and atom 3 to be in      */
  /*  atom 1 so that we can search the tree reliably.    */
  if (strcasecmp(atom1name, atom3name) < 0)
  {
    strcpy(new_node->atom1name, atom1name);
    strcpy(new_node->atom2name, atom2name);
    strcpy(new_node->atom3name, atom3name);
  }
  else
  {
    strcpy(new_node->atom3name, atom1name);
    strcpy(new_node->atom2name, atom2name);
    strcpy(new_node->atom1name, atom3name);
  }

  /*  Assign the constants and pointer values      */
  new_node->forceconstant = forceconstant;
  new_node->angle = angle;

  if (cosAngles) {
    if (strcmp("cos",norm)==0) {
//      iout << "Info: Using cos mode for angle " << buf << endl;
      NumCosAngles++;
      new_node->normal = 0;
    } else {
//      iout << "Info: Using x^2 mode for angle " << buf << endl;
      new_node->normal = 1;
    }
  } else {
    new_node->normal = 1;
  }

  if (read_count == 7)
  {
    //  Urey-Bradley constants
    if (new_node->normal == 0) {
      NAMD_die("ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
    }
    new_node->k_ub = k_ub;
    new_node->r_ub = r_ub;
  }
  else
  {
    new_node->k_ub = 0.0;
    new_node->r_ub = 0.0;
  }

  new_node->left = NULL;
  new_node->right = NULL;

  /*  Insert it into the tree          */
  anglep = add_to_angle_tree(new_node, anglep);

  return;
}
/*      END OF FUNCTION add_angle_param      */

/************************************************************************/
/*                  */
/*      FUNCTION add_to_angle_tree      */
/*                  */
/*   INPUTS:                */
/*  new_node - new node to add to the angle tree      */
/*  tree - tree to add the node to          */
/*                  */
/*   OUTPUTS:                */
/*  the function returns a pointer to the new tree with the node    */
/*   added.  Most of the time, this will be the same as the value passed*/
/*   in, but not in the case where the tree is empty.      */
/*                  */
/*  this is a recursive function that adds an angle parameter  */
/*   to the binary tree storing the angle parameters.  If a duplicate   */
/*   is found, a warning message is printed, the current values in the  */
/*   tree are replaced with the new values, and the new node is free'd  */
/*                  */
/************************************************************************/

struct angle_params *Parameters::add_to_angle_tree(struct angle_params *new_node,
             struct angle_params *tree)

{
  int compare_code;  //  Return code from strcasecmp

  /*  If the tree is empty, then the new_node is the tree    */
  if (tree == NULL)
    return(new_node);

  /*  Compare atom 1 from the new node and the head of the tree   */
  compare_code = strcasecmp(new_node->atom1name, tree->atom1name);

  if (compare_code == 0)
  {
    /*  Atom 1 is the same, compare atom 2      */
    compare_code = strcasecmp(new_node->atom2name, tree->atom2name);

    if (compare_code == 0)
    {
      /*  Atoms 1 & 2 are the same, compare atom 3  */
      compare_code = strcasecmp(new_node->atom3name, 
            tree->atom3name);

      if (compare_code == 0)
      {
        /*  All three atoms were the same, this */
        /*  is a duplicate.  Print a warning    */
        /*  message, replace the current values,*/
        /*  and free the new node    */
        //****** BEGIN CHARMM/XPLOR type changes
        /* we do not care about identical replacement */
        if ((tree->forceconstant != new_node->forceconstant) ||
            (tree->angle != new_node->angle) ||
            (tree->k_ub != new_node->k_ub) ||
            (tree->r_ub != new_node->r_ub))
        {
          iout << "\n" << iWARN << "DUPLICATE ANGLE ENTRY FOR "
            << new_node->atom1name << "-"
            << new_node->atom2name << "-"
            << new_node->atom3name
            << "\nPREVIOUS VALUES  k="
            << tree->forceconstant << "  theta0="
            << tree->angle << " k_ub="
            << tree->k_ub << " r_ub="
            << tree->r_ub
            << "\n   USING VALUES  k="
            << new_node->forceconstant << "  theta0="
            << new_node->angle << " k_ub="
            << new_node->k_ub << " r_ub=" << new_node->r_ub 
            << "\n" << endi;

          tree->forceconstant=new_node->forceconstant;
          tree->angle=new_node->angle;
          tree->k_ub=new_node->k_ub;
          tree->r_ub=new_node->r_ub;
        }
        //****** END CHARMM/XPLOR type changes

        delete new_node;

        return(tree);
      }
    }
  }

  /*  Didn't find a duplicate, so if the new_node is smaller  */
  /*  than the current head, add it to the left child.  Otherwise */
  /*  add it to the right child.          */
  if (compare_code < 0)
  {
    tree->left = add_to_angle_tree(new_node, tree->left);
  }
  else
  {
    tree->right = add_to_angle_tree(new_node, tree->right);
  }

  return(tree);
}
/*      END OF FUNCTION add_to_angle_tree    */

/************************************************************************/
/*                  */
/*      FUNCTION add_dihedral_param      */
/*                  */
/*   INPUTS:                */
/*  buf - line from paramter file with dihedral parameters    */
/*                  */
/*  this function adds an dihedral parameter.  It parses up the     */
/*   input line and then adds it to the binary tree used to store the   */
/*   dihedral parameters.            */
/*                  */
/************************************************************************/

void Parameters::add_dihedral_param(char *buf, FILE *fd)

{
  char atom1name[11];       //  Type of atom 1
  char atom2name[11];       //  Type of atom 2
  char atom3name[11];       //  Type of atom 3
  char atom4name[11];       //  Type of atom 4
  Real forceconstant;       //  Force constant
  int periodicity;       //  Periodicity
  Real phase_shift;       //  Phase shift
  int read_count;         //  Count from sscanf
  struct dihedral_params *new_node;  //  New node
  int multiplicity;       //  Multiplicity for bonds
  int i;           //  Loop counter
  char buffer[513];       //  Buffer for new line
  int ret_code;         //  Return code

  //****** BEGIN CHARMM/XPLOR type changes
  /*  Parse up the input line using sscanf      */
  if (paramType == paraXplor)
  {
    /* read XPLOR format */
    read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n", 
       atom1name, atom2name, atom3name, atom4name, &multiplicity,
       &forceconstant, &periodicity, &phase_shift);
  }
  else if (paramType == paraCharmm)
  {
    /* read CHARMM format */
    read_count=sscanf(buf, "%s %s %s %s %f %d %f\n", 
       atom1name, atom2name, atom3name, atom4name,
       &forceconstant, &periodicity, &phase_shift);
    multiplicity=1; 
  }

  if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
  {
    char err_msg[512];

    sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
    NAMD_die(err_msg);
  }
  else if ( (read_count != 7) && (paramType == paraCharmm) )
  {
    char err_msg[512];

    sprintf(err_msg, "BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
    NAMD_die(err_msg);
  }

  if ( (read_count == 4) && (paramType == paraXplor) )
  //****** END CHARMM/XPLOR type changes
  {
    read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n", 
          &forceconstant, &periodicity, &phase_shift);

    /*  Check to make sure we got what we expected    */
    if (read_count != 3)
    {
      char err_msg[512];

      sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
      NAMD_die(err_msg);
    }

    multiplicity = 1;
  }

  if (multiplicity > MAX_MULTIPLICITY)
  {
    char err_msg[181];

    sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
       multiplicity, MAX_MULTIPLICITY);
    NAMD_die(err_msg);
  }

  /*  Allocate new node            */
  new_node = new dihedral_params;

  if (new_node == NULL)
  {
    NAMD_die("memory allocation failed in Parameters::add_dihedral_param\n");
  }

  /*  Assign all of the values for this node.  Notice that since  */
  /*  the dihedrals and impropers are implemented with a linked   */
  /*  list rather than a binary tree, we don't really care about  */
  /*  the order of the atoms any more        */
  strcpy(new_node->atom1name, atom1name);
  strcpy(new_node->atom2name, atom2name);
  strcpy(new_node->atom3name, atom3name);
  strcpy(new_node->atom4name, atom4name);
  new_node->atom1wild = ! strcasecmp(atom1name, "X");
  new_node->atom2wild = ! strcasecmp(atom2name, "X");
  new_node->atom3wild = ! strcasecmp(atom3name, "X");
  new_node->atom4wild = ! strcasecmp(atom4name, "X");
  new_node->multiplicity = multiplicity;
  if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
  new_node->values[0].k = forceconstant;
  new_node->values[0].n = periodicity;
  new_node->values[0].delta = phase_shift;

  new_node->next = NULL;

  //  If the multiplicity is greater than 1, then read in other parameters
  if (multiplicity > 1)
  {
    for (i=1; i<multiplicity; i++)
    {
      ret_code = NAMD_read_line(fd, buffer);

      //  Get rid of comments at the end of a line
      if (ret_code == 0)
      {
        NAMD_remove_comment(buffer);
      }

      //  Keep reading lines until we get one that isn't blank
      while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
      {
        ret_code = NAMD_read_line(fd, buffer);
      }

      if (ret_code != 0)
      {
        NAMD_die("EOF encoutner in middle of multiple dihedral");
      }

      read_count=sscanf(buffer, "%f %d %f\n", 
            &forceconstant, &periodicity, &phase_shift);

      if (read_count != 3)
      {
        char err_msg[512];

        sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
        NAMD_die(err_msg);
      }

      if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
      new_node->values[i].k = forceconstant;
      new_node->values[i].n = periodicity;
      new_node->values[i].delta = phase_shift;
    }
  }

  //****** BEGIN CHARMM/XPLOR type changes
  /*  Add this node to the list          */
  if (paramType == paraXplor)
  {
    add_to_dihedral_list(new_node); // XPLOR
  }
  else if (paramType == paraCharmm)
  {
    add_to_charmm_dihedral_list(new_node); // CHARMM
  }
 //****** END CHARMM/XPLOR type changes

  return;
}
/*      END OF FUNCTION add_dihedral_param    */

/************************************************************************/
/*                  */
/*      FUNCTION add_to_dihedral_list      */
/*                  */
/*   INPUTS:                */
/*  new_node - node that is to be added to dihedral_list    */
/*                  */
/*  this function adds a new dihedral parameter to the linked list  */
/*   of dihedral parameters.  First, it checks for duplicates.  If a    */
/*   duplicate is found, a warning message is printed, the old values   */
/*   are replaced with the new values, and the new node is freed.  If   */
/*   Otherwise, the node is added to the list.  This list is arranged   */
/*   so that bods with wildcards are placed at the tail of the list.    */
/*   This will guarantee that if we just do a linear search, we will    */
/*   always find an exact match before a wildcard match.    */
/*                  */
/************************************************************************/

void Parameters::add_to_dihedral_list(
        struct dihedral_params *new_node)

{
  static struct dihedral_params *ptr;   //  position within list
  static struct dihedral_params *tail;  //  Pointer to the end of 
                //  the list so we can add
                //  entries to the end of the
                //  list in constant time
  int i;              //  Loop counter

  /*  If the list is currently empty, then the new node is the list*/
  if (dihedralp == NULL)
  {
    dihedralp=new_node;
    tail=new_node;

    return;
  }

  /*  The list isn't empty, so check for a duplicate    */
  ptr=dihedralp;

  while (ptr != NULL)
  {
    if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
           (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
           (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
           (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
         ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
           (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
           (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
           (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
    {
      /*  Found a duplicate        */
      //****** BEGIN CHARMM/XPLOR type changes
      /* we do not care about identical replacement */
      int echoWarn=0;  // echo warning messages ?

      if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
      
      if (!echoWarn)
      {
        for (i=0; i<ptr->multiplicity; i++)
        {
          if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
          if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
          if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
        }
      }

      if (echoWarn)
      {
        iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
          << ptr->atom1name << "-"
          << ptr->atom2name << "-"
          << ptr->atom3name << "-"
          << ptr->atom4name
          << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
        
        for (i=0; i<ptr->multiplicity; i++)
        {
          iout     << "  k=" << ptr->values[i].k
                   << "  n=" << ptr->values[i].n
                   << "  delta=" << ptr->values[i].delta;
        }

        iout << "\nUSING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";

        for (i=0; i<new_node->multiplicity; i++)
        {
          iout <<     "  k=" << new_node->values[i].k
                   << "  n=" << new_node->values[i].n
                   << "  delta=" << new_node->values[i].delta;
        }

        iout << endi;

        ptr->multiplicity = new_node->multiplicity;

        for (i=0; i<new_node->multiplicity; i++)
        {
          ptr->values[i].k = new_node->values[i].k;
          ptr->values[i].n = new_node->values[i].n;
          ptr->values[i].delta = new_node->values[i].delta;
        }

      }
      //****** END CHARMM/XPLOR type changes

      delete new_node;

      return;
    }

    ptr=ptr->next;
  }

  /*  Check to see if we have any wildcards.  Since specific  */
  /*  entries are to take precedence, we'll put anything without  */
  /*  wildcards at the begining of the list and anything with     */
  /*  wildcards at the end of the list.  Then, we can just do a   */
  /*  linear search for a bond and be guaranteed to have specific */
  /*  entries take precendence over over wildcards          */
  if ( new_node->atom1wild ||
       new_node->atom2wild ||
       new_node->atom3wild ||
       new_node->atom4wild )
  {
    /*  add to the end of the list        */
    tail->next=new_node;
    tail=new_node;

    return;
  }
  else
  {
    /*  add to the head of the list        */
    new_node->next=dihedralp;
    dihedralp=new_node;

    return;
  }

}
/*    END OF FUNCTION add_to_dihedral_list      */

//****** BEGIN CHARMM/XPLOR type changes
/************************************************************************/
/*                                                                        */
/*                        FUNCTION add_to_charmm_dihedral_list                */
/*                                                                        */
/*   INPUTS:                                                                */
/*        new_node - node that is to be added to dihedral_list                */
/*                                                                        */
/*        this function adds a new dihedral parameter to the linked list  */
/*   of dihedral parameters in CHARMM format.                           */
/*   First, it checks for duplicates.  If a duplicate is found, a       */
/*   warning message is printed. If the periodicity is the same as of   */
/*   a previous dihedral the old values are replaced with the new       */
/*   values, otherwise, the dihedral is added and the multiplicity is   */
/*   increased.                                                         */
/*   Otherwise, the node is added to the list.  This list is arranged   */
/*   so that bonds with wildcards are placed at the tail of the list.   */
/*   This will guarantee that if we just do a linear search, we will    */
/*   always find an exact match before a wildcard match.                */
/*                                                                        */
/************************************************************************/

void Parameters::add_to_charmm_dihedral_list(
                                struct dihedral_params *new_node)

{
        static struct dihedral_params *ptr;   //  position within list
        static struct dihedral_params *tail;  //  Pointer to the end of 
                                              //  the list so we can add
                                              //  entries to the end of the
                                              //  list in constant time
        int i;                                      //  Loop counter
        int replace;                          //  replace values?

        // keep track of the last dihedral param read to avoid spurious
        // error messages.
        static struct dihedral_params last_dihedral; 

        /*  If the list is currently empty, then the new node is the list*/
        if (dihedralp == NULL)
        {
                dihedralp=new_node;
                tail=new_node;
                memcpy(&last_dihedral, new_node, sizeof(dihedral_params));

                return;
        }

        /*  The list isn't empty, so check for a duplicate                */
        ptr=dihedralp;

        while (ptr != NULL)
        {
                int same_as_last = 0;
                if (  ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
                       (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
                       (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
                       (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
                     ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
                       (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
                       (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
                       (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) )
                       )
                {
                        /*  Found a duplicate                                */
                        
                        // check for same_as_last.  Note: don't believe the echoWarn crap; it controls
                        // not just whether we print warning messages, but whether we actually change
                        // values or not!  

                        if ( ( !strcmp(ptr->atom1name, last_dihedral.atom1name) && 
                               !strcmp(ptr->atom2name, last_dihedral.atom2name) &&
                               !strcmp(ptr->atom3name, last_dihedral.atom3name) &&
                               !strcmp(ptr->atom4name, last_dihedral.atom4name)))
                          same_as_last = 1;

                        //****** BEGIN CHARMM/XPLOR type changes
                        /* we do not care about identical replacement */
                        int echoWarn=1;  // echo warning messages ?

                        // ptr->multiplicity will always be >= new_node->multiplicity
                        for (i=0; i<ptr->multiplicity; i++)
                        {
                          if ((ptr->values[i].k == new_node->values[0].k) && 
                              (ptr->values[i].n == new_node->values[0].n) &&
                              (ptr->values[i].delta == new_node->values[0].delta)) 
                          {
                            // found an identical replacement
                            echoWarn=0; 
                            break;
                          }

                        }
                  
                        if (echoWarn)
                        {
                          if (!same_as_last) {
                            iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
                                 << ptr->atom1name << "-"
                                 << ptr->atom2name << "-"
                                 << ptr->atom3name << "-"
                                 << ptr->atom4name
                                 << "\nPREVIOUS VALUES MULTIPLICITY: " << ptr->multiplicity << "\n";
                          }
                          replace=0;
                          
                          for (i=0; i<ptr->multiplicity; i++)
                          {
                            if (!same_as_last) {
                              iout << "  k=" << ptr->values[i].k
                                   << "  n=" << ptr->values[i].n
                                   << "  delta=" << ptr->values[i].delta << "\n";
                            }
                            if (ptr->values[i].n == new_node->values[0].n)
                            {
                              iout << iWARN << "IDENTICAL PERIODICITY! REPLACING OLD VALUES BY: \n";
                              ptr->values[i].k = new_node->values[0].k;
                              ptr->values[i].delta = new_node->values[0].delta;
                              iout << "  k=" << ptr->values[i].k
                                   << "  n=" << ptr->values[i].n
                                   << "  delta=" << ptr->values[i].delta<< "\n";
                              replace=1;
                              break;
                            }
                          }

                          if (!replace)
                          {
                            ptr->multiplicity += 1;

                            if (ptr->multiplicity > MAX_MULTIPLICITY)
                            {
                              char err_msg[181];

                              sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
                                      ptr->multiplicity, MAX_MULTIPLICITY);
                              NAMD_die(err_msg);
                            }
                            if (!same_as_last) 
                              iout << "INCREASING MULTIPLICITY TO: " << ptr->multiplicity << "\n";

                            i= ptr->multiplicity - 1; 
                            ptr->values[i].k = new_node->values[0].k;
                            ptr->values[i].n = new_node->values[0].n;
                            ptr->values[i].delta = new_node->values[0].delta;

                            if (!same_as_last) 
                              iout << "  k=" << ptr->values[i].k
                                   << "  n=" << ptr->values[i].n
                                   << "  delta=" << ptr->values[i].delta<< "\n";
                          }
                        
                          iout << endi;
                        } 
                        //****** END CHARMM/XPLOR type changes

                        memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
                        delete new_node;

                        return;
                }

                ptr=ptr->next;
        }

        /*  CHARMM and XPLOR wildcards for dihedrals are luckily the same */
        /*  Check to see if we have any wildcards.  Since specific        */
        /*  entries are to take precedence, we'll put anything without  */
        /*  wildcards at the begining of the list and anything with     */
        /*  wildcards at the end of the list.  Then, we can just do a   */
        /*  linear search for a bond and be guaranteed to have specific */
        /*  entries take precendence over over wildcards                */
        if ( new_node->atom1wild ||
             new_node->atom2wild ||
             new_node->atom3wild ||
             new_node->atom4wild )
        {
                /*  add to the end of the list                                */
                tail->next=new_node;
                tail=new_node;

                memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
                return;
        }
        else
        {
                /*  add to the head of the list                                */
                new_node->next=dihedralp;
                dihedralp=new_node;

                memcpy(&last_dihedral, new_node, sizeof(dihedral_params));