Molecule.C

Go to the documentation of this file.

/*
   The class Molecule is used to hold all of the structural information
   for a simulation.  This information is read in from a .psf file and
   cross checked with the Parameters object passed in.  All of the structural
   information is then stored in arrays for use.
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>

#include "InfoStream.h"
#include "Molecule.h"
#include "strlib.h"
#include "MStream.h"
#include "Communicate.h"
// #include "Node.h"
#include "ObjectArena.h"
#include "Parameters.h"
#include "PDB.h"
#include "SimParameters.h"
#include "Hydrogen.h"
#include "UniqueSetIter.h"
#include "charm++.h"
/* BEGIN gf */
#include "GridForceGrid.h"
/* END gf */

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

#include "CompressPsf.h"
#include <deque>
#include <algorithm>
using namespace std;

#ifdef MEM_OPT_VERSION
template int lookupCstPool<AtomSignature>(const vector<AtomSignature>&, const AtomSignature&);
template int lookupCstPool<ExclusionSignature>(const vector<ExclusionSignature>&, const ExclusionSignature&);
#endif

class ResidueLookupElem
{
public:
  char mySegid[11];
  ResidueLookupElem *next;      // stored as a linked list
  int firstResid;               // valid resid is >= firstResid
  int lastResid;                // valid resid is <= lastResid
  ResizeArray<int> atomIndex;   // 0-based index for first atom in residue

  ResidueLookupElem(void) { next = 0; firstResid = -1; lastResid = -1; }
  ~ResidueLookupElem(void) { delete next; }
  int lookup(const char *segid, int resid, int *begin, int *end) const;
  ResidueLookupElem* append(const char *segid, int resid, int aid);
};

int ResidueLookupElem::lookup(
        const char *segid, int resid, int *begin, int *end) const {
    const ResidueLookupElem *elem = this;
    int rval = -1;  // error
    while ( elem && strcasecmp(elem->mySegid,segid) ) elem = elem->next;
    if ( elem && (resid >= elem->firstResid) && (resid <= elem->lastResid) ) {
      *begin = elem->atomIndex[resid - elem->firstResid];
      *end = elem->atomIndex[resid - elem->firstResid + 1];
      rval = 0;  // no error
    }
    return rval;
}

ResidueLookupElem* ResidueLookupElem::append(
        const char *segid, int resid, int aid) {
    ResidueLookupElem *rval = this;
    if ( firstResid == -1 ) {  // nothing added yet
      strcpy(mySegid,segid);
      firstResid = resid;
      lastResid = resid;
      atomIndex.add(aid);
      atomIndex.add(aid+1);
    } else if ( ! strcasecmp(mySegid,segid) ) {  // same segid
      if ( resid == lastResid ) {  // same resid
        atomIndex[lastResid - firstResid + 1] = aid + 1;
      } else if ( resid < lastResid ) {  // error
        // We can work around this by creating a new segment.
        iout << iWARN << "Residue " << resid <<
          " out of order in segment " << segid <<
          ", lookup for additional residues in this segment disabled.\n" << endi;
        rval = next = new ResidueLookupElem;
        next->append(segid,resid,aid);
      } else {  // new resid
        for ( ; lastResid < resid; ++lastResid ) atomIndex.add(aid);
        atomIndex[lastResid - firstResid + 1] = aid + 1;
      }
    } else {  // new segid
      rval = next = new ResidueLookupElem;
      next->append(segid,resid,aid);
    }
    return rval;
}


//  Lookup atom id from segment, residue, and name
int Molecule::get_atom_from_name(
        const char *segid, int resid, const char *aname) const {

  if (atomNames == NULL || resLookup == NULL)
  {
    NAMD_die("Tried to find atom from name on node other than node 0");
  }

  int i = 0;
  int end = 0;
  if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
  for ( ; i < end; ++i ) {
    #ifdef MEM_OPT_VERSION    
    Index idx = atomNames[i].atomnameIdx;
    if(!strcasecmp(aname, atomNamePool[idx])) return i;
    #else
    if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i;
    #endif
  }
  return -1;
}

//  Lookup number of atoms in residue from segment and residue
int Molecule::get_residue_size(
        const char *segid, int resid) const {

  if (atomNames == NULL || resLookup == NULL)
  {
    NAMD_die("Tried to find atom from name on node other than node 0");
  }
  int i = 0;
  int end = 0;
  if ( resLookup->lookup(segid,resid,&i,&end) ) return 0;
  return ( end - i );
}

//  Lookup atom id from segment, residue, and index in residue
int Molecule::get_atom_from_index_in_residue(
        const char *segid, int resid, int index) const {

  if (atomNames == NULL || resLookup == NULL)
  {
    NAMD_die("Tried to find atom from name on node other than node 0");
  }
  int i = 0;
  int end = 0;
  if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
  if ( index >= 0 && index < ( end - i ) ) return ( index + i );
  return -1;
}

/************************************************************************/
/*                  */
/*      FUNCTION initialize  */
/*                  */
/*  This is the initializer for the Molecule class.  It simply sets */
/*  the counts for all the various parameters to 0 and sets the pointers*/
/*  to the arrays that will store these parameters to NULL, since they  */
/*  have not been allocated yet.          */
/*                  */
/************************************************************************/

void Molecule::initialize(SimParameters *simParams, Parameters *param)
{
  if ( sizeof(int32) != 4 ) { NAMD_bug("sizeof(int32) != 4"); }
  this->simParams = simParams;
  this->params = param;
  /*  Initialize array pointers to NULL  */
  atoms=NULL;
  atomNames=NULL;
  resLookup=NULL;
  if ( simParams->globalForcesOn ) {
    resLookup = new ResidueLookupElem;
  }

  #ifdef MEM_OPT_VERSION
  eachAtomSig = NULL;
  atomSigPoolSize = 0;
  atomSigPool = NULL;
  massPoolSize = 0;
  atomMassPool = NULL;
  eachAtomMass = NULL;
  chargePoolSize = 0;
  atomChargePool = NULL;
  eachAtomCharge = NULL;
  #else
  bonds=NULL;
  angles=NULL;
  dihedrals=NULL;
  impropers=NULL;
  crossterms=NULL;
  #endif

  donors=NULL;
  acceptors=NULL;
  

  #ifndef MEM_OPT_VERSION      
  tmpArena=NULL;
  exclusions=NULL;
  bondsWithAtom=NULL;
  bondsByAtom=NULL;
  anglesByAtom=NULL;
  dihedralsByAtom=NULL;
  impropersByAtom=NULL;
  crosstermsByAtom=NULL;
  #endif

  #ifdef MEM_OPT_VERSION
  exclSigPool = NULL;
  exclChkSigPool = NULL;
  exclSigPoolSize = 0;
  eachAtomExclSig = NULL;
  #else
  exclusionsByAtom=NULL;
  fullExclusionsByAtom=NULL;
  modExclusionsByAtom=NULL;
  all_exclusions=NULL;
  #endif

  langevinParams=NULL;
  fixedAtomFlags=NULL;

  #ifdef MEM_OPT_VERSION  
  clusterSigs=NULL;
  #else
  cluster=NULL;
  clusterSize=NULL;
  #endif

  exPressureAtomFlags=NULL;
  rigidBondLengths=NULL;
  consIndexes=NULL;
  consParams=NULL;
  /* BEGIN gf */
  gridfrcIndexes=NULL;
  gridfrcParams=NULL;
  gridfrcGrid=NULL;
  /* END gf */
  stirIndexes=NULL;
  stirParams=NULL;
  movDragIndexes=NULL;
  movDragParams=NULL;
  rotDragIndexes=NULL;
  rotDragParams=NULL;
  consTorqueIndexes=NULL;
  consTorqueParams=NULL;
  consForceIndexes=NULL;
  consForce=NULL;
//fepb
  fepAtomFlags=NULL;
//fepe

  nameArena = new ObjectArena<char>;
  // nameArena->setAlignment(8);
  // arena->setAlignment(32);
  #ifndef MEM_OPT_VERSION
  arena = new ObjectArena<int32>;
  exclArena = new ObjectArena<char>;
  #endif
  // exclArena->setAlignment(32);

  /*  Initialize counts to 0 */
  numAtoms=0;
  numRealBonds=0;
  numBonds=0;
  numAngles=0;
  numDihedrals=0;
  numImpropers=0;
  numCrossterms=0;
  numDonors=0;
  numAcceptors=0;
  numExclusions=0;
  numLP=0;
  numConstraints=0;
  numStirredAtoms=0;
  numMovDrag=0;
  numRotDrag=0;
  numConsTorque=0;
  numConsForce=0;
  numFixedAtoms=0;
  numFixedGroups=0;
  numExPressureAtoms=0;
  numRigidBonds=0;
  numFixedRigidBonds=0;
  numMultipleDihedrals=0;
  numMultipleImpropers=0;
  numCalcBonds=0;
  numCalcAngles=0;
  numCalcDihedrals=0;
  numCalcImpropers=0;
  numCalcCrossterms=0;
  numCalcExclusions=0;

//fepb
  numFepInitial = 0;
  numFepFinal = 0;
//fepe

}

/*      END OF FUNCTION initialize */

/************************************************************************/
/*                  */
/*      FUNCTION Molecule        */
/*                  */
/*  This is the constructor for the Molecule class. */
/*                  */
/************************************************************************/

Molecule::Molecule(SimParameters *simParams, Parameters *param)
{
  initialize(simParams,param);
}

/************************************************************************/
/*                  */
/*      FUNCTION Molecule        */
/*                  */
/*  This is the constructor for the Molecule class from CHARMM/XPLOR files. */
/*                  */
/************************************************************************/

Molecule::Molecule(SimParameters *simParams, Parameters *param, char *filename, ConfigList *cfgList)
{
  initialize(simParams,param);

  if(simParams->useCompressedPsf)
      read_compressed_psf_file(filename, param);
  else if(simParams->genCompressedPsf){      
      compress_psf_file(this, filename, param, simParams, cfgList);
  }      
  else
      read_psf_file(filename, param);
}

/*      END OF FUNCTION Molecule      */

/************************************************************************/
/*                  */
/*        FUNCTION Molecule      */
/*                  */
/*  This is the destructor for the class Molecule.  It simply frees */
/*  the memory allocated for each of the arrays used to store the       */
/*  structure information.            */
/*                  */
/************************************************************************/

Molecule::~Molecule()
{
  /*  Check to see if each array was ever allocated.  If it was   */
  /*  then free it            */
  if (atoms != NULL)
    delete [] atoms;

  if (atomNames != NULL)
  {
    // subarrarys allocated from arena - automatically deleted
    delete [] atomNames;
  }
  delete nameArena;

  if (resLookup != NULL)
    delete resLookup;

  #ifdef MEM_OPT_VERSION
  if(eachAtomSig) delete [] eachAtomSig;
  if(atomSigPool) delete [] atomSigPool;
  #else
  if (bonds != NULL)
    delete [] bonds;

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

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

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

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

  if (exclusions != NULL)
    delete [] exclusions;
  #endif

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

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

  #ifdef MEM_OPT_VERSION
  if(exclSigPool) delete [] exclSigPool;
  if(exclChkSigPool) delete [] exclChkSigPool;
  if(eachAtomExclSig) delete [] eachAtomExclSig;
  #else
  if (bondsByAtom != NULL)
       delete [] bondsByAtom;
  
  if (anglesByAtom != NULL)
       delete [] anglesByAtom;
  
  if (dihedralsByAtom != NULL)
       delete [] dihedralsByAtom;
  
  if (impropersByAtom != NULL)
       delete [] impropersByAtom;
  
  if (crosstermsByAtom != NULL)
       delete [] crosstermsByAtom;  

  if (exclusionsByAtom != NULL)
       delete [] exclusionsByAtom;
  
  if (fullExclusionsByAtom != NULL)
       delete [] fullExclusionsByAtom;
  
  if (modExclusionsByAtom != NULL)
       delete [] modExclusionsByAtom;
  
  if (all_exclusions != NULL)
       delete [] all_exclusions;
  #endif


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

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


  #ifdef MEM_OPT_VERSION
  if(clusterSigs != NULL){      
      delete [] clusterSigs;
  }  
  #else
  if (cluster != NULL)
       delete [] cluster;
  if (clusterSize != NULL)
       delete [] clusterSize;
  #endif

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

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

//fepb
  if (fepAtomFlags != NULL)
       delete [] fepAtomFlags;
//fepe


  #ifndef MEM_OPT_VERSION
  delete arena;
  delete exclArena;
  #endif
}
/*      END OF FUNCTION Molecule      */

/************************************************************************/
/*                  */
/*        FUNCTION read_psf_file      */
/*                  */
/*   INPUTS:                */
/*  fname - Name of the .psf file to read        */
/*  params - pointer to Parameters object to use to obtain          */
/*     parameters for vdWs, bonds, etc.      */
/*                  */
/*  This function reads a .psf file in.  This is where just about   */
/*   all of the structural information for this class comes from.  The  */
/*   .psf file contains descriptions of the atom, bonds, angles,        */
/*   dihedrals, impropers, and exclusions.  The parameter object is     */
/*   used to look up parameters for each of these entities.    */
/*                  */
/************************************************************************/

void Molecule::read_psf_file(char *fname, Parameters *params)

{
#ifdef MEM_OPT_VERSION
    return;
#else
  char err_msg[512];  //  Error message for NAMD_die
  char buffer[512];  //  Buffer for file reading
  int i;      //  Loop counter
  int NumTitle;    //  Number of Title lines in .psf file
  FILE *psf_file;    //  pointer to .psf file
  int ret_code;    //  ret_code from NAMD_read_line calls

  /* Try and open the .psf file           */
  if ( (psf_file = Fopen(fname, "r")) == NULL)
  {
    sprintf(err_msg, "UNABLE TO OPEN .psf FILE %s", fname);
    NAMD_die(err_msg);
  }

  /*  Read till we have the first non-blank line of file    */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to see if we dropped out of the loop because of a     */
  /*  read error.  This shouldn't happen unless the file is empty */
  if (ret_code!=0)
  {
    sprintf(err_msg, "EMPTY .psf FILE %s", fname);
    NAMD_die(err_msg);
  }

  /*  The first non-blank line should contain the word "psf".    */
  /*  If we can't find it, die.               */
  if (!NAMD_find_word(buffer, "psf"))
  {
    sprintf(err_msg, "UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
       fname);
    NAMD_die(err_msg);
  }

  /*  Read until we find the next non-blank line      */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to see if we dropped out of the loop because of a     */
  /*  read error.  This shouldn't happen unless there is nothing  */
  /*  but the PSF line in the file        */
  if (ret_code!=0)
  {
    sprintf(err_msg, "MISSING EVERYTHING BUT PSF FROM %s", fname);
    NAMD_die(err_msg);
  }

  /*  This line should have the word "NTITLE" in it specifying    */
  /*  how many title lines there are        */
  if (!NAMD_find_word(buffer, "NTITLE"))
  {
    sprintf(err_msg,"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
       fname);
    NAMD_die(err_msg);
  }

  sscanf(buffer, "%d", &NumTitle);

  /*  Now skip the next NTITLE non-blank lines and then read in the*/
  /*  line which should contain NATOM        */
  i=0;

  while ( ((ret_code=NAMD_read_line(psf_file, buffer)) == 0) && 
    (i<NumTitle) )
  {
    if (!NAMD_blank_string(buffer))
      i++;
  }

  /*  Make sure we didn't exit because of a read error    */
  if (ret_code!=0)
  {
    sprintf(err_msg, "FOUND EOF INSTEAD OF NATOM IN PSF FILE %s", 
       fname);
    NAMD_die(err_msg);
  }

  while (NAMD_blank_string(buffer))
  {
    NAMD_read_line(psf_file, buffer);
  }

  /*  Check to make sure we have the line we want      */
  if (!NAMD_find_word(buffer, "NATOM"))
  {
    sprintf(err_msg, "DIDN'T FIND \"NATOM\" IN PSF FILE %s",
       fname);
    NAMD_die(err_msg);
  }

  /*  Read in the number of atoms, and then the atoms themselves  */
  sscanf(buffer, "%d", &numAtoms);

  read_atoms(psf_file, params);

  /*  Read until we find the next non-blank line      */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to make sure we didn't hit the EOF      */
  if (ret_code != 0)
  {
    NAMD_die("EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
  }

  /*  Look for the string "NBOND"          */
  if (!NAMD_find_word(buffer, "NBOND"))
  {
    NAMD_die("DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
  }

  /*  Read in the number of bonds and then the bonds themselves  */
  sscanf(buffer, "%d", &numBonds);

  if (numBonds)
    read_bonds(psf_file, params);

  /*  Read until we find the next non-blank line      */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to make sure we didn't hit the EOF      */
  if (ret_code != 0)
  {
    NAMD_die("EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
  }

  /*  Look for the string "NTHETA"        */
  if (!NAMD_find_word(buffer, "NTHETA"))
  {
    NAMD_die("DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
  }

  /*  Read in the number of angles and then the angles themselves */
  sscanf(buffer, "%d", &numAngles);

  if (numAngles)
    read_angles(psf_file, params);

  /*  Read until we find the next non-blank line      */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to make sure we didn't hit the EOF      */
  if (ret_code != 0)
  {
    NAMD_die("EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
  }

  /*  Look for the string "NPHI"          */
  if (!NAMD_find_word(buffer, "NPHI"))
  {
    NAMD_die("DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
  }

  /*  Read in the number of dihedrals and then the dihedrals      */
  sscanf(buffer, "%d", &numDihedrals);

  if (numDihedrals)
    read_dihedrals(psf_file, params);

  /*  Read until we find the next non-blank line      */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to make sure we didn't hit the EOF      */
  if (ret_code != 0)
  {
    NAMD_die("EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
  }

  /*  Look for the string "NIMPHI"        */
  if (!NAMD_find_word(buffer, "NIMPHI"))
  {
    NAMD_die("DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
  }

  /*  Read in the number of Impropers and then the impropers  */
  sscanf(buffer, "%d", &numImpropers);

  if (numImpropers)
    read_impropers(psf_file, params);

  /*  Read until we find the next non-blank line      */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to make sure we didn't hit the EOF      */
  if (ret_code != 0)
  {
    NAMD_die("EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
  }

  /*  Look for the string "NDON"        */
  if (!NAMD_find_word(buffer, "NDON"))
  {
    NAMD_die("DID NOT FIND NDON AFTER ATOM LIST IN PSF");
  }

  /*  Read in the number of hydrogen bond donors and then the donors */
  sscanf(buffer, "%d", &numDonors);

  if (numDonors)
    read_donors(psf_file);

  /*  Read until we find the next non-blank line      */
  ret_code = NAMD_read_line(psf_file, buffer);

  while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
  {
    ret_code = NAMD_read_line(psf_file, buffer);
  }

  /*  Check to make sure we didn't hit the EOF      */
  if (ret_code != 0)
  {
    NAMD_die("EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
  }

  /*  Look for the string "NACC"        */
  if (!NAMD_find_word(buffer, "NACC"))
  {
    NAMD_die("DID NOT FIND NACC AFTER ATOM LIST IN PSF");
  }

  /*  Read in the number of hydrogen bond donors and then the donors */
  sscanf(buffer, "%d", &numAcceptors);

  if (numAcceptors)
    read_acceptors(psf_file);

  /*  look for the explicit non-bonded exclusion section.     */
  while (!NAMD_find_word(buffer, "NNB"))
  {
    ret_code = NAMD_read_line(psf_file, buffer);

    if (ret_code != 0)
    {
      NAMD_die("EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
    }
  }

  /*  Read in the number of exclusions and then the exclusions    */
  sscanf(buffer, "%d", &numExclusions);

  if (numExclusions)
    read_exclusions(psf_file);

  /*  look for the cross-term section.     */
  int crossterms_present = 1;
  while (!NAMD_find_word(buffer, "NCRTERM"))
  {
    ret_code = NAMD_read_line(psf_file, buffer);

    if (ret_code != 0)
    {
      // hit EOF before finding cross-term section
      crossterms_present = 0;
      break;
    }
  }

  if ( crossterms_present) {

    /*  Read in the number of cross-terms and then the cross-terms*/
    sscanf(buffer, "%d", &numCrossterms);

    if (numCrossterms)
      read_crossterms(psf_file, params);

  }

  /*  Close the .psf file.  */
  Fclose(psf_file);

  //  analyze the data and find the status of each atom
  numRealBonds = numBonds;
  build_atom_status();

  return;
#endif
}
/*      END OF FUNCTION read_psf_file      */

/************************************************************************/
/*                  */
/*        FUNCTION read_compressed_psf_file      */
/*                  */
/*   INPUTS:                */
/*  fname - Name of the compressed .psf file to read        */
/*  params - pointer to Parameters object to use to obtain          */
/*     parameters for vdWs, bonds, etc.      */
/*                  */
/*  This function reads a compressed .psf file in.  This is where just about   */
/*   all of the structural information for this class comes from.  The  */
/*   .psf file contains descriptions of the atom, bonds, angles,        */
/*   dihedrals, impropers, and exclusions.  The parameter object is     */
/*   used to look up parameters for each of these entities.    */
/*                  */
/************************************************************************/
void Molecule::read_compressed_psf_file(char *fname, Parameters *params){
#ifndef MEM_OPT_VERSION
    return;
#else
    FILE *psf_file;    //  pointer to .psf file
    int ret_code;    //  ret_code from NAMD_read_line calls
    char buffer[512];

    /* Try and open the .psf file           */
    if ( (psf_file = Fopen(fname, "r")) == NULL)
    {
        char err_msg[512];
        sprintf(err_msg, "UNABLE TO OPEN THE COMPRESSED .psf FILE %s", fname);
        NAMD_die(err_msg);
    }     

    char strBuf[12];

    //Begin reading constant pools for atoms' basic information
    //1. segment names
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NSEGMENTNAMES"))
        NAMD_die("UNABLE TO FIND NSEGMENTNAMES");
    sscanf(buffer, "%d", &segNamePoolSize);
    if(segNamePoolSize!=0)
        segNamePool = new char *[segNamePoolSize];
    for(int i=0; i<segNamePoolSize; i++){
        NAMD_read_line(psf_file, buffer);
        sscanf(buffer, "%s", strBuf);
        segNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
        strcpy(segNamePool[i], strBuf);        
    }    

    //2. residue names
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NRESIDUENAMES"))
        NAMD_die("UNABLE TO FIND NRESIDUENAMES");
    sscanf(buffer, "%d", &resNamePoolSize);
    if(resNamePoolSize!=0)
        resNamePool = new char *[resNamePoolSize];
    for(int i=0; i<resNamePoolSize; i++){
        NAMD_read_line(psf_file, buffer);
        sscanf(buffer, "%s", strBuf);
        resNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
        strcpy(resNamePool[i], strBuf);
    }

    //3. atom names
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NATOMNAMES"))
        NAMD_die("UNABLE TO FIND NATOMNAMES");
    sscanf(buffer, "%d", &atomNamePoolSize);
    if(atomNamePoolSize!=0)
        atomNamePool = new char *[atomNamePoolSize];
    for(int i=0; i<atomNamePoolSize; i++){
        NAMD_read_line(psf_file, buffer);
        sscanf(buffer, "%s", strBuf);
        atomNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
        strcpy(atomNamePool[i], strBuf);
    }

    //4. atom types
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NATOMTYPES"))
        NAMD_die("UNABLE TO FIND NATOMTYPES");
    sscanf(buffer, "%d", &atomTypePoolSize);
    if(atomTypePoolSize!=0)
        atomTypePool = new char *[atomTypePoolSize];
    for(int i=0; i<atomTypePoolSize; i++){
        NAMD_read_line(psf_file, buffer);
        sscanf(buffer, "%s", strBuf);
        atomTypePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
        strcpy(atomTypePool[i], strBuf);
    }

    //5. charges
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NCHARGES"))
        NAMD_die("UNABLE TO FIND NCHARGES");
    sscanf(buffer, "%d", &chargePoolSize);
    if(chargePoolSize!=0)
        atomChargePool = new Real[chargePoolSize];
    for(int i=0; i<chargePoolSize; i++){
        NAMD_read_line(psf_file, buffer);
        sscanf(buffer, "%f", atomChargePool+i);        
    }

    //6. masses
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NMASSES"))
        NAMD_die("UNABLE TO FIND NMASSES");
    sscanf(buffer, "%d", &massPoolSize);
    if(massPoolSize!=0)
        atomMassPool = new Real[massPoolSize];
    for(int i=0; i<massPoolSize; i++){
        NAMD_read_line(psf_file, buffer);
        sscanf(buffer, "%f", atomMassPool+i);        
    }

    //Begin reading of atom signatures
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "ATOMSIGS"))
        NAMD_die("UNABLE TO FIND ATOMSIGS");
    sscanf(buffer, "%d", &atomSigPoolSize);
    atomSigPool = new AtomSignature[atomSigPoolSize];
    int typeCnt;
    int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
    int tisReal;
    int ttype;
    for(int i=0; i<atomSigPoolSize; i++){
        //BOND SIGS
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NBONDSIGS"))
            NAMD_die("UNABLE TO FIND NBONDSIGS");
        sscanf(buffer, "%d", &typeCnt);
        if(typeCnt!=0){
            atomSigPool[i].bondCnt = typeCnt;
            atomSigPool[i].bondSigs = new TupleSignature[typeCnt];
        }            
        for(int j=0; j<typeCnt; j++){
            NAMD_read_line(psf_file, buffer);            
            sscanf(buffer, "%d | %d | %d", &tmp1, &ttype, &tisReal);
            TupleSignature oneSig(1, BOND, (Index)ttype, (char)tisReal);
            oneSig.offset[0] = tmp1;            
            atomSigPool[i].bondSigs[j]=oneSig;
            if(tisReal) numRealBonds++;
        }

        //ANGLE SIGS
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NTHETASIGS"))
            NAMD_die("UNABLE TO FIND NTHETASIGS");
        sscanf(buffer, "%d", &typeCnt);
        if(typeCnt!=0){
            atomSigPool[i].angleCnt = typeCnt;
            atomSigPool[i].angleSigs = new TupleSignature[typeCnt];
        }            
        for(int j=0; j<typeCnt; j++){
            NAMD_read_line(psf_file, buffer);            
            sscanf(buffer, "%d %d | %d | %d", &tmp1, &tmp2, &ttype, &tisReal);
            TupleSignature oneSig(2,ANGLE,(Index)ttype, (char)tisReal);
            oneSig.offset[0] = tmp1;
            oneSig.offset[1] = tmp2;            
            atomSigPool[i].angleSigs[j] = oneSig;
        }
        //DIHEDRAL SIGS
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NPHISIGS"))
            NAMD_die("UNABLE TO FIND NPHISIGS");
        sscanf(buffer, "%d", &typeCnt);
        if(typeCnt!=0){
            atomSigPool[i].dihedralCnt = typeCnt;
            atomSigPool[i].dihedralSigs = new TupleSignature[typeCnt];
        }            
        for(int j=0; j<typeCnt; j++){
            NAMD_read_line(psf_file, buffer);            
            sscanf(buffer, "%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
            TupleSignature oneSig(3,DIHEDRAL,(Index)ttype, (char)tisReal);
            oneSig.offset[0] = tmp1;
            oneSig.offset[1] = tmp2;
            oneSig.offset[2] = tmp3;
            atomSigPool[i].dihedralSigs[j] = oneSig;
        }
        //IMPROPER SIGS
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NIMPHISIGS"))
            NAMD_die("UNABLE TO FIND NIMPHISIGS");
        sscanf(buffer, "%d", &typeCnt);
        if(typeCnt!=0){
            atomSigPool[i].improperCnt = typeCnt;
            atomSigPool[i].improperSigs = new TupleSignature[typeCnt];
        }            
        for(int j=0; j<typeCnt; j++){
            NAMD_read_line(psf_file, buffer);            
            sscanf(buffer, "%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
            TupleSignature oneSig(3,IMPROPER,(Index)ttype, (char)tisReal);
            oneSig.offset[0] = tmp1;
            oneSig.offset[1] = tmp2;
            oneSig.offset[2] = tmp3;
            atomSigPool[i].improperSigs[j] = oneSig;
        }        
        //CROSSTERM SIGS
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NCRTERMSIGS"))
            NAMD_die("UNABLE TO FIND NCRTERMSIGS");
        sscanf(buffer, "%d", &typeCnt);
        if(typeCnt!=0){
            atomSigPool[i].crosstermCnt = typeCnt;
            atomSigPool[i].crosstermSigs = new TupleSignature[typeCnt];
        }            
        for(int j=0; j<typeCnt; j++){
            NAMD_read_line(psf_file, buffer);            
            sscanf(buffer, "%d %d %d %d %d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &ttype, &tisReal);
            TupleSignature oneSig(7,CROSSTERM,(Index)ttype, (char)tisReal);
            oneSig.offset[0] = tmp1;
            oneSig.offset[1] = tmp2;
            oneSig.offset[2] = tmp3;
            oneSig.offset[3] = tmp4;
            oneSig.offset[4] = tmp5;
            oneSig.offset[5] = tmp6;
            oneSig.offset[6] = tmp7;
            atomSigPool[i].crosstermSigs[j] = oneSig;
        }        
    }

    //Reading exclusions' signatures    
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NEXCLSIGS")){
        NAMD_die("UNABLE TO FIND NEXCLSIGS");
    }
    sscanf(buffer, "%d", &exclSigPoolSize);
    if(exclSigPoolSize>0) exclSigPool = new ExclusionSignature[exclSigPoolSize];
    vector<int> fullExcls;
    vector<int> modExcls;
    for(int i=0; i<exclSigPoolSize; i++){
        int fullExclCnt = NAMD_read_int(psf_file, buffer);
        for(int j=0; j<fullExclCnt; j++)
            fullExcls.push_back(NAMD_read_int(psf_file, buffer));
        int modExclCnt = NAMD_read_int(psf_file, buffer);
        for(int j=0; j<modExclCnt; j++)
            modExcls.push_back(NAMD_read_int(psf_file, buffer));

        //The integers in both vectors should be in the increasing the order
        //since they have been sorted during the compression of psf file
        exclSigPool[i].setOffsets(fullExcls, modExcls);

        fullExcls.clear();
        modExcls.clear();
    }

    //Now begin to read atoms
    //This part could be read in the batch mode. All the above information
    //can be stored in memory which only occupies several KBs.
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NATOM"))
        NAMD_die("UNABLE TO FIND NATOM");
    sscanf(buffer, "%d", &numAtoms);
    
    if(numAtoms>0){
        atoms = new AtomCstInfo[numAtoms];
        atomNames = new AtomNameIdx[numAtoms];
        eachAtomMass = new Index[numAtoms];
        eachAtomCharge = new Index[numAtoms];
        eachAtomSig = new Index[numAtoms];
        eachAtomExclSig = new Index[numAtoms];

        clusterSigs = new int[numAtoms];

        hydrogenGroup.resize(0);
        ResidueLookupElem *tmpResLookup = resLookup;

        int residue_number; //for residue number
        char *segment_name;
        int read_count;
        int idx[10];
        

        for(int i=0; i<numAtoms; i++){
            NAMD_read_line(psf_file, buffer);
            read_count = sscanf(buffer, "%d %d %d %d %d %d %d %d %d %d",
                                idx, idx+1, idx+2, idx+3, idx+4,
                                idx+5, idx+6, idx+7, idx+8, idx+9);
            if(read_count!=10){
                char err_msg[128];
                sprintf(err_msg, "BAD ATOM LINE FORMAT IN COMPRESSED PSF FILE IN ATOM LINE %d\nLINE=%s",i+1, buffer);
                NAMD_die(err_msg);
            }

            segment_name = segNamePool[idx[0]];
            residue_number = idx[1];            
            atomNames[i].resnameIdx = (Index)idx[2];
            atomNames[i].atomnameIdx = (Index)idx[3];
            atomNames[i].atomtypeIdx = (Index)idx[4];
            eachAtomCharge[i] = (Index)idx[5];
            eachAtomMass[i] = (Index)idx[6];
            eachAtomSig[i] = (Index)idx[7];
            eachAtomExclSig[i] = (Index)idx[8];

            clusterSigs[i] = idx[9];

            //debugExclNum += (exclSigPool[idx[8]].fullExclCnt+exclSigPool[idx[8]].modExclCnt);

            //Add this atom to residue lookup table
            if(tmpResLookup) tmpResLookup =
                tmpResLookup->append(segment_name, residue_number, i);
            //Determine the type of the atom (H or O)
            Real thisAtomMass = atomMassPool[eachAtomMass[i]];
            if (thisAtomMass < 0.1) {
              atoms[i].status = LonepairAtom;
            } else if(thisAtomMass <= 3.5){
                atoms[i].status = HydrogenAtom;
            }else if(atomNamePool[atomNames[i].atomnameIdx][0]=='O' &&
                     (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)){
                atoms[i].status = OxygenAtom;
            }

            //Look up the vdw constants for this atom
            params->assign_vdw_index(atomTypePool[atomNames[i].atomtypeIdx],
                                     &(atoms[i]));
        }

        //printf("debugExclNum: %d\n", debugExclNum);
    }
    //build up information for numBonds/numDihedrals.. etc
    numBonds=0;
    numAngles=0;
    numDihedrals=0;
    numImpropers=0;
    numCrossterms=0;
    numTotalExclusions=0;
    for(int i=0; i<numAtoms; i++){
        AtomSignature *thisSig = &atomSigPool[eachAtomSig[i]];
        numBonds += thisSig->bondCnt;
        numAngles += thisSig->angleCnt;
        numDihedrals += thisSig->dihedralCnt;
        numImpropers += thisSig->improperCnt;
        numCrossterms += thisSig->crosstermCnt;

        ExclusionSignature *exclSig = &exclSigPool[eachAtomExclSig[i]];
        numTotalExclusions += (exclSig->fullExclCnt + exclSig->modExclCnt);
    }
    
    numTotalExclusions /= 2;

    //read extra bond parameters if there is an input of extra bonds (extraBondsOn is true)
    if(simParams->extraBondsOn){
        int numExtraParams=0;

        //read extra bond params
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NEXTRABONDPARAMS")){
            NAMD_die("UNABLE TO FIND NEXTRABONDPARAMS");
        }
        sscanf(buffer, "%d", &numExtraParams);
        if(numExtraParams>0){
            BondValue *newParams = new BondValue[params->NumBondParams+numExtraParams];
            memcpy(newParams, params->bond_array, params->NumBondParams*sizeof(BondValue));
            delete [] params->bond_array;            

            int curNumPs = params->NumBondParams;
            for(int i=0; i<numExtraParams; i++){
                Real k, x0;
                NAMD_read_line(psf_file, buffer);
                sscanf(buffer, "%f %f", &k, &x0);
                newParams[curNumPs+i].k = k;
                newParams[curNumPs+i].x0 = x0;
            }
            params->bond_array = newParams;
            params->NumBondParams += numExtraParams;
        }

        //read extra angle params
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NEXTRAANGLEPARAMS")){
            NAMD_die("UNABLE TO FIND NEXTRAANGLEPARAMS");
        }
        sscanf(buffer, "%d", &numExtraParams);
        if(numExtraParams>0){
            NAMD_die("CURRENTLY NOT SUPPORT EXTRA ANGLES");
        }

        //read extra diheral params
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NEXTRADIHEDRALPARAMS")){
            NAMD_die("UNABLE TO FIND NEXTRADIHEDRALPARAMS");
        }
        sscanf(buffer, "%d", &numExtraParams);
        if(numExtraParams>0){
            NAMD_die("CURRENTLY NOT SUPPORT EXTRA DIHEDRALS");
        }

        //read extra improper params
        NAMD_read_line(psf_file, buffer);
        if(!NAMD_find_word(buffer, "NEXTRAIMPROPERPARAMS")){
            NAMD_die("UNABLE TO FIND NEXTRAIMPROPERPARAMS");
        }
        sscanf(buffer, "%d", &numExtraParams);
        if(numExtraParams>0){
            NAMD_die("CURRENTLY NOT SUPPORT EXTRA IMPROPERS");
        }
    }

    //read DIHEDRALPARAMARRAY and IMPROPERPARAMARRAY    
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "DIHEDRALPARAMARRAY"))
        NAMD_die("UNABLE TO FIND DIHEDRALPARAMARRAY");
    for(int i=0; i<params->NumDihedralParams; i++){
        params->dihedral_array[i].multiplicity = NAMD_read_int(psf_file, buffer);
    }

    NAMD_read_line(psf_file, buffer); //to read a simple single '\n' line 
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "IMPROPERPARAMARRAY"))
        NAMD_die("UNABLE TO FIND IMPROPERPARAMARRAY");
    for(int i=0; i<params->NumImproperParams; i++){
        params->improper_array[i].multiplicity = NAMD_read_int(psf_file, buffer);
    }

    Fclose(psf_file);

    //numRealBonds is set when reading bondSignatures from compressed psf file
    
    build_atom_status();
#endif
}
/*      END OF FUNCTION read_compressed_psf_file      */

/************************************************************************/
/*                  */
/*        FUNCTION read_atoms      */
/*                  */
/*   INPUTS:                */
/*  fd - file pointer to the .psf file        */
/*  params - Parameters object to use for parameters    */
/*                  */
/*  this function reads in the Atoms section of the .psf file.      */
/*   This section consists of numAtoms lines that are of the form:  */
/*     <atom#> <mol> <seg#> <res> <atomname> <atomtype> <charge> <mass> */
/*   Each line is read into the appropriate entry in the atoms array.   */
/*   The parameters object is then used to determine the vdW constants  */
/*   for this atom.              */
/*                  */
/************************************************************************/

void Molecule::read_atoms(FILE *fd, Parameters *params)

{
#ifdef MEM_OPT_VERSION
    return;
#else
  char buffer[512];  // Buffer for reading from file
  int atom_number=0;  // Atom number 
  int last_atom_number=0; // Last atom number, used to assure
        // atoms are in order
  char segment_name[11]; // Segment name
  char residue_number[11]; // Residue number
  char residue_name[11];  // Residue name
  char atom_name[11];  // Atom name
  char atom_type[11];  // Atom type
  Real charge;    // Charge for the current atom
  Real mass;    // Mass for the current atom
  int read_count;    // Number of fields read by sscanf

  /*  Allocate the atom arrays          */
  atoms     = new Atom[numAtoms];
  atomNames = new AtomNameInfo[numAtoms];
  hydrogenGroup.resize(0);

  if (atoms == NULL || atomNames == NULL )
  {
    NAMD_die("memory allocation failed in Molecule::read_atoms");
  }

  ResidueLookupElem *tmpResLookup = resLookup;

  /*  Loop and read in numAtoms atom lines.      */
  while (atom_number < numAtoms)
  {
    /*  Get the line from the file        */
    NAMD_read_line(fd, buffer);

    /*  If its blank or a comment, skip it      */
    if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') )
      continue;

    /*  Parse up the line          */
    read_count=sscanf(buffer, "%d %s %s %s %s %s %f %f",
       &atom_number, segment_name, residue_number,
       residue_name, atom_name, atom_type, &charge, &mass);

    /*  Check to make sure we found what we were expecting  */
    if (read_count != 8)
    {
      char err_msg[128];

      sprintf(err_msg, "BAD ATOM LINE FORMAT IN PSF FILE IN ATOM LINE %d\nLINE=%s",
         last_atom_number+1, buffer);
      NAMD_die(err_msg);
    }

    /*  Check if this is in XPLOR format  */
    int atom_type_num;
    if ( sscanf(atom_type, "%d", &atom_type_num) > 0 )
    {
      NAMD_die("Structure (psf) file is either in CHARMM format (with numbers for atoms types, the X-PLOR format using names is required) or the segment name field is empty.");
    }

    /*  Make sure the atoms were in sequence    */
    if (atom_number != last_atom_number+1)
    {
      char err_msg[128];

      sprintf(err_msg, "ATOM NUMBERS OUT OF ORDER AT ATOM #%d OF PSF FILE",
         last_atom_number+1);
      NAMD_die(err_msg);
    }

    last_atom_number++;

    /*  Dynamically allocate strings for atom name, atom    */
    /*  type, etc so that we only allocate as much space    */
    /*  for these strings as we really need      */
    int reslength = strlen(residue_name)+1;
    int namelength = strlen(atom_name)+1;
    int typelength = strlen(atom_type)+1;

    atomNames[atom_number-1].resname = nameArena->getNewArray(reslength);
    atomNames[atom_number-1].atomname = nameArena->getNewArray(namelength);
    atomNames[atom_number-1].atomtype = nameArena->getNewArray(typelength);
  
    if (atomNames[atom_number-1].resname == NULL)
    {
      NAMD_die("memory allocation failed in Molecule::read_atoms");
    }

    /*  Put the values from this atom into the atoms array  */
    strcpy(atomNames[atom_number-1].resname, residue_name);
    strcpy(atomNames[atom_number-1].atomname, atom_name);
    strcpy(atomNames[atom_number-1].atomtype, atom_type);
    atoms[atom_number-1].mass = mass;
    atoms[atom_number-1].charge = charge;
    atoms[atom_number-1].status = UnknownAtom;

    /*  Add this atom to residue lookup table */
    if ( tmpResLookup ) tmpResLookup =
        tmpResLookup->append(segment_name, atoi(residue_number), atom_number-1);

    /*  Determine the type of the atom (H or O) */
    if (atoms[atom_number-1].mass < 0.1) {
      atoms[atom_number-1].status |= LonepairAtom;
    } else if (atoms[atom_number-1].mass <=3.5) {
      atoms[atom_number-1].status |= HydrogenAtom;
    } else if ((atomNames[atom_number-1].atomname[0] == 'O') && 
         (atoms[atom_number-1].mass >= 14.0) && 
         (atoms[atom_number-1].mass <= 18.0)) {
      atoms[atom_number-1].status |= OxygenAtom;
    }

    /*  Look up the vdw constants for this atom    */
    params->assign_vdw_index(atomNames[atom_number-1].atomtype, 
       &(atoms[atom_number-1]));
        }

  return;
#endif
}
/*      END OF FUNCTION read_atoms      */

/************************************************************************/
/*                  */
/*      FUNCTION read_bonds        */
/*                  */
/*  read_bonds reads in the bond section of the .psf file.  This    */
/*  section contains a list of pairs of numbers where each pair is      */
/*  represents two atoms that are bonded together.  Each atom pair is   */
/*  read in.  Then that parameter object is queried to determine the    */
/*  force constant and rest distance for the bond.      */
/*                  */
/************************************************************************/

void Molecule::read_bonds(FILE *fd, Parameters *params)

{
#ifdef MEM_OPT_VERSION
    return;
#else
  int atom_nums[2];  // Atom indexes for the bonded atoms
  char atom1name[11];  // Atom type for atom #1
  char atom2name[11];  // Atom type for atom #2
  register int j;      // Loop counter
  int num_read=0;    // Number of bonds read so far
  int origNumBonds = numBonds;   // number of bonds in file header

  /*  Allocate the array to hold the bonds      */
  bonds=new Bond[numBonds];

  if (bonds == NULL)
  {
    NAMD_die("memory allocations failed in Molecule::read_bonds");
  }

  /*  Loop through and read in all the bonds      */
  while (num_read < numBonds)
  {
    /*  Loop and read in the two atom indexes    */
    for (j=0; j<2; j++)
    {
      /*  Read the atom number from the file.         */
      /*  Subtract 1 to convert the index from the    */
      /*  1 to NumAtoms used in the file to the       */
      /*  0 to NumAtoms-1 that we need    */
      atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;

      /*  Check to make sure the index isn't too big  */
      if (atom_nums[j] >= numAtoms)
      {
        char err_msg[128];

        sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
        NAMD_die(err_msg);
      }
    }

    /*  Get the atom type for the two atoms.  When we query */
    /*  the parameter object, we need to send the atom type */
    /*  that is alphabetically first as atom 1.    */
    if (strcasecmp(atomNames[atom_nums[0]].atomtype, 
         atomNames[atom_nums[1]].atomtype) < 0)
    {
      strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
      strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
    }
    else
    {
      strcpy(atom2name, atomNames[atom_nums[0]].atomtype);
      strcpy(atom1name, atomNames[atom_nums[1]].atomtype);
    }

    /*  Query the parameter object for the constants for    */
    /*  this bond            */
    Bond *b = &(bonds[num_read]);
    params->assign_bond_index(atom1name, atom2name, b);

    /*  Assign the atom indexes to the array element  */
    b->atom1=atom_nums[0];
    b->atom2=atom_nums[1];

    /*  Make sure this isn't a fake bond meant for shake in x-plor.  */
    Real k, x0;
    params->get_bond_params(&k,&x0,b->bond_type);
    if ( k == 0. ) --numBonds;  // fake bond
    else ++num_read;  // real bond
  }

  /*  Tell user about our subterfuge  */
  if ( numBonds != origNumBonds ) {
    iout << iWARN << "Ignored " << origNumBonds - numBonds <<
            " bonds with zero force constants.\n" << endi;
    iout << iWARN <<
        "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
  }

  return;
#endif
}
/*      END OF FUNCTION read_bonds      */

/************************************************************************/
/*                  */
/*      FUNCTION read_angles        */
/*                  */
/*   INPUTS:                */
/*  fd - File descriptor for .psf file        */
/*  params - Parameters object to query for parameters    */
/*                  */
/*  read_angles reads the angle parameters from the .psf file.      */
/*   This section of the .psf file consists of a list of triplets of    */
/*   atom indexes.  Each triplet represents three atoms connected via   */
/*   an angle bond.  The parameter object is queried to obtain the      */
/*   constants for each bond.            */
/*                  */
/************************************************************************/

void Molecule::read_angles(FILE *fd, Parameters *params)

{
#ifdef MEM_OPT_VERSION
    return;
#else
  int atom_nums[3];  //  Atom numbers for the three atoms
  char atom1name[11];  //  Atom type for atom 1
  char atom2name[11];  //  Atom type for atom 2
  char atom3name[11];  //  Atom type for atom 3
  register int j;      //  Loop counter
  int num_read=0;    //  Number of angles read so far
  int origNumAngles = numAngles;  // Number of angles in file
  /*  Alloc the array of angles          */
  angles=new Angle[numAngles];

  if (angles == NULL)
  {
    NAMD_die("memory allocation failed in Molecule::read_angles");
  }

  /*  Loop through and read all the angles      */
  while (num_read < numAngles)
  {
    /*  Loop through the 3 atom indexes in the current angle*/
    for (j=0; j<3; j++)
    {
      /*  Read the atom number from the file.         */
      /*  Subtract 1 to convert the index from the    */
      /*  1 to NumAtoms used in the file to the       */
      /*  0 to NumAtoms-1 that we need    */
      atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;

      /*  Check to make sure the atom index doesn't   */
      /*  exceed the Number of Atoms      */
      if (atom_nums[j] >= numAtoms)
      {
        char err_msg[128];

        sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
        NAMD_die(err_msg);
      }
    }

    /*  Place the bond name that is alphabetically first  */
    /*  in the atom1name.  This is OK since the order of    */
    /*  atom1 and atom3 are interchangable.  And to search  */
    /*  the tree of angle parameters, we need the order     */
    /*  to be predictable.          */
    if (strcasecmp(atomNames[atom_nums[0]].atomtype, 
         atomNames[atom_nums[2]].atomtype) < 0)
    {
      strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
      strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
      strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
    }
    else
    {
      strcpy(atom1name, atomNames[atom_nums[2]].atomtype);
      strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
      strcpy(atom3name, atomNames[atom_nums[0]].atomtype);
    }

    /*  Get the constant values for this bond from the  */
    /*  parameter object          */
    params->assign_angle_index(atom1name, atom2name, 
       atom3name, &(angles[num_read]));

    /*  Assign the three atom indices      */
    angles[num_read].atom1=atom_nums[0];
    angles[num_read].atom2=atom_nums[1];
    angles[num_read].atom3=atom_nums[2];

    /*  Make sure this isn't a fake angle meant for shake in x-plor.  */
    Real k, t0, k_ub, r_ub;
    params->get_angle_params(&k,&t0,&k_ub,&r_ub,angles[num_read].angle_type);
    if ( k == 0. && k_ub == 0. ) --numAngles;  // fake angle
    else ++num_read;  // real angle
  }

  /*  Tell user about our subterfuge  */
  if ( numAngles != origNumAngles ) {
    iout << iWARN << "Ignored " << origNumAngles - numAngles <<
            " angles with zero force constants.\n" << endi;
  }

  return;
#endif
}
/*      END OF FUNCTION read_angles      */

/************************************************************************/
/*                  */
/*        FUNCTION read_dihedrals      */
/*                  */
/*   INPUTS:                */
/*  fd - file descriptor for the .psf file        */
/*  params - pointer to parameter object        */
/*                  */
/*  read_dihedreals reads the dihedral section of the .psf file.    */
/*   This section of the file contains a list of quartets of atom       */
/*   numbers.  Each quartet represents a group of atoms that form a     */
/*   dihedral bond.              */
/*                  */
/************************************************************************/

void Molecule::read_dihedrals(FILE *fd, Parameters *params)
{
#ifdef MEM_OPT_VERSION
    return;
#else
  int atom_nums[4];  // The 4 atom indexes
  int last_atom_nums[4];  // Atom numbers from previous bond
  char atom1name[11];  // Atom type for atom 1
  char atom2name[11];  // Atom type for atom 2
  char atom3name[11];  // Atom type for atom 3
  char atom4name[11];  // Atom type for atom 4
  register int j;      // loop counter
  int num_read=0;    // number of dihedrals read so far
  int multiplicity=1;  // multiplicity of the current bond
  Bool duplicate_bond;  // Is this a duplicate of the last bond
  int num_unique=0;   // Number of unique dihedral bonds

  //  Initialize the array used to check for duplicate dihedrals
  for (j=0; j<4; j++)
    last_atom_nums[j] = -1;

  /*  Allocate an array to hold the Dihedrals      */
  dihedrals = new Dihedral[numDihedrals];

  if (dihedrals == NULL)
  {
    NAMD_die("memory allocation failed in Molecule::read_dihedrals");
  }

  /*  Loop through and read all the dihedrals      */
  while (num_read < numDihedrals)
  {
    duplicate_bond = TRUE;

    /*  Loop through and read the 4 indexes for this bond   */
    for (j=0; j<4; j++)
    {
      /*  Read the atom number from the file.         */
      /*  Subtract 1 to convert the index from the    */
      /*  1 to NumAtoms used in the file to the       */
      /*  0 to NumAtoms-1 that we need    */
      atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;

      /*  Check for an atom index that is too large  */
      if (atom_nums[j] >= numAtoms)
      {
        char err_msg[128];

        sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
        NAMD_die(err_msg);
      }

      //  Check to see if this atom matches the last bond
      if (atom_nums[j] != last_atom_nums[j])
      {
        duplicate_bond = FALSE;
      }

      last_atom_nums[j] = atom_nums[j];
    }

    /*  Get the atom types for the 4 atoms so we can look  */
    /*  up the constants in the parameter object    */
    strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
    strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
    strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
    strcpy(atom4name, atomNames[atom_nums[3]].atomtype);

    //  Check to see if this is really a new bond or just
    //  a repeat of the last one
    if (duplicate_bond)
    {
      //  This is a duplicate, so increase the multiplicity
      multiplicity++;

      if (multiplicity == 2)
      {
        numMultipleDihedrals++;
      }
    }
    else
    {
      multiplicity=1;
      num_unique++;
    }

    /*  Get the constants for this dihedral bond    */
    params->assign_dihedral_index(atom1name, atom2name, 
       atom3name, atom4name, &(dihedrals[num_unique-1]),
       multiplicity);

    /*  Assign the atom indexes        */
    dihedrals[num_unique-1].atom1=atom_nums[0];
    dihedrals[num_unique-1].atom2=atom_nums[1];
    dihedrals[num_unique-1].atom3=atom_nums[2];
    dihedrals[num_unique-1].atom4=atom_nums[3];

    num_read++;
  }

  numDihedrals = num_unique;

  return;
#endif
}
/*      END OF FUNCTION read_dihedral      */

/************************************************************************/
/*                  */
/*        FUNCTION read_impropers      */
/*                  */
/*   INPUTS:                */
/*  fd - file descriptor for .psf file        */
/*  params - parameter object          */
/*                  */
/*  read_impropers reads the improper section of the .psf file.  */
/*   This section is identical to the dihedral section in that it is    */
/*   made up of a list of quartets of atom indexes that define the      */
/*   atoms that are bonded together.          */
/*                  */
/************************************************************************/

void Molecule::read_impropers(FILE *fd, Parameters *params)
{
#ifdef MEM_OPT_VERSION
    return;
#else
  int atom_nums[4];  //  Atom indexes for the 4 atoms
  int last_atom_nums[4];  //  Atom indexes from previous bond
  char atom1name[11];  //  Atom type for atom 1
  char atom2name[11];  //  Atom type for atom 2
  char atom3name[11];  //  Atom type for atom 3
  char atom4name[11];  //  Atom type for atom 4
  register int j;      //  Loop counter
  int num_read=0;    //  Number of impropers read so far
  int multiplicity=1;  // multiplicity of the current bond
  Bool duplicate_bond;  // Is this a duplicate of the last bond
  int num_unique=0;   // Number of unique dihedral bonds

  //  Initialize the array used to look for duplicate improper
  //  entries.  Set them all to -1 so we know nothing will match
  for (j=0; j<4; j++)
    last_atom_nums[j] = -1;

  /*  Allocate the array to hold the impropers      */
  impropers=new Improper[numImpropers];

  if (impropers == NULL)
  {
    NAMD_die("memory allocation failed in Molecule::read_impropers");
  }

  /*  Loop through and read all the impropers      */