Molecule.C

Go to the documentation of this file.

// Molecule.C is compiled twice!
// MOLECULE2_C undefined only compiles first half of file
// MOLECULE2_C defined only compiles second half of file
// This is shameful but it works.  Molecule needs refactoring badly.

/*
   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 "largefiles.h"  // must be first!

#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 "ComputeGridForce.h"
#include "GridForceGrid.h"

#include "MGridforceParams.h"
/* END gf */

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

#include "CompressPsf.h"
#include "ParallelIOMgr.h"
#include <deque>
#include <algorithm>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#ifndef GROMACS_PAIR
#define GROMACS_PAIR 1
#endif

#ifndef GROMACS_EXCLUSIONS
#define GROMACS_EXCLUSIONS 1
#endif

using namespace std;

#ifndef MOLECULE2_C  // first object file

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

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;

  // DRUDE
  is_lonepairs_psf = 1;  // anticipate having lone pair hosts
  is_drude_psf = 0;  // assume not Drude model
  drudeConsts=NULL;
  lphosts=NULL;  // might have lone pair hosts without Drude!
  anisos=NULL;
  tholes=NULL;
  lphostIndexes=NULL;
  // DRUDE

  //LCPO
  lcpoParamType = NULL;

  //for compressing molecule info
  atomSegResids=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;
  // JLai
  gromacsPairByAtom=NULL;
  // End of JLai
  // DRUDE
  tholesByAtom=NULL;
  anisosByAtom=NULL;
  // DRUDE
  #endif

  #ifdef MEM_OPT_VERSION
  exclSigPool = NULL;
  exclChkSigPool = NULL;
  exclSigPoolSize = 0;
  eachAtomExclSig = NULL;

  fixedAtomsSet = NULL;
  constrainedAtomsSet = 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;  
  #endif
  clusterSize=NULL;

  exPressureAtomFlags=NULL;
  rigidBondLengths=NULL;
  consIndexes=NULL;
  consParams=NULL;
  /* BEGIN gf */
  gridfrcIndexes=NULL;
  gridfrcParams=NULL;
  gridfrcGrid=NULL;
  numGridforces=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;
  alch_unpert_bonds = NULL;
  alch_unpert_angles = NULL;
  alch_unpert_dihedrals = NULL;
//fepe
//soluteScaling
  ssAtomFlags=NULL;
  ss_vdw_type=NULL;
  ss_index=NULL;
//soluteScaling
  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;
  // JLai
  numLJPair=0;
  // End of JLai
  numDonors=0;
  numAcceptors=0;
  numExclusions=0;

  // DRUDE
  numLonepairs=0;
  numDrudeAtoms=0;
  numTholes=0;
  numAnisos=0;
  numLphosts=0;
  numZeroMassAtoms=0;
  // DRUDE

  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;
  numCalcTholes=0;
  numCalcAnisos=0;
  numCalcCrossterms=0;
  numCalcExclusions=0;
  numCalcFullExclusions=0;
  // JLai
  numCalcLJPair=0;
  // End of JLai

//fepb
  numFepInitial = 0;
  numFepFinal = 0;
  num_alch_unpert_Bonds = 0;
  num_alch_unpert_Angles = 0;
  num_alch_unpert_Dihedrals = 0;
//fepe

  //fields related with pluginIO-based loading molecule structure
  occupancy = NULL;
  bfactor = NULL;

  qmElementArray=0;
  qmDummyElement=0;
  qmGrpSizes=0;
  qmAtomGroup=0;
  qmAtmChrg=0;
  qmAtmIndx=0;
  qmGrpID=0;
  qmGrpChrg=0;
  qmGrpMult=0;
  qmGrpNumBonds=0;
  qmMMBond=0;
  qmGrpBonds=0;
  qmMMBondedIndx=0;
  qmMMChargeTarget=0;
  qmMMNumTargs=0;
  qmDummyBondVal=0;
  qmMeMMindx=0;
  qmMeQMGrp=0;
  qmCustomPCIdxs=0;
  qmCustPCSizes=0;
  qmLSSSize=0;
  qmLSSIdxs=0;
  qmLSSMass=0;
  qmLSSRefIDs=0;
  qmLSSRefMass=0;
  qmLSSRefSize=0;
  qmNumBonds=0;
  cSMDindex=0;
  cSMDindxLen=0;
  cSMDpairs=0;
  cSMDKs=0;
  cSMDVels=0;
  cSMDcoffs=0;
  cSMDnumInst=0;
  
  goInit();
}

/*      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);

#ifdef MEM_OPT_VERSION
  if(simParams->useCompressedPsf)
      read_mol_signatures(filename, param, cfgList);
#else
        read_psf_file(filename, param); 
 //LCPO
  if (simParams->LCPOOn)
    assignLCPOTypes( 0 );
#endif      
 }

/************************************************************************/
/*                                                                      */
/*      FUNCTION Molecule                                               */
/*                                                                      */
/*  This is the constructor for the Molecule class from plugin IO.      */
/*                                                                      */
/************************************************************************/
Molecule::Molecule(SimParameters *simParams, Parameters *param, molfile_plugin_t *pIOHdl, void *pIOFileHdl, int natoms)
{
#ifdef MEM_OPT_VERSION
  NAMD_die("Sorry, plugin IO is not supported in the memory optimized version.");
#else
    initialize(simParams, param);
    numAtoms = natoms;
    int optflags = MOLFILE_BADOPTIONS;
    molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*sizeof(molfile_atom_t));
    memset(atomarray, 0, natoms*sizeof(molfile_atom_t));

    //1a. read basic atoms information
    int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
    if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
        free(atomarray);
        NAMD_die("ERROR: plugin failed reading structure data");
    }
    if(optflags == MOLFILE_BADOPTIONS) {
        free(atomarray);
        NAMD_die("ERROR: plugin didn't initialize optional data flags");
    }
    if(optflags & MOLFILE_OCCUPANCY) {
        setOccupancyData(atomarray);
    }
    if(optflags & MOLFILE_BFACTOR) {
        setBFactorData(atomarray);
    }
    //1b. load basic atoms information to the molecule object
    plgLoadAtomBasics(atomarray);    
    free(atomarray);

    //2a. read bonds
    //indices are one-based in read_bonds
    int *from, *to;
    float *bondorder;
    int *bondtype, nbondtypes;
    char **bondtypename;
    if(pIOHdl->read_bonds!=NULL) {
        if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
                                 &bondtype, &nbondtypes, &bondtypename)){
            NAMD_die("ERROR: failed reading bond information.");
        }
    }    
    //2b. load bonds information to the molecule object
    if(numBonds!=0) {
        plgLoadBonds(from,to);
    }

    //3a. read other bonded structures
    int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
    int ctermcols, ctermrows;
    int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
    int *impropertypes, numimpropertypes; 
    char **angletypenames, **dihedraltypenames, **impropertypenames;

    plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
    if(pIOHdl->read_angles!=NULL) {
        if(pIOHdl->read_angles(pIOFileHdl,
                  &numAngles, &plgAngles,
                  &angletypes, &numangletypes, &angletypenames,
                  &numDihedrals, &plgDihedrals,
                  &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
                  &numImpropers, &plgImpropers,
                  &impropertypes, &numimpropertypes, &impropertypenames,
                  &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
            NAMD_die("ERROR: failed reading angle information.");
        }
    }
    //3b. load other bonded structures to the molecule object
    if(numAngles!=0) plgLoadAngles(plgAngles);
    if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
    if(numImpropers!=0) plgLoadImpropers(plgImpropers);
    if(numCrossterms!=0) plgLoadCrossterms(plgCterms);

  numRealBonds = numBonds;
  build_atom_status();
  //LCPO
  if (simParams->LCPOOn)
    assignLCPOTypes( 2 );
#endif
}

/*      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;

  // DRUDE: free arrays read from PSF
  if (drudeConsts != NULL) delete [] drudeConsts;
  if (lphosts != NULL) delete [] lphosts;
  if (anisos != NULL) delete [] anisos;
  if (tholes != NULL) delete [] tholes;
  if (lphostIndexes != NULL) delete [] lphostIndexes;
  // DRUDE

  //LCPO
  if (lcpoParamType != NULL) delete [] lcpoParamType;

  #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;
  if(fixedAtomsSet) delete fixedAtomsSet;
  if(constrainedAtomsSet) delete constrainedAtomsSet;
  #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;

  // JLai
  if (gromacsPairByAtom != NULL)
      delete [] gromacsPairByAtom;
  // End of JLai

  // DRUDE
  if (tholesByAtom != NULL)
       delete [] tholesByAtom;
  if (anisosByAtom != NULL)
       delete [] anisosByAtom;
  // DRUDE
  #endif

  //LCPO
  if (lcpoParamType != NULL)
    delete [] lcpoParamType;

  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;  
  #endif
  if (clusterSize != NULL)
       delete [] clusterSize;

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

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

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

//soluteScaling
  if (ssAtomFlags != NULL)
       delete [] ssAtomFlags;
  if (ss_vdw_type != NULL)
       delete [] ss_vdw_type;
  if (ss_index != NULL)
       delete [] ss_index;
//soluteScaling

  if (qmAtomGroup != NULL)
       delete [] qmAtomGroup;
  
  if (qmAtmIndx != NULL)
       delete [] qmAtmIndx;
  
  if (qmAtmChrg != NULL)
       delete [] qmAtmChrg;
  
  
  if (qmGrpNumBonds != NULL)
       delete [] qmGrpNumBonds;
  
  if (qmGrpSizes != NULL)
       delete [] qmGrpSizes;
  
  if (qmDummyBondVal != NULL)
       delete [] qmDummyBondVal;
  
  if (qmMMNumTargs != NULL)
       delete [] qmMMNumTargs;
  
  if (qmGrpID != NULL)
       delete [] qmGrpID;
  
  if (qmGrpChrg != NULL)
       delete [] qmGrpChrg;
  
  if (qmGrpMult != NULL)
       delete [] qmGrpMult;
  
  if (qmMeMMindx != NULL)
       delete [] qmMeMMindx;
  
  if (qmMeQMGrp != NULL)
       delete [] qmMeQMGrp;
  
  if (qmLSSSize != NULL)
       delete [] qmLSSSize;
  
  if (qmLSSIdxs != NULL)
       delete [] qmLSSIdxs;
  
  if (qmLSSMass != NULL)
       delete [] qmLSSMass;
  
  if (qmLSSRefSize != NULL)
       delete [] qmLSSRefSize;
  
  if (qmLSSRefIDs != NULL)
       delete [] qmLSSRefIDs;
  
  if (qmLSSRefMass != NULL)
       delete [] qmLSSRefMass;
  
  if (qmMMBond != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
          if (qmMMBond[grpIndx] != NULL)
              delete [] qmMMBond[grpIndx];
      }
      delete [] qmMMBond;
  }
  
  if (qmGrpBonds != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
          if (qmGrpBonds[grpIndx] != NULL)
              delete [] qmGrpBonds[grpIndx];
      }
      delete [] qmGrpBonds;
  }
  
  if (qmMMBondedIndx != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
          if (qmMMBondedIndx[grpIndx] != NULL)
              delete [] qmMMBondedIndx[grpIndx];
      }
      delete [] qmMMBondedIndx;
  }
  
  if (qmMMChargeTarget != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
          if (qmMMChargeTarget[grpIndx] != NULL)
              delete [] qmMMChargeTarget[grpIndx];
      }
      delete [] qmMMChargeTarget;
  }
  
    if (qmElementArray != NULL){
        for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
            if (qmElementArray[atmInd] != NULL)
                delete [] qmElementArray[atmInd];
        }
        delete [] qmElementArray;
    }
    
    if (qmDummyElement != NULL){
        for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
            if (qmDummyElement[atmInd] != NULL)
                delete [] qmDummyElement[atmInd];
        }
        delete [] qmDummyElement;
    }

    if (qmCustPCSizes != NULL){
        delete [] qmCustPCSizes;
    }
    
    if (qmCustomPCIdxs != NULL){
        delete [] qmCustomPCIdxs;
    }
    
    if (cSMDindex != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
          if (cSMDindex[grpIndx] != NULL)
              delete [] cSMDindex[grpIndx];
      }
      delete [] cSMDindex;
    }
    
    if (cSMDpairs != NULL) {
      for(int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
          if (cSMDpairs[instIndx] != NULL)
              delete [] cSMDpairs[instIndx];
      }
      delete [] cSMDpairs;
    }
    
    if (cSMDindxLen != NULL)
        delete [] cSMDindxLen;
    if (cSMDKs != NULL)
        delete [] cSMDKs;
    if (cSMDVels != NULL)
        delete [] cSMDVels;
    if (cSMDcoffs != NULL)
        delete [] cSMDcoffs;
    
  #ifndef MEM_OPT_VERSION
  delete arena;
  delete exclArena;
  #endif
}
/*      END OF FUNCTION Molecule      */

#ifndef MEM_OPT_VERSION

//===Non-memory optimized version of functions that read Molecule file===//

/************************************************************************/
/*                  */
/*        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)
{
  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);
  }

  // DRUDE: set flag if we discover Drude PSF
  if (NAMD_find_word(buffer, "drude"))
  {
    if ( ! simParams->drudeOn ) {
      iout << iWARN << "Reading PSF supporting DRUDE without "
        "enabling the Drude model in the simulation config file\n" << endi;
    }
    is_drude_psf = 1;
  }
  // DRUDE

  /*  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);

  /*  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);

  //
  // The remaining sections that we can read might be optional to the PSF.
  //
  // Possibilities are:
  // - Drude simulation:
  //   + NUMLP (probably, but necessarily?)
  //   + NUMANISO (required)
  //   + NCRTERM (optional)
  // - non-Drude simulation
  //   + NUMLP (optional)
  //   + NCRTERM (optional)
  //
  // Also, there might be unrecognized PSF sections that we should skip.
  //

  // Keep reading until we recognize another section of PSF or find EOF.
  int is_found = 0;
  int is_found_numlp = 0;
  int is_found_numaniso = 0;
  int is_found_ncrterm = 0;
  while ( ! is_found ) {
    // Read until we find the next non-blank line
    do {
      ret_code = NAMD_read_line(psf_file, buffer);
    } while (ret_code == 0 && NAMD_blank_string(buffer) != 0 );
    // we either found EOF or we will try to match word in buffer
    if (ret_code != 0) {
      is_found = -1;  // found end of file
    }
    else if ( (is_found_numlp = NAMD_find_word(buffer, "NUMLP")) > 0) {
      is_found = is_found_numlp;
    }
    else if ( (is_found_numaniso = NAMD_find_word(buffer, "NUMANISO")) > 0) {
      is_found = is_found_numaniso;
    }
    else if ( (is_found_ncrterm = NAMD_find_word(buffer, "NCRTERM")) > 0) {
      is_found = is_found_ncrterm;
    }
  }

  if (is_found_numlp) {
    // We found lone pair hosts.
    // Read the number and then read the lone pair hosts.
    sscanf(buffer, "%d", &numLphosts);
  }

  if (numLphosts == 0) {
    // We either had no NUMLP section or we did and zero were listed.
    // Nevertheless, we have no lone pair hosts
    // so reset the simparams flag before simulation
    simParams->lonepairs = FALSE;
    is_lonepairs_psf = 0;
  }
  else if (simParams->lonepairs == FALSE /* but numLphosts > 0 */) {
    // Config file "lonepairs" option is (now) enabled by default.
    // Bad things will happen when lone pair hosts exist but "lonepairs"
    // in simparams is disabled. In this event, terminate with an error.
    NAMD_die("FOUND LONE PAIR HOSTS IN PSF WITH \"LONEPAIRS\" DISABLED IN CONFIG FILE");
  }

  // Process previously read bonds now that we know if lonepairs for TIP4 are explicit
  if (numBonds) process_bonds(params);

  if (is_found_numlp) {
    // Must follow process_bonds() since lone pairs are placed at end of bonds array
    if (numLphosts > 0) read_lphosts(psf_file);

    // Keep reading for next keyword.
    is_found = 0;
    is_found_numaniso = 0;
    is_found_ncrterm = 0;
    while ( ! is_found ) {
      // Read until we find the next non-blank line
      do {
        ret_code = NAMD_read_line(psf_file, buffer);
      } while (ret_code == 0 && NAMD_blank_string(buffer) != 0 );
      // we either found EOF or we will try to match word in buffer
      if (ret_code != 0) {
        is_found = -1;  // found end of file
      }
      else if ( (is_found_numaniso = NAMD_find_word(buffer, "NUMANISO")) > 0) {
        is_found = is_found_numaniso;
      }
      else if ( (is_found_ncrterm = NAMD_find_word(buffer, "NCRTERM")) > 0) {
        is_found = is_found_ncrterm;
      }
    }
  }

  if (is_found_numaniso && is_drude_psf) {
    // We are reading a Drude PSF and found the anisotropic terms.
    // Read the number and then read those terms.
    sscanf(buffer, "%d", &numAnisos);
    if (numAnisos > 0) read_anisos(psf_file);

    // Keep reading for next keyword.
    is_found = 0;
    is_found_ncrterm = 0;
    while ( ! is_found ) {
      // Read until we find the next non-blank line
      do {
        ret_code = NAMD_read_line(psf_file, buffer);
      } while (ret_code == 0 && NAMD_blank_string(buffer) != 0 );
      // we either found EOF or we will try to match word in buffer
      if (ret_code != 0) {
        is_found = -1;  // found end of file
      }
      else if ( (is_found_ncrterm = NAMD_find_word(buffer, "NCRTERM")) > 0) {
        is_found = is_found_ncrterm;
      }
    }
  }
  else if (is_drude_psf /* but not is_found_numaniso */) {
    NAMD_die("DID NOT FIND REQUIRED NUMANISO IN DRUDE PSF FILE");
  }
  else if (is_found_numaniso /* but not is_drude_file */) {
    NAMD_die("FOUND NUMANISO IN PSF FILE MISSING DRUDE DESIGNATION");
  }

  if (is_found_ncrterm) {
    // We found crossterms section of PSF.
    // Read the number and then read the crossterms.
    sscanf(buffer, "%d", &numCrossterms);
    if (numCrossterms > 0) read_crossterms(psf_file, params);
  }

  // Nothing else for us to read.

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

  if (is_drude_psf && numDrudeAtoms) {
    // Automatically build Drude bonds that were previously ignored.
    Bond *newbonds = new Bond[numBonds+numDrudeAtoms];
    memcpy(newbonds, bonds, numBonds*sizeof(Bond));
    delete [] bonds;
    bonds = newbonds;
    int origNumBonds = numBonds;
    for (i=0; i < numAtoms; i++) {
      if (!is_drude(i)) continue;
      Bond *b = &(bonds[numBonds++]);
      b->atom1 = i-1;
      b->atom2 = i;
      params->assign_bond_index(
          atomNames[i-1].atomtype, atomNames[i].atomtype, b
      );
    }
    if (numBonds-origNumBonds != numDrudeAtoms) {
      NAMD_die("must have same number of Drude particles and parents");
    }
  }

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

/************************************************************************/
/*                  */
/*        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)

{
  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];
  if(simParams->genCompressedPsf) {
      atomSegResids = new AtomSegResInfo[numAtoms];
  }

  // DRUDE: supplement Atom data
  if (is_drude_psf) {
    drudeConsts = new DrudeConst[numAtoms];
  }
  // DRUDE

  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)
  {
    // Standard PSF format has 8 columns:
    // ATOMNUM  SEGNAME  RESIDUE  RESNAME  ATOMNAME  ATOMTYPE  CHARGE  MASS

    /*  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);
    if (mass <= 0.05) ++numZeroMassAtoms;

    /*  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);
    }

    // DRUDE: read alpha and thole parameters from atom line
    if (is_drude_psf)
    {
      // Drude model PSF format has 11 columns, the 8 above plus 3 more:
      //   (unknown integer)  ALPHA  THOLE
      // These constants are used for the Thole interactions
      // (dipole interactions occurring between excluded non-bonded terms).

      Real alpha, thole;
      read_count=sscanf(buffer,
//          "%*d %*s %*s %*s %*s %*s %*f %*f %*d %*f %*f %f %f", &alpha, &thole);
                // the two columns preceding alpha and thole will disappear
          "%*d %*s %*s %*s %*s %*s %*f %*f %*d %f %f", &alpha, &thole);
      if (read_count != 2)
      {
        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);
      }
      drudeConsts[atom_number-1].alpha = alpha;
      drudeConsts[atom_number-1].thole = thole;
      if (fabs(alpha) >= 1e-6) ++numDrudeAtoms;
    }
    // DRUDE

    /*  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);

    if(atomSegResids) { //for compressing molecule information
        AtomSegResInfo *one = atomSegResids + (atom_number - 1);
        memcpy(one->segname, segment_name, strlen(segment_name)+1);
        one->resid = atoi(residue_number);
    }

    /*  Determine the type of the atom (H or O) */
    if ( simParams->ignoreMass ) {
    } else if (atoms[atom_number-1].mass <= 0.05) {
      atoms[atom_number-1].status |= LonepairAtom;
    } else if (atoms[atom_number-1].mass < 1.0) {
      atoms[atom_number-1].status |= DrudeAtom;
    } 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;
}
/*      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)

{
  int atom_nums[2];  // Atom indexes for the bonded atoms
  register int j;      // Loop counter
  int num_read=0;    // Number of bonds read so far

  /*  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);
      }
    }

    /*  Assign the atom indexes to the array element  */
    Bond *b = &(bonds[num_read]);
    b->atom1=atom_nums[0];
    b->atom2=atom_nums[1];

    ++num_read;
  }

}
/*      END OF FUNCTION read_bonds      */

/************************************************************************/
void Molecule::process_bonds(Parameters *params) {

  int atom_nums[2];  // Atom indexes for the bonded atoms
  int origNumBonds = numBonds;   // number of bonds in file header
  int num_read=0;    // Number of bonds read so far
  int numZeroFrcBonds = 0;
  int numLPBonds = 0;
  int numDrudeBonds = 0;

  for (int old_read=0; old_read < origNumBonds; ++old_read)
  {
    /*  Assign the atom indexes to the new array element  */
    Bond *b = &(bonds[num_read]);
    *b = bonds[old_read];
    atom_nums[0] = b->atom1;
    atom_nums[1] = b->atom2;

    /*  Query the parameter object for the constants for    */
    /*  this bond            */
    params->assign_bond_index(
        atomNames[atom_nums[0]].atomtype,
        atomNames[atom_nums[1]].atomtype,
        b);

    /*  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);

    Bool is_lp_bond = (is_lp(b->atom1) || is_lp(b->atom2));
    Bool is_drude_bond = (is_drude(b->atom1) || is_drude(b->atom2));
    numZeroFrcBonds += (k == 0.);
    numLPBonds += is_lp_bond;
    numDrudeBonds += is_drude_bond;
    /* haochuan: bonds involving lone pairs are still explicitly
     * defined in the parameter files of TIP4 water model (2020-06-29),
     * so we need to revert to the old behavior for including them in
     * numBonds.
     */
    switch (simParams->watmodel) {
      case WAT_TIP4: {
        // Drude force field does not use TIP4 water
        // so we can ignore is_drude_bond here
        if ( is_lonepairs_psf ) { // new format, has NUMLP section in psf, numLphosts > 0
          if (k == 0. || is_lp_bond) --numBonds;  // fake bond
          else ++num_read;  // real bond
        } else { // old format, no NUMLP section in psf
          numLPBonds -= is_lp_bond;
          if (k == 0. && !is_lp_bond) --numBonds;  // fake bond
          else ++num_read;  // real bond
        }
        break;
      }
      case WAT_SWM4:
        // intentionally fall through
      default: {
        // should be this the default behavior?
        if (k == 0. || is_lp_bond || is_drude_bond) --numBonds;  // fake bond
        else ++num_read;  // real bond
      }
    }
  }

  if ( num_read != numBonds ) {
    NAMD_bug("num_read != numBonds in Molecule::process_bonds()");
  }

  /*  Tell user about our subterfuge  */
  if ( numBonds != origNumBonds ) {
    if (numZeroFrcBonds) {
      iout << iWARN << "Ignored " << numZeroFrcBonds <<
          " bonds with zero force constants.\n" <<
          iWARN << "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
          endi;
    }
    if (numLPBonds) {
      iout << iWARN << "Ignored " << numLPBonds <<
          " bonds with lone pairs.\n" <<
          iWARN << "Will infer lonepair bonds from LPhost entries.\n" << endi;
    }
    if (numDrudeBonds) {
      iout << iWARN << "Ignored " << numDrudeBonds <<
          " bonds with Drude particles.\n" <<
          iWARN << "Will use polarizability to assign Drude bonds.\n" << endi;
    }
  }

}
/*      END OF FUNCTION process_bonds      */

/************************************************************************/
void Molecule::read_alch_unpert_bonds(FILE *fd) {
  int atom_nums[2];  // Atom indexes for the bonded atoms
  register int j;      // Loop counter
  int num_read=0;    // Number of bonds read so far

  alch_unpert_bonds=new Bond[num_alch_unpert_Bonds];

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

  while (num_read < num_alch_unpert_Bonds) {
    for (j=0; j<2; j++) {
      atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;

      if (atom_nums[j] >= numAtoms) {
        char err_msg[128];

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

    Bond *b = &(alch_unpert_bonds[num_read]);
    b->atom1=atom_nums[0];
    b->atom2=atom_nums[1];

    ++num_read;
  }
  return;
}
/*      END OF FUNCTION read_alch_unpert_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)

{
  int atom_nums[3];  //  Atom numbers for the three atoms
  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);
      }
    }

    /*  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];

    /*  Get the constant values for this bond from the  */
    /*  parameter object          */
    params->assign_angle_index(
      atomNames[atom_nums[0]].atomtype,
      atomNames[atom_nums[1]].atomtype,
      atomNames[atom_nums[2]].atomtype,
      &(angles[num_read]), simParams->alchOn ? -1 : 0);
    if ( angles[num_read].angle_type == -1 ) {
      iout << iWARN << "ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
           << endi;
    }

    /*  Make sure this isn't a fake angle meant for shake in x-plor.  */
    Real k, t0, k_ub, r_ub;
    if ( angles[num_read].angle_type == -1 ) { k = -1.;  k_ub = -1.; } else
    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;
}
/*      END OF FUNCTION read_angles      */

/************************************************************************/
void Molecule::read_alch_unpert_angles(FILE *fd) {
  int atom_nums[3];  //  Atom numbers for the three atoms
  register int j;    //  Loop counter
  int num_read=0;    //  Number of angles read so far

  alch_unpert_angles=new Angle[num_alch_unpert_Angles];

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

  while (num_read < num_alch_unpert_Angles) {
    for (j=0; j<3; j++) {
      atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;

      if (atom_nums[j] >= numAtoms) {
        char err_msg[128];
        sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN ALCH UNPERT PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
        NAMD_die(err_msg);
      }
    }

    alch_unpert_angles[num_read].atom1=atom_nums[0];
    alch_unpert_angles[num_read].atom2=atom_nums[1];
    alch_unpert_angles[num_read].atom3=atom_nums[2];

    ++num_read;
  }
  return;
}
/*      END OF FUNCTION read_alch_unpert_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)
{
  int atom_nums[4];  // The 4 atom indexes
  int last_atom_nums[4];  // Atom numbers from previous bond
  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];
    }

    //  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++;
    }

    /*  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];

    /*  Get the constants for this dihedral bond    */
    params->assign_dihedral_index(
       atomNames[atom_nums[0]].atomtype,
       atomNames[atom_nums[1]].atomtype,
       atomNames[atom_nums[2]].atomtype,
       atomNames[atom_nums[3]].atomtype,
       &(dihedrals[num_unique-1]),
       multiplicity, simParams->alchOn ? -1 : 0);
    if ( dihedrals[num_unique-1].dihedral_type == -1 ) {
      iout << iWARN << "ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
           << endi;
    }

    num_read++;
  }

  numDihedrals = num_unique;

  return;
}
/*      END OF FUNCTION read_dihedral      */

/*************************************************************************/
void Molecule::read_alch_unpert_dihedrals(FILE *fd) {
  int atom_nums[4];  // The 4 atom indexes
  int num_read=0;    // number of dihedrals read so far
  register int j;    // loop counter

  alch_unpert_dihedrals = new Dihedral[num_alch_unpert_Dihedrals];

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

  while (num_read < num_alch_unpert_Dihedrals) {
    for (j=0; j<4; j++) {
      atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;

      if (atom_nums[j] >= numAtoms) {
        char err_msg[128];

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

    alch_unpert_dihedrals[num_read].atom1=atom_nums[0];
    alch_unpert_dihedrals[num_read].atom2=atom_nums[1];
    alch_unpert_dihedrals[num_read].atom3=atom_nums[2];
    alch_unpert_dihedrals[num_read].atom4=atom_nums[3];

    num_read++;
  }
  return;
}
/*     END OF FUNCTION read_alch_unpert_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)
{
  int atom_nums[4];  //  Atom indexes for the 4 atoms
  int last_atom_nums[4];  //  Atom indexes from previous bond
  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      */
  while (num_read < numImpropers)
  {
    duplicate_bond = TRUE;

    /*  Loop through the 4 indexes for this improper  */
    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, "IMPROPERS")-1;

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

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

      if (atom_nums[j] != last_atom_nums[j])
      {
        duplicate_bond = FALSE;
      }

      last_atom_nums[j] = atom_nums[j];
    }

    //  Check to see if this is a duplicate improper
    if (duplicate_bond)
    {
      //  This is a duplicate improper.  So we don't
      //  really count this entry, we just update
      //  the parameters object
      multiplicity++;

      if (multiplicity == 2)
      {
        //  Count the number of multiples.
        numMultipleImpropers++;
      }
    }
    else
    {
      //  Not a duplicate
      multiplicity = 1;
      num_unique++;
    }

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

    /*  Look up the constants for this bond      */
    params->assign_improper_index(
       atomNames[atom_nums[0]].atomtype,
       atomNames[atom_nums[1]].atomtype,
       atomNames[atom_nums[2]].atomtype,
       atomNames[atom_nums[3]].atomtype,
       &(impropers[num_unique-1]),
       multiplicity);

    num_read++;
  }

  //  Now reset the numImpropers value to the number of UNIQUE
  //  impropers.  Sure, we waste a few entries in the improper_array
  //  on the master node, but it is very little space . . .
  numImpropers = num_unique;

  return;
}
/*      END OF FUNCTION read_impropers      */

/************************************************************************/
/*                  */
/*        FUNCTION read_crossterms      */
/*                  */
/*   INPUTS:                */
/*  fd - file descriptor for .psf file        */
/*  params - parameter object          */
/*                  */
/*   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_crossterms(FILE *fd, Parameters *params)

{
  int atom_nums[8];  //  Atom indexes for the 4 atoms
  int last_atom_nums[8];  //  Atom indexes from previous bond
  register int j;      //  Loop counter
  int num_read=0;    //  Number of items read so far
  Bool duplicate_bond;  // Is this a duplicate of the last bond

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

  /*  Allocate the array to hold the cross-terms */
  crossterms=new Crossterm[numCrossterms];

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

  /*  Loop through and read all the cross-terms      */
  while (num_read < numCrossterms)
  {
    duplicate_bond = TRUE;

    /*  Loop through the 8 indexes for this cross-term */
    for (j=0; j<8; 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, "CROSS-TERMS")-1;

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

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

      if (atom_nums[j] != last_atom_nums[j])
      {
        duplicate_bond = FALSE;
      }

      last_atom_nums[j] = atom_nums[j];
    }

    //  Check to see if this is a duplicate term
    if (duplicate_bond)
    {
      iout << iWARN << "Duplicate cross-term detected.\n" << endi;
    }

    /*  Assign the atom indexes        */
    crossterms[num_read].atom1=atom_nums[0];
    crossterms[num_read].atom2=atom_nums[1];
    crossterms[num_read].atom3=atom_nums[2];
    crossterms[num_read].atom4=atom_nums[3];
    crossterms[num_read].atom5=atom_nums[4];
    crossterms[num_read].atom6=atom_nums[5];
    crossterms[num_read].atom7=atom_nums[6];
    crossterms[num_read].atom8=atom_nums[7];

    /*  Look up the constants for this bond      */
    params->assign_crossterm_index(
       atomNames[atom_nums[0]].atomtype,
       atomNames[atom_nums[1]].atomtype,
       atomNames[atom_nums[2]].atomtype,
       atomNames[atom_nums[3]].atomtype,
       atomNames[atom_nums[4]].atomtype,
       atomNames[atom_nums[5]].atomtype,
       atomNames[atom_nums[6]].atomtype,
       atomNames[atom_nums[7]].atomtype,
       &(crossterms[num_read]));

    if(!duplicate_bond) num_read++;
  }

  numCrossterms = num_read;

  return;
}
/*      END OF FUNCTION read_impropers      */

/************************************************************************/
/*                  */
/*      FUNCTION read_donors        */
/*                  */
/*  read_donors 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 part of an H-bond.  Each atom pair is */
/*  read in.                                                            */
/*                  */
/*  Donor atoms are the heavy atoms to which hydrogens are bonded.      */
/*  There will always be a donor atom for each donor pair.  However,    */
/*  for a united-atom model there may not be an explicit hydrogen       */
/*  present, in which case the second atom index in the pair will be    */
/*  given as 0 in the PSF (and stored as -1 in this program's internal  */
/*  storage).                                                           */
/************************************************************************/

void Molecule::read_donors(FILE *fd)
{
  int d[2];               // temporary storage of donor atom index
  register int j;      // Loop counter
  int num_read=0;    // Number of bonds read so far
  int num_no_hydr=0;      // Number of bonds with no hydrogen given

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

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

  /*  Loop through and read in all the donors      */
  while (num_read < numDonors)
  {
    /*  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    */
      d[j]=NAMD_read_int(fd, "DONORS")-1;

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

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

      /*  Check if there is a hydrogen given */
      if (d[j] < 0)
                          num_no_hydr++;
    }

    /*  Assign the atom indexes to the array element  */
    Bond *b = &(donors[num_read]);
    b->atom1=d[0];
    b->atom2=d[1];

    num_read++;
  }

  return;
}
/*      END OF FUNCTION read_donors      */


/************************************************************************/
/*                  */
/*      FUNCTION read_acceptors        */
/*                  */
/*  read_acceptors 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 part of an H-bond.  Each atom pair is */
/*  read in.                                                            */
/*                  */
/*  Acceptor atoms are the heavy atoms to which hydrogens directly      */
/*  orient in a hydrogen bond interaction.  There will always be an     */
/*  acceptor atom for each acceptor pair.  The antecedent atom, to      */
/*  which the acceptor is bound, may not be given in the structure,     */
/*  however, in which case the second atom index in the pair will be    */
/*  given as 0 in the PSF (and stored as -1 in this program's internal  */
/*  storage).                                                           */
/************************************************************************/

void Molecule::read_acceptors(FILE *fd)
{
  int d[2];               // temporary storage of atom index
  register int j;      // Loop counter
  int num_read=0;    // Number of bonds read so far
        int num_no_ante=0;      // number of pairs with no antecedent

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

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

  /*  Loop through and read in all the acceptors      */
  while (num_read < numAcceptors)
  {
    /*  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    */
      d[j]=NAMD_read_int(fd, "ACCEPTORS")-1;

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

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

      /*  Check if there is an antecedent given */
      if (d[j] < 0)
                          num_no_ante++;
    }

    /*  Assign the atom indexes to the array element  */
    Bond *b = &(acceptors[num_read]);
    b->atom1=d[0];
    b->atom2=d[1];

    num_read++;
  }

  return;
}
/*      END OF FUNCTION read_acceptors      */


/************************************************************************/
/*                  */
/*      FUNCTION read_exclusions      */
/*                  */
/*   INPUTS:                */
/*  fd - file descriptor for .psf file        */
/*                  */
/*  read_exclusions reads in the explicit non-bonded exclusions     */
/*  from the .psf file.  This section is a little funky, so hang on.    */
/*  Ok, first there is a list of atom indexes that is NumExclusions     */
/*  long.  These are in some sense the atoms that will be exlcuded.     */
/*  Following this list is a list of NumAtoms length that is a list     */
/*  of indexes into the list of excluded atoms.  So an example.  Suppose*/
/*  we have a 5 atom simulation with 3 explicit exclusions.  The .psf   */
/*  file could look like:            */
/*                  */
/*  3!NNB                */
/*  3 4 5                */
/*  0 1 3 3 3              */
/*                  */
/*  This would mean that atom 1 has no explicit exclusions.  Atom 2     */
/*  has an explicit exclusion with atom 3.  Atom 3 has an explicit      */
/*  exclusion with atoms 4 AND 5.  And atoms 4 and 5 have no explicit   */
/*  exclusions.  Got it!?!  I'm not sure who dreamed this up . . .      */
/*                  */
/************************************************************************/

void Molecule::read_exclusions(FILE *fd)

{
  int *exclusion_atoms;  //  Array of indexes of excluded atoms
  register int num_read=0;    //  Number fo exclusions read in
  int current_index;  //  Current index value
  int last_index;    //  the previous index value
  register int insert_index=0;  //  index of where we are in exlcusions array

  /*  Allocate the array of exclusion structures and the array of */
  /*  exlcuded atom indexes          */
  exclusions      = new Exclusion[numExclusions];
  exclusion_atoms = new int[numExclusions];

  if ( (exclusions == NULL) || (exclusion_atoms == NULL) )
  {
    NAMD_die("memory allocation failed in Molecule::read_exclusions");
  }

  /*  First, read in the excluded atoms list      */
  for (num_read=0; num_read<numExclusions; num_read++)
  {
    /*  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    */
    exclusion_atoms[num_read]=NAMD_read_int(fd, "IMPROPERS")-1;

    /*  Check for an illegal index        */
    if (exclusion_atoms[num_read] >= numAtoms)
    {
      char err_msg[128];

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

  /*  Now, go through and read the list of NumAtoms pointers into */
  /*  the array that we just read in        */
  last_index=0;

  for (num_read=0; num_read<numAtoms; num_read++)
  {
    /*  Read in the current index value      */
    current_index=NAMD_read_int(fd, "EXCLUSIONS");

    /*  Check for an illegal pointer      */
    if (current_index>numExclusions)
    {
      char err_msg[128];

      sprintf(err_msg, "EXCLUSION INDEX %d LARGER THAN NUMBER OF EXLCUSIONS %d IN PSF FILE, EXCLUSION #%d\n", 
         current_index+1, numExclusions, num_read);
      NAMD_die(err_msg);
    }

    /*  Check to see if it matches the last index.  If so   */
    /*  than this atom has no exclusions.  If not, then     */
    /*  we have to build some exclusions      */
    if (current_index != last_index)
    {
      /*  This atom has some exclusions.  Loop from   */
      /*  the last_index to the current index.  This  */
      /*  will include how ever many exclusions this  */
      /*  atom has          */
      for (insert_index=last_index; 
           insert_index<current_index; insert_index++)
      {
        /*  Assign the two atoms involved.      */
        /*  The first one is our position in    */
        /*  the list, the second is based on    */
        /*  the pointer into the index list     */
        int a1 = num_read;
        int a2 = exclusion_atoms[insert_index];
        if ( a1 < a2 ) {
          exclusions[insert_index].atom1 = a1;
          exclusions[insert_index].atom2 = a2;
        } else if ( a2 < a1 ) {
          exclusions[insert_index].atom1 = a2;
          exclusions[insert_index].atom2 = a1;
        } else {
          char err_msg[128];
          sprintf(err_msg, "ATOM %d EXCLUDED FROM ITSELF IN PSF FILE\n", a1+1);
          NAMD_die(err_msg);
        }
      }

      last_index=current_index;
    }
  }

  /*  Free our temporary list of indexes        */
  delete [] exclusion_atoms;

  return;
}
/*      END OF FUNCTION read_exclusions      */

/************************************************************************/
/*      FUNCTION read_exclusions      */
/*                  */
/*   INPUTS:                */
/*   int* atom_i - array of atom i indices    */
/*   int* atom_j - array of atom j indices    */
/*   int num_exclusion - length of array           */
/*                  */
/* JLai August 16th, 2012 */
/************************************************************************/
void Molecule::read_exclusions(int* atom_i, int* atom_j, int num_exclusion)
{
    /*  Allocate the array of exclusion structures and the array of */
  /*  exlcuded atom indexes          */
  exclusions       = new Exclusion[num_exclusion];
  int loop_counter = 0;  
  int a=0;
  int b=0;

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

  /* The following code only guarantees that exclusion.atom1 is < exclusion.atom2 */
  for (loop_counter = 0; loop_counter < num_exclusion; loop_counter++) {
        
        if ( (atom_i == NULL) || (atom_j == NULL) ) {
          NAMD_die("null pointer expection in Molecule::read_exclusions");
        }

        a = atom_i[loop_counter];
        b = atom_j[loop_counter];
        if(a < b) {
                exclusions[loop_counter].atom1 = a;
                exclusions[loop_counter].atom2 = b;
        } else {
                exclusions[loop_counter].atom1 = b;
                exclusions[loop_counter].atom2 = a;
        }
        exclusionSet.add(Exclusion(exclusions[loop_counter].atom1,exclusions[loop_counter].atom2));
  }

  if ( ! CkMyPe() ) {
    iout << iINFO << "ADDED " << num_exclusion << " EXPLICIT EXCLUSIONS: THIS VALUE WILL *NOT* BE ADDED TO THE STRUCTURE SUMMARY\n" << endi;
  }

   return;
}
/*      END OF FUNCTION read_exclusions      */

/************************************************************************
 *
 *        FUNCTION read_lphosts
 *
 *   INPUTS:
 *  fd - file pointer to the .psf file
 *
 *  This function reads in the lone pair host section of the .psf file.
 *  Each lone pair is supported by 2 or 3 host atoms.
 *
 *  All lonepair specifications in a PSF are expected to have three
 *  associated values.  However, the meaning of these values depends
 *  on the lonepair type.  Nonetheless, for simplicity with the old
 *  code, the struct values are still called "distance", "angle",
 *  and "dihedral", even when this is not what the value signifies.
 *
 ************************************************************************/
void Molecule::read_lphosts(FILE *fd)
{
  char buffer[512];  // Buffer for reading from file
  char weight[8];    // Weighting type identified by string --
                     // unused/unsupported outside of CHARMM RTF
                     // so we read it, but ignore it.
  int numhosts;  // Refers to the number of atoms that support the
                 // given lone pair, either 2 or 3.
  int index;  // 1-based index into the stream of numbers (8 per line)
              // that indicates the start of each LP host record.
  int next_index = 1;  // Forecast next index value as an error check.
  int i, read_count;
  Real value1, value2, value3; 

  // called only if numLphosts > 0
  lphosts = new Lphost[numLphosts];
  if (lphosts == NULL) {
    NAMD_die("memory allocation failed in Molecule::read_lphosts");
  }
  for (i = 0;  i < numLphosts;  i++) {
    NAMD_read_line(fd, buffer);
    if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
    read_count=sscanf(buffer, "%d %d %6s %f %f %f",
        &numhosts, &index, weight, &value1, &value2, &value3);
    // The "weight" specification in PSF remains hardcoded as "F"
    // (for false) inside NAMD.  Again, no extant force field uses
    // lonepairs with weighted placement, so this can really only be
    // specified inside an RTF, which NAMD never reads anyway.
    if (read_count != 6 || index != next_index ||
        strcmp(weight,"F") != 0 || numhosts < 2 || 3 < numhosts) {
      char err_msg[128];
      sprintf(err_msg, "BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n"
          "LINE=%s\n", i+1, buffer);
      NAMD_die(err_msg);
    }
    lphosts[i].numhosts = numhosts; // currently must be 2 or 3
    next_index += numhosts + 1;  // add 1 to account for LP index
    if (numhosts == 2) {
      lphosts[i].distance = value1;
      lphosts[i].angle = value2;
      lphosts[i].dihedral = 0.0;  // ignore value3
    }
    else {  // numhosts == 3
      lphosts[i].distance = value1;
      lphosts[i].angle = value2 * (M_PI/180);     // convert to radians
      lphosts[i].dihedral = value3 * (M_PI/180);  // convert to radians
    }
  }
  // Resize bonds to accommodate the lonepairs.
  Bond *newbonds = new Bond[numBonds+numLphosts];
  memcpy(newbonds, bonds, numBonds*sizeof(Bond));
  delete [] bonds;
  bonds = newbonds;
  for (i = 0;  i < numLphosts;  i++) {
    // Subtract 1 to get 0-based atom index
    lphosts[i].atom1 = NAMD_read_int(fd, "LPHOSTS")-1;
    lphosts[i].atom2 = NAMD_read_int(fd, "LPHOSTS")-1;
    lphosts[i].atom3 = NAMD_read_int(fd, "LPHOSTS")-1;
    // For numhosts==2, set unused atom4 to atom1
    lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
        NAMD_read_int(fd, "LPHOSTS")-1 : lphosts[i].atom1 );
    // Add dummy bond entry for connectivity.
    Bond *b = &(bonds[numBonds++]);
    b->atom1 = lphosts[i].atom1;
    b->atom2 = lphosts[i].atom2;
    b->bond_type = -1; // dummy index, never used
  }
}
/*      END OF FUNCTION read_lphosts    */

/************************************************************************/
/*                  */
/*        FUNCTION read_anisos     */
/*                  */
/*   INPUTS:                */
/*  fd - file pointer to the .psf file        */
/*                  */
/*  this function reads in the anisotropic terms section of .psf file. */
/*                  */
void Molecule::read_anisos(FILE *fd)
{
  char buffer[512];  // Buffer for reading from file
  int numhosts, index, i, read_count;
  Real k11, k22, k33;

  anisos = new Aniso[numAnisos];
  if (anisos == NULL)
  {
    NAMD_die("memory allocation failed in Molecule::read_anisos");
  }
  for (i = 0;  i < numAnisos;  i++)
  {
    NAMD_read_line(fd, buffer);
    if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
    read_count=sscanf(buffer, "%f %f %f", &k11, &k22, &k33);
    if (read_count != 3)
    {
      char err_msg[128];
      sprintf(err_msg, "BAD FORMAT FOR ANISO LINE %d IN PSF FILE LINE\n"
          "LINE=%s\n", i+1, buffer);
      NAMD_die(err_msg);
    }
    anisos[i].k11 = k11;
    anisos[i].k22 = k22;
    anisos[i].k33 = k33;
  }
  for (i = 0;  i < numAnisos;  i++) {
    anisos[i].atom1 = NAMD_read_int(fd, "ANISOS")-1;
    anisos[i].atom2 = NAMD_read_int(fd, "ANISOS")-1;
    anisos[i].atom3 = NAMD_read_int(fd, "ANISOS")-1;
    anisos[i].atom4 = NAMD_read_int(fd, "ANISOS")-1;
  }
}
/*      END OF FUNCTION read_anisos     */

//LCPO
inline int getLCPOTypeAmber(char atomType[11], int numBonds) {

  //Hydrogen
  if (atomType[0] == 'H' || atomType[0] == 'h') {
    return 0;

  //Carbon
  } else if (atomType[0] == 'C' || atomType[0] == 'c') {
    if (//Sp3 Carbon
      //atomType[1] == 'T')// ||
      strcmp(atomType, "CT" )==0 )
      //strcmp(atomType, "CP1" )==0 ||
      //strcmp(atomType, "CP2" )==0 ||
      //strcmp(atomType, "CP3" )==0 ||
      //strcmp(atomType, "CS"  )==0 )
      {
      if (numBonds == 1)
        return 1;
      else if (numBonds == 2)
        return 2;
      else if (numBonds == 3)
        return 3;
      else if (numBonds == 4)
        return 4;
      else
        return 1;

    } else {//Sp2 or other
      if (numBonds == 2)
        return 5;
      else if (numBonds == 3)
        return 6;
      else
        return 1;
    }

  //Nitrogen
  } else if (atomType[0] == 'N' || atomType[0] == 'n') {
    if ( strcmp(atomType, "N3"  ) == 0 ) { //Sp3 Nitrogen
      if (numBonds == 1)
        return 11;
      else if (numBonds == 2)
        return 12;
      else if (numBonds == 3)
        return 13;
      else
        return 11;

    } else {//SP2 Nitrogen
      if (numBonds == 1)
        return 14;
      else if (numBonds == 2)
        return 15;
      else if (numBonds == 3)
        return 16;
      else
        return 11;
    }

  //Oxygen
  } else if (atomType[0] == 'O' || atomType[0] == 'o') {

    if ( strcmp(atomType, "O" )==0) {//Sp2 Oxygen
      return 9;
    } else if (strcmp(atomType, "O2" )==0) {//Carboxylate Oxygen
      return 10;
    } else { // Sp3 Oxygen
      if (numBonds == 1)
        return 7;
      else if (numBonds == 2)
        return 8;
      else
        return 7;
    }

  //Sulfur
  } else if (atomType[0] == 'S' || atomType[0] == 's') {
    if ( strcmp(atomType, "SH" )==0) { //Sulfur 1 neighbor
      return 17;
    } else {
      return 18;
    }

  //Phosphorus
  } else if (atomType[0] == 'P' || atomType[0] == 'p') {
      if (numBonds == 3)
        return 19;
      else if (numBonds == 4)
        return 20;
      else
        return 19;
  } else if (atomType[0] == 'Z') { // ? just to agree with Amber mdread.f
    return 0;
  } else  if ( strcmp(atomType, "MG" )==0) { //Mg
    return 22;
  } else { // unknown atom type
    return 5;
  }
  return 5;
} // getLCPOTypeAmber

inline int getLCPOTypeCharmm(char atomType[11], int numBonds) {

  //Hydrogen
  if (atomType[0] == 'H') {
    return 0;

  //Carbon
  } else if (atomType[0] == 'C') {
    if (//Sp3 Carbon
      atomType[1] == 'T' ||
      strcmp(atomType, "CP1" )==0 ||
      strcmp(atomType, "CP2" )==0 ||
      strcmp(atomType, "CP3" )==0 ||
      strcmp(atomType, "CS"  )==0 ) {
      if (numBonds == 1)
        return 1;
      else if (numBonds == 2)
        return 2;
      else if (numBonds == 3)
        return 3;
      else if (numBonds == 4)
        return 4;
      else
        return 1;

    } else if (//Sp2
      strcmp(atomType, "C"   )==0 ||
      strcmp(atomType, "CA"  )==0 ||
      strcmp(atomType, "CC"  )==0 ||
      strcmp(atomType, "CD"  )==0 ||
      strcmp(atomType, "CN"  )==0 ||
      strcmp(atomType, "CY"  )==0 ||
      strcmp(atomType, "C3"  )==0 ||
      strcmp(atomType, "CE1" )==0 ||
      strcmp(atomType, "CE2" )==0 ||
      strcmp(atomType, "CST" )==0 ||
      strcmp(atomType, "CAP" )==0 ||
      strcmp(atomType, "COA" )==0 ||
      strcmp(atomType, "CPT" )==0 ||
      strcmp(atomType, "CPH1")==0 ||
      strcmp(atomType, "CPH2")==0
      ) {
      if (numBonds == 2)
        return 5;
      else if (numBonds == 3)
        return 6;
      else
        return 1;
    } else { // other Carbon
        return 1;
    }

  //Nitrogen
  } else if (atomType[0] == 'N') {
    if (//Sp3 Nitrogen
      //strcmp(atomType, "N"   )==0 ||
      //strcmp(atomType, "NH1" )==0 ||
      //strcmp(atomType, "NH2" )==0 ||
      strcmp(atomType, "NH3" )==0 ||
      //strcmp(atomType, "NC2" )==0 ||
      //strcmp(atomType, "NY"  )==0 ||
      strcmp(atomType, "NP"  )==0
      ) {
      if (numBonds == 1)
        return 11;
      else if (numBonds == 2)
        return 12;
      else if (numBonds == 3)
        return 13;
      else
        return 11;

    } else if (//SP2 Nitrogen
      strcmp(atomType, "NY"  )==0 || //
      strcmp(atomType, "NC2" )==0 || //
      strcmp(atomType, "N"   )==0 || //
      strcmp(atomType, "NH1" )==0 || //
      strcmp(atomType, "NH2" )==0 || //
      strcmp(atomType, "NR1" )==0 ||
      strcmp(atomType, "NR2" )==0 ||
      strcmp(atomType, "NR3" )==0 ||
      strcmp(atomType, "NPH" )==0 ||
      strcmp(atomType, "NC"  )==0
      ) {
      if (numBonds == 1)
        return 14;
      else if (numBonds == 2)
        return 15;
      else if (numBonds == 3)
        return 16;
      else
        return 11;
    } else { // other Nitrogen
      return 11;
    }

  //Oxygen
  } else if (atomType[0] == 'O') {
    if (//Sp3 Oxygen
      strcmp(atomType, "OH1" )==0 ||
      strcmp(atomType, "OS"  )==0 ||
      strcmp(atomType, "OC"  )==0 || //
      strcmp(atomType, "OT"  )==0
      ) {
      if (numBonds == 1)
        return 7;
      else if (numBonds == 2)
        return 8;
      else
        return 7;
    } else if ( // Sp2 Oxygen
      strcmp(atomType, "O"   )==0 ||
      strcmp(atomType, "OB"  )==0 ||
      strcmp(atomType, "OST" )==0 ||
      strcmp(atomType, "OCA" )==0 ||
      strcmp(atomType, "OM"  )==0
      ) {
      return 9;
    } else if ( // SP1 Oxygen
      strcmp(atomType, "OC"  )==0
      ) {
      return 10;
    } else { // other Oxygen
      return 7;
    }

  //Sulfur
  } else if (atomType[0] == 'S') {
      if (numBonds == 1)
        return 17;
      else
        return 18;

  //Phosphorus
  } else if (atomType[0] == 'P') {
      if (numBonds == 3)
        return 19;
      else if (numBonds == 4)
        return 20;
      else
        return 19;
  } else { // unknown atom type
    return 5;
  }
  return 5;
} // getLCPOTypeCharmm

//input type is Charmm/Amber/other
//0 - Charmm/Xplor
//1 - Amber
//2 - Plugin
//3 - Gromacs
void Molecule::assignLCPOTypes(int inputType) {
  int *heavyBonds = new int[numAtoms];
  for (int i = 0; i < numAtoms; i++)
    heavyBonds[i] = 0;
  for (int i = 0; i < numBonds; i++ ) {
    Bond curBond = bonds[i];
    int a1 = bonds[i].atom1;
    int a2 = bonds[i].atom2;
    if (atoms[a1].mass > 2.f && atoms[a2].mass > 2.f) {
      heavyBonds[a1]++;
      heavyBonds[a2]++;
    }
  }

  lcpoParamType = new int[numAtoms];

  int warning = 0;
  for (int i = 0; i < numAtoms; i++) {
    //print vdw_type and numbonds

    if (inputType == 1) { // Amber
      lcpoParamType[i] = getLCPOTypeAmber(atomNames[i].atomtype, heavyBonds[i]);
    } else { // Charmm
      lcpoParamType[i] = getLCPOTypeCharmm(atomNames[i].atomtype, heavyBonds[i]);
    }
/*
    CkPrintf("%d MOL: ATOM[%05d] = { %4s %d } : %d\n",
      inputType,
      i+1,
      atomNames[i].atomtype,
      heavyBonds[i],
      lcpoParamType[i]
      );
*/
    if ( atoms[i].mass < 1.5 && lcpoParamType[i] != 0 ) {
      if (atoms[i].status & LonepairAtom) {
        warning |= LonepairAtom;
        lcpoParamType[i] = 0;  // reset LCPO type for LP
      }
      else if (atoms[i].status & DrudeAtom) {
        warning |= DrudeAtom;
        lcpoParamType[i] = 0;  // reset LCPO type for Drude
      }
      else {
        CkPrintf("ERROR in Molecule::assignLCPOTypes(): "
            "Light atom given heavy atom LCPO type.\n");
      }
    }

    //CkPrintf("VDW_TYPE %02d %4s\n", atoms[i].vdw_type, atomNames[i].atomtype);
  } // for atoms

  if (warning & LonepairAtom) {
    iout << iWARN << "LONE PAIRS TO BE IGNORED BY SASA\n" << endi;
  }
  if (warning & DrudeAtom) {
    iout << iWARN << "DRUDE PARTICLES TO BE IGNORED BY SASA\n" << endi;
  }

  delete [] heavyBonds;

} // buildLCPOTable

void Molecule::plgLoadAtomBasics(molfile_atom_t *atomarray){
    atoms = new Atom[numAtoms];
    atomNames = new AtomNameInfo[numAtoms];
    if(simParams->genCompressedPsf) {
        atomSegResids = new AtomSegResInfo[numAtoms];
    }    
    hydrogenGroup.resize(0);

    ResidueLookupElem *tmpResLookup = resLookup;

    for(int i=0; i<numAtoms; i++) {
        int reslength = strlen(atomarray[i].resname)+1;
        int namelength = strlen(atomarray[i].name)+1;
        int typelength = strlen(atomarray[i].type)+1;
        atomNames[i].resname = nameArena->getNewArray(reslength);
        atomNames[i].atomname = nameArena->getNewArray(namelength);
        atomNames[i].atomtype = nameArena->getNewArray(typelength);
        strcpy(atomNames[i].resname, atomarray[i].resname);
        strcpy(atomNames[i].atomname, atomarray[i].name);
        strcpy(atomNames[i].atomtype, atomarray[i].type);

        atoms[i].mass = atomarray[i].mass;
        atoms[i].charge = atomarray[i].charge;
        atoms[i].status = UnknownAtom;

        //add this atom to residue lookup table
        if(tmpResLookup) {
            tmpResLookup = tmpResLookup->append(atomarray[i].segid, atomarray[i].resid, i);
        }

        if(atomSegResids) { //for compressing molecule information
            AtomSegResInfo *one = atomSegResids + i;
            memcpy(one->segname, atomarray[i].segid, strlen(atomarray[i].segid)+1);
            one->resid = atomarray[i].resid;
        }
        //Determine the type of the atom
        if ( simParams->ignoreMass ) {
        }else if(atoms[i].mass <= 0.05) {
            atoms[i].status |= LonepairAtom;
        }else if(atoms[i].mass < 1.0) {
            atoms[i].status |= DrudeAtom;
        }else if(atoms[i].mass <= 3.5) {
            atoms[i].status |= HydrogenAtom;
        }else if((atomNames[i].atomname[0] == 'O') &&
                 (atoms[i].mass>=14.0) && (atoms[i].mass<=18.0)){
            atoms[i].status |= OxygenAtom;
        }
        //Look up the vdw constants for this atom
        params->assign_vdw_index(atomNames[i].atomtype, &atoms[i]);
    }
}

void Molecule::plgLoadBonds(int *from, int *to){
    bonds = new Bond[numBonds];
    int realNumBonds = 0;
    for(int i=0; i<numBonds; i++) {
        Bond *thisBond = bonds+realNumBonds;
        thisBond->atom1 = from[i]-1;
        thisBond->atom2 = to[i]-1;

        params->assign_bond_index(
            atomNames[thisBond->atom1].atomtype,
            atomNames[thisBond->atom2].atomtype,
            thisBond);

        //Make sure this isn't a fake bond meant for shake in x-plor
        Real k, x0;
        params->get_bond_params(&k, &x0, thisBond->bond_type);
        if (is_lonepairs_psf) {
            //need to retain Lonepair bonds for Drude
            if(k!=0. || is_lp(thisBond->atom1) || 
               is_lp(thisBond->atom2)) {               
                realNumBonds++;
            }
        }else{
            if(k != 0.) realNumBonds++;
        }
    }

    if(numBonds != realNumBonds) {
        iout << iWARN << "Ignored" << numBonds-realNumBonds <<
            "bonds with zero force constants.\n" <<endi;
        iout << iWARN << "Will get H-H distance in rigid H20 from H-O-H angle.\n" <<endi;
    }
    numBonds = realNumBonds;
}

void Molecule::plgLoadAngles(int *plgAngles)
{    
    angles=new Angle[numAngles];
    int *atomid = plgAngles;
    int numRealAngles = 0;
    for(int i=0; i<numAngles; i++) {
        Angle *thisAngle = angles+numRealAngles;
        thisAngle->atom1 = atomid[0]-1;
        thisAngle->atom2 = atomid[1]-1;
        thisAngle->atom3 = atomid[2]-1;
        atomid += 3;

        params->assign_angle_index(
            atomNames[thisAngle->atom1].atomtype,
            atomNames[thisAngle->atom2].atomtype,
            atomNames[thisAngle->atom3].atomtype,
                                thisAngle, simParams->alchOn ? -1 : 0);
        if ( thisAngle->angle_type == -1 ) {
          iout << iWARN << "ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
               << endi;
        }

        Real k, t0, k_ub, r_ub;
        if ( thisAngle->angle_type == -1 ) { k = -1.;  k_ub = -1.; } else
        params->get_angle_params(&k, &t0, &k_ub, &r_ub, thisAngle->angle_type);
        if(k!=0. || k_ub!=0.) numRealAngles++;
    }

    if(numAngles != numRealAngles) {
        iout << iWARN << "Ignored" << numAngles-numRealAngles << 
            " angles with zero force constants.\n" << endi; 
    }
    numAngles = numRealAngles;
}

void Molecule::plgLoadDihedrals(int *plgDihedrals)
{
    std::map< std::string, int > cache;

    int lastAtomIds[4];
    int multiplicity = 1; //multiplicity of the current bond

    lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
    dihedrals = new Dihedral[numDihedrals];
    int numRealDihedrals = 0;
    int *atomid = plgDihedrals;
    for(int i=0; i<numDihedrals; i++, atomid+=4) {
        Dihedral *thisDihedral = dihedrals + numRealDihedrals;
        Bool duplicate_bond = TRUE;
        for(int j=0; j<4; j++) {
            if(atomid[j] != lastAtomIds[j]) {
                duplicate_bond = FALSE;
            }
            lastAtomIds[j] = atomid[j];
        }

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

        thisDihedral->atom1 = atomid[0]-1;
        thisDihedral->atom2 = atomid[1]-1;
        thisDihedral->atom3 = atomid[2]-1;
        thisDihedral->atom4 = atomid[3]-1;

      char query[128];
      sprintf(query,"%10s :: %10s :: %10s :: %10s :: %d",
        atomNames[atomid[0]-1].atomtype,
        atomNames[atomid[1]-1].atomtype,
        atomNames[atomid[2]-1].atomtype,
        atomNames[atomid[3]-1].atomtype,
        multiplicity);
      auto search = cache.find(query);
      if ( search != cache.end() ) { 
        thisDihedral->dihedral_type = search->second;
      } else {
        params->assign_dihedral_index(
            atomNames[atomid[0]-1].atomtype,
            atomNames[atomid[1]-1].atomtype,
            atomNames[atomid[2]-1].atomtype,
            atomNames[atomid[3]-1].atomtype,
            thisDihedral, multiplicity, simParams->alchOn ? -1 : 0);
        if ( thisDihedral->dihedral_type == -1 ) {
          iout << iWARN << "ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
               << endi;
        }
        cache[query] = thisDihedral->dihedral_type;
      }
    }

    numDihedrals = numRealDihedrals;
}

void Molecule::plgLoadImpropers(int *plgImpropers)
{
    int lastAtomIds[4];
    int multiplicity = 1; //multiplicity of the current bond

    lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
    impropers = new Improper[numImpropers];
    int numRealImpropers = 0;
    int *atomid = plgImpropers;
    for(int i=0; i<numImpropers; i++, atomid+=4) {
        Improper *thisImproper = impropers + numRealImpropers;
        Bool duplicate_bond = TRUE;
        for(int j=0; j<4; j++) {
            if(atomid[j] != lastAtomIds[j]) {
                duplicate_bond = FALSE;
            }
            lastAtomIds[j] = atomid[j];
        }

        if(duplicate_bond) {
            multiplicity++;
            if(multiplicity==2) {
                numMultipleImpropers++;
            }
        }else{
            multiplicity=1;
            numRealImpropers++;
        }

        thisImproper->atom1 = atomid[0]-1;
        thisImproper->atom2 = atomid[1]-1;
        thisImproper->atom3 = atomid[2]-1;
        thisImproper->atom4 = atomid[3]-1;

        params->assign_improper_index(
            atomNames[atomid[0]-1].atomtype,
            atomNames[atomid[1]-1].atomtype,
            atomNames[atomid[2]-1].atomtype,
            atomNames[atomid[3]-1].atomtype,
            thisImproper, multiplicity);
    }

    numImpropers = numRealImpropers;
}

void Molecule::plgLoadCrossterms(int *plgCterms)
{
    int lastAtomIds[8];    

    for(int i=0; i<8; i++)
        lastAtomIds[i]=-1;
    
    crossterms = new Crossterm[numCrossterms];
    int numRealCrossterms = 0;
    int *atomid = plgCterms;
    for(int i=0; i<numCrossterms; i++, atomid+=8) {
        Crossterm *thisCrossterm = crossterms + numRealCrossterms;
        Bool duplicate_bond = TRUE;
        for(int j=0; j<8; j++) {
            if(atomid[j] != lastAtomIds[j]) {
                duplicate_bond = FALSE;
            }
            lastAtomIds[j] = atomid[j];
        }

        if(duplicate_bond) {
            iout << iWARN <<"Duplicate cross-term detected.\n" << endi;
        } else
            numRealCrossterms++;

        thisCrossterm->atom1 = atomid[0]-1;
        thisCrossterm->atom2 = atomid[1]-1;
        thisCrossterm->atom3 = atomid[2]-1;
        thisCrossterm->atom4 = atomid[3]-1;
        thisCrossterm->atom5 = atomid[4]-1;
        thisCrossterm->atom6 = atomid[5]-1;
        thisCrossterm->atom7 = atomid[6]-1;
        thisCrossterm->atom8 = atomid[7]-1;

        params->assign_crossterm_index(
            atomNames[atomid[0]-1].atomtype,
            atomNames[atomid[1]-1].atomtype,
            atomNames[atomid[2]-1].atomtype,
            atomNames[atomid[3]-1].atomtype,
            atomNames[atomid[4]-1].atomtype,
            atomNames[atomid[5]-1].atomtype,
            atomNames[atomid[6]-1].atomtype,
            atomNames[atomid[7]-1].atomtype,
            thisCrossterm);
    }

    numCrossterms = numRealCrossterms;
}

void Molecule::setOccupancyData(molfile_atom_t *atomarray){
    occupancy = new float[numAtoms];
    for(int i=0; i<numAtoms; i++) {
        occupancy[i] = atomarray[i].occupancy;
    }
}

void Molecule::setBFactorData(molfile_atom_t *atomarray){
    bfactor = new float[numAtoms];
    for(int i=0; i<numAtoms; i++) {
        bfactor[i] = atomarray[i].bfactor;
    }
}

    /************************************************************************/
    /*                  */
    /*      FUNCTION build_lists_by_atom      */
    /*                  */
    /*  This function builds O(NumAtoms) arrays that store the bonds,   */
    /*  angles, dihedrals, and impropers, that each atom is involved in.    */
    /*  This is a space hog, but VERY fast.  This will certainly have to    */
    /*  change to make things scalable in memory, but for now, speed is the */
    /*  thing!                */
    /*                  */
    /************************************************************************/
    void Molecule::build_lists_by_atom()       
    {
       register int i;      //  Loop counter

       register int numFixedAtoms = this->numFixedAtoms;
       // if we want forces on fixed atoms then just pretend
       // there are none for the purposes of this routine;
       if ( simParams->fixedAtomsForces ) numFixedAtoms = 0;

//fepb
//     int numFepInitial = this->numFepInitial;
//     int numFepFinal = this->numFepFinal;
//fepe
       tmpArena = new ObjectArena<int32>;
       bondsWithAtom = new int32 *[numAtoms];
       cluster = new int32 [numAtoms];
       clusterSize = new int32 [numAtoms];

       bondsByAtom = new int32 *[numAtoms];
       anglesByAtom = new int32 *[numAtoms];
       dihedralsByAtom = new int32 *[numAtoms];
       impropersByAtom = new int32 *[numAtoms];
       crosstermsByAtom = new int32 *[numAtoms];

       exclusionsByAtom = new int32 *[numAtoms];
       fullExclusionsByAtom = new int32 *[numAtoms];
       modExclusionsByAtom = new int32 *[numAtoms];

       // JLai
       gromacsPairByAtom = new int32 *[numAtoms];
       // End of JLai

       int32 *byAtomSize = new int32[numAtoms];

       const int pair_self = 
         simParams->pairInteractionOn ? simParams->pairInteractionSelf : 0;

       DebugM(3,"Building bond lists.\n");
    
       //  Build the bond lists
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       for (i=0; i<numRealBonds; i++)
       {
         byAtomSize[bonds[i].atom1]++;
         byAtomSize[bonds[i].atom2]++;
       }
       for (i=0; i<numAtoms; i++)
       {
         bondsWithAtom[i] = tmpArena->getNewArray(byAtomSize[i]+1);
         bondsWithAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numRealBonds; i++)
       {
         int a1 = bonds[i].atom1;
         int a2 = bonds[i].atom2;
         bondsWithAtom[a1][byAtomSize[a1]++] = i;
         bondsWithAtom[a2][byAtomSize[a2]++] = i;
       }

        
       // Updates all bond, angle, dihedral, improper and crossterm
       // to reflect the QM region (which can't have any of there terms)
       if (simParams->qmForcesOn) {
           
           DebugM(3,"Calculating exclusions for QM simulation.\n");
           build_exclusions();
           
           delete_qm_bonded() ;
           
           DebugM(3,"Re-Building bond lists.\n");
           
           // We re-calculate the bondsWithAtom list for cluster 
           // info calculation below.
           for (i=0; i<numAtoms; i++)
           {
             byAtomSize[i] = 0;
           }
           for (i=0; i<numRealBonds; i++)
           {
             byAtomSize[bonds[i].atom1]++;
             byAtomSize[bonds[i].atom2]++;
           }
           for (i=0; i<numAtoms; i++)
           {
             bondsWithAtom[i][byAtomSize[i]] = -1;
             byAtomSize[i] = 0;
           }
           for (i=0; i<numRealBonds; i++)
           {
             int a1 = bonds[i].atom1;
             int a2 = bonds[i].atom2;
             bondsWithAtom[a1][byAtomSize[a1]++] = i;
             bondsWithAtom[a2][byAtomSize[a2]++] = i;
           }
       }
        
       //  Build cluster information (contiguous clusters)
       for (i=0; i<numAtoms; i++) {
         cluster[i] = i;
       }
       for (i=0; i<numAtoms; i++) {
         int ci = i;
         while ( cluster[ci] != ci ) ci = cluster[ci];
         for ( int32 *b = bondsWithAtom[i]; *b != -1; ++b ) {
           int a = bonds[*b].atom1;
           if ( a == i ) a = bonds[*b].atom2;
           if ( a > i ) {
             int ca = a;
             while ( cluster[ca] != ca ) ca = cluster[ca];
             if ( ca > ci ) cluster[ca] = cluster[ci];
             else cluster[ci] = cluster[ca];
           }
         }
       }
       while ( 1 ) {
         int allok = 1;
         for (i=0; i<numAtoms; i++) {
           int ci = cluster[i];
           if ( cluster[ci] != ci ) {
             allok = 0;
             cluster[i] = cluster[ci];
           }
         }
         if ( allok ) break;
       }
       
       for (i=0; i<numAtoms; i++) {
         clusterSize[i] = 0;
       }       
       for (i=0; i<numAtoms; i++) {           
         clusterSize[cluster[i]] += 1;
       }

/*
       //Getting number of clusters for debugging
       int numClusters=0;
       for(int i=0; i<numAtoms; i++){
           if(clusterSize[i]!=0) numClusters++;
       }
       printf("Num of clusters: %d\n", numClusters);
*/

       //  Build the bond lists
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       numCalcBonds = 0;
       for (i=0; i<numBonds; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
                            && fixedAtomFlags[bonds[i].atom2] ) continue;
   
         if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
         byAtomSize[bonds[i].atom1]++;
         numCalcBonds++;
       }
       for (i=0; i<numAtoms; i++)
       {
         bondsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         bondsByAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numBonds; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
                            && fixedAtomFlags[bonds[i].atom2] ) continue;
         if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
         int a1 = bonds[i].atom1;
         bondsByAtom[a1][byAtomSize[a1]++] = i;
       }
       for (i=0; i<numBonds; ++i) {
         int a1 = bonds[i].atom1;
         int a2 = bonds[i].atom2;
         int j;
         if ( a1 == a2 ) {
           char buff[512];
           sprintf(buff,"Atom %d is bonded to itself", a1+1);
           NAMD_die(buff);
         }
         for ( j = 0; j < byAtomSize[a1]; ++j ) {
           int b = bondsByAtom[a1][j];
           int ba1 = bonds[b].atom1;
           int ba2 = bonds[b].atom2;
           if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
             char buff[512];
             sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
             NAMD_die(buff);
           }
         }
         for ( j = 0; j < byAtomSize[a2]; ++j ) {
           int b = bondsByAtom[a2][j];
           int ba1 = bonds[b].atom1;
           int ba2 = bonds[b].atom2;
           if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
             char buff[512];
             sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
             NAMD_die(buff);
           }
         }
       }
       
       DebugM(3,"Building angle lists.\n");
    
       //  Build the angle lists
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       numCalcAngles = 0;
       for (i=0; i<numAngles; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
                            && fixedAtomFlags[angles[i].atom2]
                            && fixedAtomFlags[angles[i].atom3] ) continue;
         if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
         byAtomSize[angles[i].atom1]++;
         numCalcAngles++;
       }
       for (i=0; i<numAtoms; i++)
       {
         anglesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         anglesByAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numAngles; i++)
       {
         if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
         if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
                            && fixedAtomFlags[angles[i].atom2]
                            && fixedAtomFlags[angles[i].atom3] ) continue;
         int a1 = angles[i].atom1;
         anglesByAtom[a1][byAtomSize[a1]++] = i;
       }
       
       DebugM(3,"Building improper lists.\n");
    
       //  Build the improper lists
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       numCalcImpropers = 0;
       for (i=0; i<numImpropers; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
                            && fixedAtomFlags[impropers[i].atom2]
                            && fixedAtomFlags[impropers[i].atom3]
                            && fixedAtomFlags[impropers[i].atom4] ) continue;
         if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
         byAtomSize[impropers[i].atom1]++;
         numCalcImpropers++;
       }
       for (i=0; i<numAtoms; i++)
       {
         impropersByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         impropersByAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numImpropers; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
                            && fixedAtomFlags[impropers[i].atom2]
                            && fixedAtomFlags[impropers[i].atom3]
                            && fixedAtomFlags[impropers[i].atom4] ) continue;
         if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
         int a1 = impropers[i].atom1;
         impropersByAtom[a1][byAtomSize[a1]++] = i;
       }
       
       DebugM(3,"Building dihedral lists.\n");
    
       //  Build the dihedral lists
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       numCalcDihedrals = 0;
       for (i=0; i<numDihedrals; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
                            && fixedAtomFlags[dihedrals[i].atom2]
                            && fixedAtomFlags[dihedrals[i].atom3]
                            && fixedAtomFlags[dihedrals[i].atom4] ) continue;
         if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
         byAtomSize[dihedrals[i].atom1]++;
         numCalcDihedrals++;
       }
       for (i=0; i<numAtoms; i++)
       {
         dihedralsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         dihedralsByAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numDihedrals; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
                            && fixedAtomFlags[dihedrals[i].atom2]
                            && fixedAtomFlags[dihedrals[i].atom3]
                            && fixedAtomFlags[dihedrals[i].atom4] ) continue;
         if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
         int a1 = dihedrals[i].atom1;
         dihedralsByAtom[a1][byAtomSize[a1]++] = i;
       }
    
       DebugM(3,"Building crossterm lists.\n");
    
       //  Build the crossterm lists
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       numCalcCrossterms = 0;
       for (i=0; i<numCrossterms; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
                            && fixedAtomFlags[crossterms[i].atom2]
                            && fixedAtomFlags[crossterms[i].atom3]
                            && fixedAtomFlags[crossterms[i].atom4]
                            && fixedAtomFlags[crossterms[i].atom5]
                            && fixedAtomFlags[crossterms[i].atom6]
                            && fixedAtomFlags[crossterms[i].atom7]
                            && fixedAtomFlags[crossterms[i].atom8] ) continue;
         if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
         byAtomSize[crossterms[i].atom1]++;
         numCalcCrossterms++;
       }
       for (i=0; i<numAtoms; i++)
       {
         crosstermsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         crosstermsByAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numCrossterms; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
                            && fixedAtomFlags[crossterms[i].atom2]
                            && fixedAtomFlags[crossterms[i].atom3]
                            && fixedAtomFlags[crossterms[i].atom4]
                            && fixedAtomFlags[crossterms[i].atom5]
                            && fixedAtomFlags[crossterms[i].atom6]
                            && fixedAtomFlags[crossterms[i].atom7]
                            && fixedAtomFlags[crossterms[i].atom8] ) continue;
         if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
         int a1 = crossterms[i].atom1;
         crosstermsByAtom[a1][byAtomSize[a1]++] = i;
       }

       // DRUDE: init lphostIndexes array
       if (is_lonepairs_psf) {
         // allocate lone pair host index array only if we need it!
         DebugM(3,"Initializing lone pair host index array.\n");
         lphostIndexes = new int32[numAtoms];
         for (i = 0;  i < numAtoms;  i++) {
           lphostIndexes[i] = -1;
         }
         for (i = 0;  i < numLphosts;  i++) {
           int32 index = lphosts[i].atom1;
           lphostIndexes[index] = i;
         }
       }
       // DRUDE
    
       // JLai
       DebugM(3,"Building gromacsPair lists.\n");
    
       //  Build the gromacsPair lists
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       numCalcLJPair = 0;
       for (i=0; i<numLJPair; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[gromacsPair[i].atom1]
                            && fixedAtomFlags[gromacsPair[i].atom2] ) continue;
         if ( pair_self && fepAtomFlags[gromacsPair[i].atom1] != 1) continue;
         byAtomSize[gromacsPair[i].atom1]++;
         numCalcLJPair++;
       }
       for (i=0; i<numAtoms; i++)
       {
         gromacsPairByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         gromacsPairByAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numLJPair; i++)
       {
           if ( numFixedAtoms && fixedAtomFlags[gromacsPair[i].atom1]
                && fixedAtomFlags[gromacsPair[i].atom2] ) continue;
           if ( pair_self && fepAtomFlags[gromacsPair[i].atom1] != 1) continue;
           int a1 = gromacsPair[i].atom1;
           gromacsPairByAtom[a1][byAtomSize[a1]++] = i;
           }

       // End of JLai

       DebugM(3,"Building exclusion data.\n");
    
       //  Build the arrays of exclusions for each atom
       if (! simParams->qmForcesOn)
       build_exclusions();

       //  Remove temporary structures
       delete [] bondsWithAtom;  bondsWithAtom = 0;
       delete tmpArena;  tmpArena = 0;

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

       // 1-4 exclusions which are also fully excluded were eliminated by hash table
       numTotalExclusions = exclusionSet.size();
       if ( ! CkMyPe() ) {
         iout << iINFO << "ADDED " << (numTotalExclusions - numExclusions) << " IMPLICIT EXCLUSIONS\n" << endi;
       }
       exclusions = new Exclusion[numTotalExclusions];
       UniqueSetIter<Exclusion> exclIter(exclusionSet);
       for ( exclIter=exclIter.begin(),i=0; exclIter != exclIter.end(); exclIter++,i++ )
       {
         exclusions[i] = *exclIter;
       }
       // Free exclusionSet storage
       // exclusionSet.clear(1);
       exclusionSet.clear();

       DebugM(3,"Building exclusion lists.\n");
    
       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
       }
       numCalcExclusions = 0;
       numCalcFullExclusions = 0;
       for (i=0; i<numTotalExclusions; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
                            && fixedAtomFlags[exclusions[i].atom2] ) continue;
         byAtomSize[exclusions[i].atom1]++;
         numCalcExclusions++;
         if ( ! exclusions[i].modified ) numCalcFullExclusions++;
       }

       for (i=0; i<numAtoms; i++)
       {
         exclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         exclusionsByAtom[i][byAtomSize[i]] = -1;
         byAtomSize[i] = 0;
       }
       for (i=0; i<numTotalExclusions; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
                            && fixedAtomFlags[exclusions[i].atom2] ) continue;
         int a1 = exclusions[i].atom1;
         exclusionsByAtom[a1][byAtomSize[a1]++] = i;
       }

       int32 *byAtomSize2 = new int32[numAtoms];

       for (i=0; i<numAtoms; i++)
       {
         byAtomSize[i] = 0;
         byAtomSize2[i] = 0;
       }

       for (i=0; i<numTotalExclusions; i++)
       {
         if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
                            && fixedAtomFlags[exclusions[i].atom2] ) continue;
         if ( exclusions[i].modified ) {
           byAtomSize2[exclusions[i].atom1]++;
           byAtomSize2[exclusions[i].atom2]++;
         } else {
           byAtomSize[exclusions[i].atom1]++;
           byAtomSize[exclusions[i].atom2]++;
         }
       }

       for (i=0; i<numAtoms; i++)
       {
         fullExclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
         fullExclusionsByAtom[i][0] = 0;
         modExclusionsByAtom[i] = arena->getNewArray(byAtomSize2[i]+1);
         modExclusionsByAtom[i][0] = 0;
       }

       for (i=0; i<numTotalExclusions; i++)
       {
         int a1 = exclusions[i].atom1;
         int a2 = exclusions[i].atom2;
         if ( numFixedAtoms && fixedAtomFlags[a1]
                            && fixedAtomFlags[a2] ) continue;
         int32 *l1, *l2;
         if ( exclusions[i].modified ) {
           l1 = modExclusionsByAtom[a1];
           l2 = modExclusionsByAtom[a2];
         } else {
           l1 = fullExclusionsByAtom[a1];
           l2 = fullExclusionsByAtom[a2];
         }
         l1[++(*l1)] = a2;
         l2[++(*l2)] = a1;
       }

       if ( ! CkMyPe() && simParams->printExclusions ) {
         for (i=0; i<numAtoms; i++) {
           int32 *lf = fullExclusionsByAtom[i];
           iout << "EXCL " << i << " FULL";
           int nf = *(lf++);
           for ( int j = 0; j < nf; ++j ) {
             iout << " " << *(lf++);
           }
           iout << "\n";
           int32 *lm = modExclusionsByAtom[i];
           iout << "EXCL " << i << " MOD";
           int nm = *(lm++);
           for ( int j = 0; j < nm; ++j ) {
             iout << " " << *(lm++);
           }
           iout << "\n" << endi;
         }
       }

       // DRUDE
       if (is_drude_psf || simParams->drudeOn) {

         // build Thole (screened Coulomb) correction terms;
         // they are constructed implicitly from exclusions

         // free the previous Thole array if already allocated
         if (tholes != NULL) delete[] tholes;
         numTholes = 0;

         // count the number of Thole terms
         for (i = 0;  i < numTotalExclusions;  i++) {
           /* skip over the modified exclusions */
           if (exclusions[i].modified) continue;
           int a1 = exclusions[i].atom1;
           int a2 = exclusions[i].atom2;
           if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
             numTholes++;
           }
         }

         // allocate space for Thole terms
         if (numTholes != 0) tholes = new Thole[numTholes];
         else tholes = NULL;
         int nt = 0;

         Real c = COULOMB*simParams->nonbondedScaling/simParams->dielectric;

         // store Thole terms
         for (i = 0;  i < numTotalExclusions;  i++) {
           /* skip over the modified exclusions */
           if (exclusions[i].modified) continue;
           int a1 = exclusions[i].atom1;
           int a2 = exclusions[i].atom2;
           // exclusions are stored with a1 < a2
           if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
             Real thsum = drudeConsts[a1].thole + drudeConsts[a2].thole;
             Real aprod = drudeConsts[a1].alpha * drudeConsts[a2].alpha;
             // guard against having alpha==0
             Real apower = (aprod <= 0 ? 0 : powf(aprod, -1.f/6));
             tholes[nt].atom1 = a1;
             tholes[nt].atom2 = a1+1;
             tholes[nt].atom3 = a2;
             tholes[nt].atom4 = a2+1;
             tholes[nt].aa = apower * thsum;
             tholes[nt].qq = c * atoms[a1+1].charge * atoms[a2+1].charge;
             nt++;
           }
         }

         // build Thole lists by atom
         DebugM(3, "Building Thole correction term lists.\n");
         tholesByAtom = new int32 *[numAtoms];

         for (i = 0;  i < numAtoms;  i++) {
           byAtomSize[i] = 0;
         }
         numCalcTholes = 0;
         for (i = 0;  i < numTholes;  i++) {
           if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
                              && fixedAtomFlags[tholes[i].atom2]
                              && fixedAtomFlags[tholes[i].atom3]
                              && fixedAtomFlags[tholes[i].atom4] ) continue;
           if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1) continue;
           byAtomSize[tholes[i].atom1]++;
           numCalcTholes++;
         }
         for (i = 0;  i < numAtoms;  i++) {
           tholesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
           tholesByAtom[i][byAtomSize[i]] = -1;
           byAtomSize[i] = 0;
         }
         for (i = 0;  i < numTholes;  i++) {
           if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
                              && fixedAtomFlags[tholes[i].atom2]
                              && fixedAtomFlags[tholes[i].atom3]
                              && fixedAtomFlags[tholes[i].atom4] ) continue;
           if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1) continue;
           int a1 = tholes[i].atom1;
           tholesByAtom[a1][byAtomSize[a1]++] = i;
         }

         // build anisotropic lists by atom
         DebugM(3, "Building anisotropic term lists.\n");
         anisosByAtom = new int32 *[numAtoms];

         for (i = 0;  i < numAtoms;  i++) {
           byAtomSize[i] = 0;
         }
         numCalcAnisos = 0;
         for (i = 0;  i < numAnisos;  i++) {
           if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
                              && fixedAtomFlags[anisos[i].atom2]
                              && fixedAtomFlags[anisos[i].atom3]
                              && fixedAtomFlags[anisos[i].atom4] ) continue;
           if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1) continue;
           byAtomSize[anisos[i].atom1]++;
           numCalcAnisos++;
         }
         for (i = 0;  i < numAtoms;  i++) {
           anisosByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
           anisosByAtom[i][byAtomSize[i]] = -1;
           byAtomSize[i] = 0;
         }
         for (i = 0;  i < numAnisos;  i++) {
           if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
                              && fixedAtomFlags[anisos[i].atom2]
                              && fixedAtomFlags[anisos[i].atom3]
                              && fixedAtomFlags[anisos[i].atom4] ) continue;
           if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1) continue;
           int a1 = anisos[i].atom1;
           anisosByAtom[a1][byAtomSize[a1]++] = i;
         }

       }
       // DRUDE

       delete [] byAtomSize;  byAtomSize = 0;
       delete [] byAtomSize2;  byAtomSize2 = 0;


       //  Allocate an array to hold the exclusions for each atom
       all_exclusions = new ExclusionCheck[numAtoms];

       for (i=0; i<numAtoms; i++)
       {
         all_exclusions[i].min = numAtoms;
         all_exclusions[i].max = -1;
       }
       for (i=0; i<numTotalExclusions; i++)
       {
         // first atom should alway have lower number!
         int a1 = exclusions[i].atom1;
         int a2 = exclusions[i].atom2;
         if ( numFixedAtoms && fixedAtomFlags[a1]
                            && fixedAtomFlags[a2] ) continue;
         if ( all_exclusions[a1].min > a2 ) all_exclusions[a1].min = a2;
         if ( all_exclusions[a2].min > a1 ) all_exclusions[a2].min = a1;
         if ( a2 > all_exclusions[a1].max ) all_exclusions[a1].max = a2;
         if ( a1 > all_exclusions[a2].max ) all_exclusions[a2].max = a1;
       }

       // build array of all full exclusions for water etc.
       int maxDenseAllFull = 0;
       int numDenseAllFull = 0;
       for (i=0; i<numAtoms; i++) {
         int iInMiddle = ( i < all_exclusions[i].max &&
                       i > all_exclusions[i].min ) ? 1 : 0;
         int s = all_exclusions[i].max - all_exclusions[i].min + 1;
         if ( s == fullExclusionsByAtom[i][0] + iInMiddle ) {
            if ( s > maxDenseAllFull ) maxDenseAllFull = s;
            all_exclusions[i].flags = (char*)-1;  // shared array
         } else {
            all_exclusions[i].flags = 0;  // individual array
         }
       }
       char *denseFullArray = exclArena->getNewArray(maxDenseAllFull);
       for ( i=0; i<maxDenseAllFull; ++i ) denseFullArray[i] = EXCHCK_FULL;

       int exclmem = maxDenseAllFull;
       int maxExclusionFlags = simParams->maxExclusionFlags;
       for (i=0; i<numAtoms; i++) {
         int s = all_exclusions[i].max - all_exclusions[i].min + 1;
         if ( all_exclusions[i].max != -1 ) {
           if ( all_exclusions[i].flags ) {
             all_exclusions[i].flags = denseFullArray;
             ++numDenseAllFull;
           } else if ( s < maxExclusionFlags ) {
             char *f = all_exclusions[i].flags = exclArena->getNewArray(s);
             for ( int k=0; k<s; ++k ) f[k] = 0;
             exclmem += s;
           } else {
             all_exclusions[i].flags = 0;  // need to build on the fly
           }
         } else {
           all_exclusions[i].flags = (char*)-1; // should never dereference
         }
       }
       if ( 0 ) {
         iout << iINFO << numTotalExclusions << " exclusions consume "
            << exclmem << " bytes.\n" << endi;
         iout << iINFO << numDenseAllFull
            << " atoms sharing one array.\n" << endi;
       }
       for (i=0; i<numTotalExclusions; i++)
       {
         int a1 = exclusions[i].atom1;
         int a2 = exclusions[i].atom2;
         if ( numFixedAtoms && fixedAtomFlags[a1]
                            && fixedAtomFlags[a2] ) continue;
         if ( exclusions[i].modified ) {
           if ( all_exclusions[a1].flags )
             all_exclusions[a1].flags[a2-all_exclusions[a1].min] = EXCHCK_MOD;
           if ( all_exclusions[a2].flags )
             all_exclusions[a2].flags[a1-all_exclusions[a2].min] = EXCHCK_MOD;
         } else {
           if ( all_exclusions[a1].flags )
             all_exclusions[a1].flags[a2-all_exclusions[a1].min] = EXCHCK_FULL;
           if ( all_exclusions[a2].flags )
             all_exclusions[a2].flags[a1-all_exclusions[a2].min] = EXCHCK_FULL;
         }
       }
    }

    /*    END OF FUNCTION build_lists_by_atom    */

    /****************************************************************/
    /*                */
    /*      FUNCTION build_exclusions    */
    /*                */
    /*  This function builds a list of all the exlcusions       */
    /*  atoms.  These lists include explicit exclusions as well as  */
    /*  exclusions that are calculated based on the bonded structure*/
    /*  and the exclusion flag.  For each pair of atoms that are    */
    /*  excluded, the larger of the 2 atom indexes is stored in the */
    /*  array of the smaller index.  All the arrays are not sorted. */
    /*  Then to determine if two atoms have an exclusion, a linear  */
    /*  search is done on the array of the atom with the smaller    */
    /*  index for the larger index.          */
    /*  If the exclusion policy is set to scaled1-4, there are  */
    /*  actually two lists built.  One contains the pairs of atoms  */
    /*  that are to be exlcuded (i.e., explicit exclusions, 1-2,    */
    /*  and 1-3 interactions) and the other contains just the 1-4   */
    /*  interactions, since they will need to be determined   */
    /*  independantly of the other exclusions.      */
    /*                */
    /****************************************************************/

    void Molecule::build_exclusions()
    {
      register int i;          //  Loop counter
      ExclusionSettings exclude_flag;    //  Exclusion policy

      exclude_flag = simParams->exclude;

      //  Go through the explicit exclusions and add them to the arrays
      for (i=0; i<numExclusions; i++)
      {
        exclusionSet.add(exclusions[i]);
      }

      // If this is AMBER force field, and readExclusions is TRUE,
      // then all the exclusions were read from parm file, and we
      // shouldn't generate any of them.
      if (!simParams->amberOn || !simParams->readExclusions)
      { //  Now calculate the bonded exlcusions based on the exclusion policy
        switch (exclude_flag)
        {
         case NONE:
           break;
         case ONETWO:
           build12excl();
           break;
          case ONETHREE:
            build12excl();
            build13excl();
            break;
          case ONEFOUR:
            build12excl();
            build13excl();
            build14excl(0);
            break;
          case SCALED14:
            build12excl();
            build13excl();
            build14excl(1);
            break;
        }
      }

      stripFepExcl();

      // DRUDE
      if (is_lonepairs_psf || is_drude_psf) {
        build_inherited_excl(SCALED14 == exclude_flag);
      }
    }
    /*      END OF FUNCTION build_exclusions    */


    // Extend exclusions for the Drude model.  The Drude model is generally
    // used with the 1-3 exclusion policy, although the code below also
    // supports the 1-2 exclusion policy.  The use of light (or massless)
    // pseudo-atoms requires the introduction of extra exclusions.
    //
    // Here is the algorithm for determining Drude model exclusions:
    // (1)  Each Drude particle and each lone pair has a single parent atom.
    //      The parent atom must be a heavy atom.
    // (2)  Each Drude particle and lone pair inherit the exclusions of its
    //      parent atom.
    // (3)  If two heavy atoms are excluded and they both have either a
    //      Drude particle or a lone pair, the these light (or massless)
    //      particles are also excluded from interacting with each other.
    void Molecule::build_inherited_excl(int modified) {
      ExclusionSettings exclude_flag = simParams->exclude;
      int32 *bond1, *bond2, *bond3, *bond4, *bond5;
      int32 i, j, mid1, mid2, mid3, mid4;

      // validate that each Drude or lone pair particle
      // has a unique parent that is a heavy atom
      for (i = 0;  i < numAtoms;  i++) {

        if (!is_drude(i) && !is_lp(i)) continue;
        // make sure that i is either Drude or LP

        // find parent (heavy) atom of particle i
        bond1 = bondsWithAtom[i];

        if (-1 == *bond1) {  // i must have one bond
          char err_msg[512];
          const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
          sprintf(err_msg, "FOUND ISOLATED %s PARTICLE %d", idescrip, i+1);
          NAMD_die(err_msg);
        }
        if (-1 != *(bond1+1)) {  // and only one bond
          char err_msg[512];
          const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
          sprintf(err_msg, "FOUND MULTIPLY LINKED %s PARTICLE %d",
              idescrip, i+1);
          NAMD_die(err_msg);
        }

        // mid1 is parent of particle i
        mid1 = bonds[*bond1].atom1;
        if (mid1 == i) mid1 = bonds[*bond1].atom2;

        // make sure that mid1 is a heavy atom
        if (is_drude(mid1) || is_lp(mid1) || is_hydrogen(mid1)) {
          char err_msg[512];
          const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
          sprintf(err_msg, "PARENT ATOM %d of %s PARTICLE %d "
              "IS NOT HEAVY ATOM", mid1+1, idescrip, i+1);
          NAMD_die(err_msg);
        }

        if (exclude_flag == NONE) {
          // add (i,mid1) as an exclusion
          if (i < mid1) {
            exclusionSet.add(Exclusion(i, mid1));
          }
          else {
            exclusionSet.add(Exclusion(mid1, i));
          }

          // also exclude any Drude particles or LPs bonded to mid1
          bond2 = bondsWithAtom[mid1];
          while (*bond2 != -1) {
            j = bonds[*bond2].atom1;
            if ((is_drude(j) || is_lp(j)) && j != mid1) {
              if      (i < j) exclusionSet.add(Exclusion(i, j));
              else if (j < i) exclusionSet.add(Exclusion(j, i));
            }
            j = bonds[*bond2].atom2;
            if ((is_drude(j) || is_lp(j)) && j != mid1) {
              if      (i < j) exclusionSet.add(Exclusion(i, j));
              else if (j < i) exclusionSet.add(Exclusion(j, i));
            }
            bond2++;
          }
        }
        else {  // if ONETWO or ONETHREE or ONEFOUR or SCALED14

          // find the next link
          bond2 = bondsWithAtom[mid1];

          // loop through all the bonds connected to atom mid1
          while (*bond2 != -1) {
            if (bonds[*bond2].atom1 == mid1) {
              mid2 = bonds[*bond2].atom2;
            }
            else {
              mid2 = bonds[*bond2].atom1;
            }

            // Make sure that we don't double back to where we started from.
            // Doing so causes strange behavior.
            if (mid2 == i) {
              bond2++;
              continue;
            }

            if (exclude_flag == ONETWO) {
              // add (i,mid2) as an exclusion
              if (i < mid2) {
                exclusionSet.add(Exclusion(i, mid2));
              }
              else {
                exclusionSet.add(Exclusion(mid2, i));
              }

              // also exclude any Drude particles or LPs bonded to mid2
              bond3 = bondsWithAtom[mid2];
              while (*bond3 != -1) {
                j = bonds[*bond3].atom1;
                if ((is_drude(j) || is_lp(j)) && j != mid2) {
                  if      (i < j) exclusionSet.add(Exclusion(i, j));
                  else if (j < i) exclusionSet.add(Exclusion(j, i));
                }
                j = bonds[*bond3].atom2;
                if ((is_drude(j) || is_lp(j)) && j != mid2) {
                  if      (i < j) exclusionSet.add(Exclusion(i, j));
                  else if (j < i) exclusionSet.add(Exclusion(j, i));
                }
                bond3++;
              }
            }
            else { // if ONETHREE or ONEFOUR or SCALED14

              // find the next link
              bond3 = bondsWithAtom[mid2];

              // loop through all the bonds connected to mid2
              while (*bond3 != -1) {

                if (bonds[*bond3].atom1 == mid2) {
                  mid3 = bonds[*bond3].atom2;
                }
                else {
                  mid3 = bonds[*bond3].atom1;
                }

                // Make sure we don't double back to where we started.
                // Doing so causes strange behavior.
                if (mid3 == mid1) {
                  bond3++;
                  continue;
                }

                // add (i,mid3) as an exclusion
                if (i < mid3) {
                  exclusionSet.add(Exclusion(i, mid3));
                }
                else if (mid3 < i) {
                  exclusionSet.add(Exclusion(mid3, i));
                }

                if (exclude_flag == ONETHREE) {
                  // also exclude any Drude particles or LPs bonded to mid3
                  bond4 = bondsWithAtom[mid3];
                  while (*bond4 != -1) {
                    j = bonds[*bond4].atom1;
                    if ((is_drude(j) || is_lp(j)) && j != mid3) {
                      if      (i < j) exclusionSet.add(Exclusion(i, j));
                      else if (j < i) exclusionSet.add(Exclusion(j, i));
                    }
                    j = bonds[*bond4].atom2;
                    if ((is_drude(j) || is_lp(j)) && j != mid3) {
                      if      (i < j) exclusionSet.add(Exclusion(i, j));
                      else if (j < i) exclusionSet.add(Exclusion(j, i));
                    }
                    bond4++;
                  }
                }
                else { // if ONEFOUR or SCALED14

                  // find next link
                  bond4 = bondsWithAtom[mid3];

                  // loop through all the bonds connected to mid3
                  while (*bond4 != -1) {

                    if (bonds[*bond4].atom1 == mid3) {
                      mid4 = bonds[*bond4].atom2;
                    }
                    else {
                      mid4 = bonds[*bond4].atom1;
                    }

                    // Make sure we don't double back to where we started.
                    // Doing so causes strange behavior.
                    if (mid4 == mid2) {
                      bond4++;
                      continue;
                    }

                    if (is_drude(mid4) || is_lp(mid4)) {
                      // (i,mid4) is 1-3 excl
                      if (i < mid4) {
                        exclusionSet.add(Exclusion(i, mid4));
                      }
                      else if (mid4 < i) {
                        exclusionSet.add(Exclusion(mid4, i));
                      }
                      bond4++;
                      continue;
                    }

                    // (mid1,mid4) is an existing heavy atom exclusion
                    // if we have modified 1-4 exclusions, make sure
                    // that (mid1,mid4) is modified 1-4 exclusion
                    // rather than something closer due to a ring
                    int modi = modified;
                    if (modified) {
                      int amin = (mid1 < mid4 ? mid1 : mid4);
                      int amax = (mid1 >= mid4 ? mid1 : mid4);
                      Exclusion *pe = exclusionSet.find(Exclusion(amin,amax));
                      if (pe==0) {
                        // since there is not an existing exclusion
                        // between (mid1,mid4), don't inherit!
                        bond4++;
                        continue;
                      }
                      modi = pe->modified;
                    }

                    if (i < mid4) {
                      exclusionSet.add(Exclusion(i, mid4, modi));
                    }
                    else if (mid4 < i) {
                      exclusionSet.add(Exclusion(mid4, i, modi));                    
                    }

                    // also exclude any Drude particles or LPs bonded to mid4
                    // using the "modi" setting of (mid1,mid4) exclusion
                    bond5 = bondsWithAtom[mid4];
                    while (*bond5 != -1) {
                      j = bonds[*bond5].atom1;
                      if ((is_drude(j) || is_lp(j)) && j != mid4) {
                        if      (i<j) exclusionSet.add(Exclusion(i,j,modi));
                        else if (j<i) exclusionSet.add(Exclusion(j,i,modi));
                      }
                      j = bonds[*bond5].atom2;
                      if ((is_drude(j) || is_lp(j)) && j != mid4) {
                        if      (i<j) exclusionSet.add(Exclusion(i,j,modi));
                        else if (j<i) exclusionSet.add(Exclusion(j,i,modi));
                      }
                      bond5++;
                    }
                    ++bond4;
                  } // while bond4

                } // else (if ONEFOUR or SCALED14)

                ++bond3;
              } // while bond3

            } // else (if ONETHREE or ONEFOUR or SCALED14)
           
            ++bond2;
          } // while bond2

        } // else (if ONETWO or ONETHREE or ONEFOUR or SCALED14)

      } // for i
    } 
    // DRUDE


    /************************************************************************/
    /*                  */
    /*      FUNCTION build12excl        */
    /*                  */
    /************************************************************************/

    void Molecule::build12excl(void)       
    {
       int32 *current_val;  //  Current value to check
       register int i;    //  Loop counter to loop through all atoms
       
       //  Loop through all the atoms marking the bonded interactions for each one
       for (i=0; i<numAtoms; i++)
       {
      current_val = bondsWithAtom[i];
       
      //  Loop through all the bonds for this atom
      while (*current_val != -1)
      {
         if (bonds[*current_val].atom1 == i)
         {
      if (i<bonds[*current_val].atom2)
      {
         exclusionSet.add(Exclusion(i,bonds[*current_val].atom2));
      }
         }
         else
         {
      if (i<bonds[*current_val].atom1)
      {
         exclusionSet.add(Exclusion(i,bonds[*current_val].atom1));
      }
         }
    
         ++current_val;
      }
       }
    }
    /*      END OF FUNCTION build12excl      */

    /************************************************************************/
    /*                  */
    /*      FUNCTION build13excl        */
    /*                  */
    /************************************************************************/

    void Molecule::build13excl(void)       
    {
       int32 *bond1, *bond2;  //  The two bonds being checked
       int middle_atom;  //  Common third atom
       register int i;    //  Loop counter to loop through all atoms
       
       //  Loop through all the atoms looking at the bonded connections
       //  for each one
       for (i=0; i<numAtoms; i++)
       {
       bond1 = bondsWithAtom[i];
       
       //  Loop through all the bonds directly connect to atom i
       while (*bond1 != -1)
       {
        if (bonds[*bond1].atom1 == i)
        {
          middle_atom=bonds[*bond1].atom2;
        }
        else
        {
          middle_atom=bonds[*bond1].atom1;
        }

        bond2 = bondsWithAtom[middle_atom];

        //  Now loop through all the bonds connect to the
        //  middle atom
        while (*bond2 != -1)
        {
          if (bonds[*bond2].atom1 == middle_atom)
          {
            if (i < bonds[*bond2].atom2)
            {
              exclusionSet.add(Exclusion(i,bonds[*bond2].atom2));
            }
          }
          else
          {
            if (i < bonds[*bond2].atom1)
            {
              exclusionSet.add(Exclusion(i,bonds[*bond2].atom1));
            }
          }

          ++bond2;
        }

        ++bond1;
      }
       }
    }
    /*      END OF FUNCTION build13excl      */

    /************************************************************************/
    /*                  */
    /*        FUNCTION build14excl      */
    /*                  */
    /************************************************************************/


    void Molecule::build14excl(int modified)       
    {
       int32 *bond1, *bond2, *bond3;  //  The two bonds being checked
       int mid1, mid2;    //  Middle atoms
       register int i;      //  Counter to loop through all atoms
       
       //  Loop through all the atoms
       for (i=0; i<numAtoms; i++)
       {  
         if (is_drude(i) || is_lp(i)) continue;  // skip Drude and LP for now

      // Get all the bonds connect directly to atom i
      bond1 = bondsWithAtom[i];
       
      while (*bond1 != -1)
      {
        if (bonds[*bond1].atom1 == i)
        {
          mid1=bonds[*bond1].atom2;
        }
        else
        {
          mid1=bonds[*bond1].atom1;
        }

        bond2 = bondsWithAtom[mid1];

        //  Loop through all the bonds connected to atom mid1
        while (*bond2 != -1)
        {
          if (bonds[*bond2].atom1 == mid1)
          {
            mid2 = bonds[*bond2].atom2;
          }
          else
          {
            mid2 = bonds[*bond2].atom1;
          }

          //  Make sure that we don't double back to where
          //  we started from.  This causes strange behavior.
          //  Trust me, I've been there . . .
          if (mid2 == i)
          {
            ++bond2;
            continue;
          }

          bond3=bondsWithAtom[mid2];

          //  Loop through all the bonds connected to mid2
          while (*bond3 != -1)
          {
            if (bonds[*bond3].atom1 == mid2)
            {
              int j = bonds[*bond3].atom2;
              //  Make sure that we don't double back to where
              //  we started from.  This causes strange behavior.
              //  Trust me, I've been there . . .
              //  I added this!!!  Why wasn't it there before?  -JCP
              if (j != mid1)
              if (i < j && !is_drude(j) && !is_lp(j))  // skip Drude and LP
              {
                 exclusionSet.add(Exclusion(i,j,modified));
              }
            }
            else
            {
              int j = bonds[*bond3].atom1;
              //  Make sure that we don't double back to where
              //  we started from.  This causes strange behavior.
              //  Trust me, I've been there . . .
              //  I added this!!!  Why wasn't it there before?  -JCP
              if (j != mid1)
              if (i < j && !is_drude(j) && !is_lp(j))  // skip Drude and LP
              {
                 exclusionSet.add(Exclusion(i,j,modified));
              }
            }

            ++bond3;
          }

          ++bond2;
        }
    
        ++bond1;
      }
       }
    }
    /*      END OF FUNCTION build14excl      */


    /************************************************************************/
    /*                                                                      */
    /*        FUNCTION stripFepExcl                                         */
    /*                                                                      */
    /************************************************************************/
  void Molecule::stripFepExcl(void)
  {   
    UniqueSet<Exclusion> fepExclusionSet;
    UniqueSetIter<Exclusion> exclIter(exclusionSet);

    if ( simParams->alchOn || simParams->lesOn ) {
       for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
       {
         int t1 = get_fep_type(exclIter->atom1);
         int t2 = get_fep_type(exclIter->atom2);
         if (t1 && t2 && t1 !=t2 && abs(t1-t2) != 2) {
           fepExclusionSet.add(*exclIter);
         }
       }
    } else if ( simParams->pairInteractionOn ) {
      for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
      {
        int ifep_type = get_fep_type(exclIter->atom1);
        int jfep_type = get_fep_type(exclIter->atom2);
        if ( simParams->pairInteractionSelf ) {
          // for pair-self, both atoms must be in group 1.
          if (ifep_type != 1 || jfep_type != 1) {
            fepExclusionSet.add(*exclIter);
          }
        } else {

          // for pair, must have one from each group.
          if (!(ifep_type == 1 && jfep_type == 2) &&
              !(ifep_type == 2 && jfep_type == 1)) {
            fepExclusionSet.add(*exclIter);
          }
        }
       }
    }

    UniqueSetIter<Exclusion> fepIter(fepExclusionSet);
    for ( fepIter=fepIter.begin(); fepIter != fepIter.end(); fepIter++ )
    {
      exclusionSet.del(*fepIter);
    }
  }
    /*      END OF FUNCTION stripFepExcl      */

#else

//===Memory optimized version of functions that read Molecule file===//
void Molecule::read_mol_signatures(char *fname, Parameters *params, ConfigList *cfgList){
    FILE *psf_file;    //  pointer to .psf file
    int ret_code;    //  ret_code from NAMD_read_line calls
    char buffer[512];

    
    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];

    
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "FORMAT VERSION")) {
        NAMD_die("The compressed psf file format is incorrect, please re-generate!\n");
    }
    float psfVer = 0.0f;
    sscanf(buffer, "FORMAT VERSION: %f\n", &psfVer);
    if(fabs(psfVer - COMPRESSED_PSF_VER)>1e-6) {
        NAMD_die("The compressed psf file format is incorrect, please re-generate!\n");
    }

    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NSEGMENTNAMES"))
        NAMD_die("UNABLE TO FIND NSEGMENTNAMES");
    sscanf(buffer, "%d", &segNamePoolSize);
#if 0
    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);
    }
#else
    for(int i=0; i<segNamePoolSize; i++) NAMD_read_line(psf_file, buffer);
#endif
    
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NRESIDUENAMES"))
        NAMD_die("UNABLE TO FIND NRESIDUENAMES");
    sscanf(buffer, "%d", &resNamePoolSize);
#if 0
    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);
    }
#else
    for(int i=0; i<resNamePoolSize; i++) NAMD_read_line(psf_file, buffer);
#endif

    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);
    }
    
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NATOMTYPES"))
        NAMD_die("UNABLE TO FIND NATOMTYPES");
    sscanf(buffer, "%d", &atomTypePoolSize);
#if 0
    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);
    }
#else
    for(int i=0; i<atomTypePoolSize; i++) NAMD_read_line(psf_file, buffer);
#endif
    
    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);
    }
    
    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);
    }
    
    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++){
        
        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++;
        }

        
        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;
        }
        
        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;
        }
        
        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;
        }
        
        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;
        }
    }

    
    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));

        
        exclSigPool[i].setOffsets(fullExcls, modExcls);

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

    
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NCLUSTERS")) {
        NAMD_die("UNABLE TO FIND NCLUSTERS");
    }
    sscanf(buffer, "%d", &numClusters);

    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NATOM"))
        NAMD_die("UNABLE TO FIND NATOM");
    sscanf(buffer, "%d", &numAtoms);

    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NHYDROGENGROUP"))
        NAMD_die("UNABLE TO FIND NHYDROGENGROUP");
    sscanf(buffer, "%d", &numHydrogenGroups);

    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "MAXHYDROGENGROUPSIZE"))
        NAMD_die("UNABLE TO FIND MAXHYDROGENGROUPSIZE");
    sscanf(buffer, "%d", &maxHydrogenGroupSize);
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "NMIGRATIONGROUP"))
        NAMD_die("UNABLE TO FIND NMIGRATIONGROUP");
    sscanf(buffer, "%d", &numMigrationGroups);
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "MAXMIGRATIONGROUPSIZE"))
        NAMD_die("UNABLE TO FIND MAXMIGRATIONGROUPSIZE");
    sscanf(buffer, "%d", &maxMigrationGroupSize);

    int inputRigidType = -1;
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "RIGIDBONDTYPE"))
      NAMD_die("UNABLE TO FIND RIGIDBONDTYPE");
    sscanf(buffer, "%d", &inputRigidType);
    if(simParams->rigidBonds != RIGID_NONE){
      //check whether the input rigid bond type matches
      if(simParams->rigidBonds != inputRigidType){
        const char *tmpstr[]={"RIGID_NONE", "RIGID_ALL", "RIGID_WATER"};
        char errmsg[125];
        sprintf(errmsg, "RIGIDBOND TYPE MISMATCH BETWEEN INPUT (%s) AND CURRENT RUN (%s)", 
                tmpstr[inputRigidType], tmpstr[simParams->rigidBonds]);
        NAMD_die(errmsg);
      }
    }
#if 0
//    int isOccupancyValid, isBFactorValid;
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "OCCUPANCYVALID"))
        NAMD_die("UNABLE TO FIND OCCUPANCYVALID");
    sscanf(buffer, "%d", &isOccupancyValid);
    NAMD_read_line(psf_file, buffer);
    if(!NAMD_find_word(buffer, "TEMPFACTORVALID"))
        NAMD_die("UNABLE TO FIND TEMPFACTORVALID");
    sscanf(buffer, "%d", &isBFactorValid);    
#endif

    //Just reading for the parameters values; extra Bonds, Dihedrals etc.
    //have been taken into account when compressing the molecule object.
    //The actual number of Bonds, Dihedrals etc. will be calculated based
    //on atom signatures.
    if(cfgList && simParams->extraBondsOn)
        build_extra_bonds(params, cfgList->find("extraBondsFile"));

    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);
}

/*
 * The following method is called on every input processors. However, in SMP mode, two 
 * input procs are likely to be inside the same SMP node. Additionally, there's only one 
 * Molecule object per SMP node. Therefore, if there are any assignments to the molecule 
 * object in this function, there's going to be  DATA RACE !!  That's why the calculation of 
 * numTupes(Bonds, Angles, etc) and numExclusions has to be done in ParallelIOMgr, and 
 * then reduce those values.   
 * -Chao Mei
 */
void Molecule::read_binary_atom_info(int fromAtomID, int toAtomID, InputAtomList& inAtoms){
    int numAtomsPar = toAtomID-fromAtomID+1;
    CmiAssert(numAtomsPar > 0);
    CmiAssert(inAtoms.size() == numAtomsPar);
    
    /*
    //all the following vars are not needed as the
    //atom info is loaded into inAtoms.
    atoms = new AtomCstInfo[numAtomsPar];
    atomNames = new AtomNameIdx[numAtomsPar];
    eachAtomMass = new Index[numAtomsPar];
    eachAtomCharge = new Index[numAtomsPar];
    eachAtomSig = new Index[numAtomsPar];
    eachAtomExclSig = new Index[numAtomsPar];
    */
    /*
    atoms = new AtomCstInfo[numAtomsPar];
    atomNames = new AtomNameIdx[numAtomsPar];
    */

    atoms = NULL;
    atomNames = NULL;
    eachAtomMass = NULL;
    eachAtomCharge = NULL;
    eachAtomSig = NULL;
    eachAtomExclSig = NULL;
    clusterSigs = NULL;

    /* HydrogenGroup is not needed here anymore
    hydrogenGroup.resize(numAtomsPar);
    ResidueLookupElem *tmpResLookup = resLookup;
    */
/*
    if (isOccupancyValid) {
        occupancy = new float[numAtomsPar];
    }
    if (isBFactorValid) {
        bfactor = new float[numAtomsPar];
    }
*/
    occupancy = NULL;
    bfactor = NULL;
    
    char *segment_name; 

    //use "fopen" instead of "Fopen" because "stat" call inside Fopen may 
    //fail on some platforms (such as BG/P) for very large file because of 
    //EOVERFLOW, say a 2GB file. -Chao Mei
    FILE *perAtomFile = fopen(simParams->binAtomFile, "rb");
    if (perAtomFile==NULL) {
        char err_msg[512];
        sprintf(err_msg, "UNABLE TO OPEN THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s", simParams->binAtomFile);
        NAMD_die(err_msg);
    }
    int needFlip = 0;
    int magicNum = COMPRESSED_PSF_MAGICNUM;
    int rMagicNum = COMPRESSED_PSF_MAGICNUM;
    flipNum((char *)&rMagicNum, sizeof(int), 1);
    int fMagicNum;
    fread(&fMagicNum, sizeof(int), 1, perAtomFile);
    if (fMagicNum==magicNum) {
        needFlip = 0;
    } else if (fMagicNum==rMagicNum) {
        needFlip = 1;
    } else {
        char err_msg[512];
        sprintf(err_msg, "THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS CORRUPTED", simParams->binAtomFile);
        NAMD_die(err_msg);
    }

    float verNum =  0.0f;
    fread(&verNum, sizeof(float), 1, perAtomFile);
    if (needFlip) flipNum((char *)&verNum, sizeof(float), 1);
    if (fabs(verNum - COMPRESSED_PSF_VER)>1e-6) {
        char err_msg[512];
        sprintf(err_msg, "THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n", simParams->binAtomFile);
        NAMD_die(err_msg);
    }

    int recSize = 0;
    fread(&recSize, sizeof(int), 1, perAtomFile);
    if(needFlip) flipNum((char *)&recSize, sizeof(int), 1);
    if(recSize != sizeof(OutputAtomRecord)){
      char err_msg[512];
      sprintf(err_msg, "THE ASSOCIATED PER-ATOM RECORD SIZE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n", simParams->binAtomFile);
      NAMD_die(err_msg);
    }
    
    const int BUFELEMS = 32*1024; //32K elems

    //remember to convert to long in case of int overflow!
    int64 startbyte=((int64)fromAtomID)*sizeof(OutputAtomRecord);
#ifdef WIN32
    if ( _fseeki64(perAtomFile,startbyte,SEEK_CUR) )
#else
    if ( fseeko(perAtomFile,startbyte,SEEK_CUR) )
#endif
    {
      char errmsg[512];
      sprintf(errmsg, "Error on seeking binary file %s", simParams->binAtomFile);
      NAMD_err(errmsg);
    }

    //reduce the number of fread calls as file I/O is expensive.
    OutputAtomRecord *elemsBuf = new OutputAtomRecord[BUFELEMS];
    int atomsCnt = numAtomsPar;
    int curIdx=0;
    OutputAtomRecord *oneRec = NULL;
    while(atomsCnt >= BUFELEMS) {
      if ( fread((char *)elemsBuf, sizeof(OutputAtomRecord), BUFELEMS, perAtomFile) != BUFELEMS ) {
        char errmsg[512];
        sprintf(errmsg, "Error on reading binary file %s", simParams->binAtomFile);
        NAMD_err(errmsg);
      }
      oneRec = elemsBuf;
      for(int i=0; i<BUFELEMS; i++, curIdx++, oneRec++) {
        InputAtom *fAtom = &(inAtoms[curIdx]);
        int aid = curIdx+fromAtomID;
        if(needFlip) oneRec->flip();
        load_one_inputatom(aid, oneRec, fAtom);        
      }
      atomsCnt -= BUFELEMS;
    }

    if ( fread(elemsBuf, sizeof(OutputAtomRecord), atomsCnt, perAtomFile) != atomsCnt ) {
      char errmsg[512];
      sprintf(errmsg, "Error on reading binary file %s", simParams->binAtomFile);
      NAMD_err(errmsg);
    }
    oneRec = elemsBuf;    
    for(int i=curIdx; i<numAtomsPar; i++, oneRec++) {
      InputAtom *fAtom = &(inAtoms[i]);
      int aid = i+fromAtomID;
      if(needFlip) oneRec->flip();
      load_one_inputatom(aid,oneRec,fAtom);      
    }

    if ( fclose(perAtomFile) ) {
      char errmsg[512];
      sprintf(errmsg, "Error on closing binary file %s", simParams->binAtomFile);
      NAMD_err(errmsg);
    }

    delete [] elemsBuf;

    //deal with fixed atoms info
    if(simParams->fixedAtomsOn){
        int listIdx=0;
        is_atom_fixed(fromAtomID, &listIdx);
        for(int i=listIdx; i<fixedAtomsSet->size(); i++){
            const AtomSet one = fixedAtomsSet->item(i);
            //set the atoms in this range to be fixed
            int sAtomId = one.aid1>fromAtomID ? one.aid1:fromAtomID;
            int eAtomId = one.aid2>toAtomID? toAtomID:one.aid2;
            for(int j=sAtomId; j<=eAtomId; j++)
                inAtoms[j-fromAtomID].atomFixed = 1;
        }
    }
}

void Molecule::load_one_inputatom(int aid, OutputAtomRecord *one, InputAtom *fAtom){

  char *thisAtomName = NULL;
  fAtom->isValid=true;
  
  //segment_name = segNamePool[sIdx[0]];
  /*
  atomNames[i].resnameIdx = sIdx[1];
  atomNames[i].atomnameIdx = sIdx[2];
  atomNames[i].atomtypeIdx = sIdx[3]; 
  */
  thisAtomName = atomNamePool[one->sSet.atomNameIdx];
         
  fAtom->charge = atomChargePool[one->sSet.chargeIdx];
  fAtom->mass = atomMassPool[one->sSet.massIdx];
  // Using double precision division for reciprocal mass.
  fAtom->recipMass = ( fAtom->mass > 0 ? (1. / fAtom->mass) : 0 );
  fAtom->sigId = one->iSet.atomSigIdx;
  fAtom->exclId = one->iSet.exclSigIdx;
  fAtom->vdwType = one->sSet.vdw_type;
  
  //atoms[i].vdw_type = sIdx[8];
  
  int residue_number; //for residue number
  residue_number = one->iSet.resID;
  
  /*
  atoms[i].partner = iIdx[2];
  atoms[i].hydrogenList= iIdx[3];
  */
  
  fAtom->id=aid;
  fAtom->atomFixed = 0;
  fAtom->hydList = one->iSet.hydrogenList;
  fAtom->hydrogenGroupSize=one->iSet.atomsInGroup;
  fAtom->GPID=one->iSet.GPID;        
  //fAtom->waterVal=one->waterVal;
  fAtom->migrationGroupSize=one->iSet.atomsInMigrationGroup;
  fAtom->MPID=one->iSet.MPID;
  fAtom->isGP=(fAtom->hydrogenGroupSize ? 1 : 0);
  fAtom->isMP=( fAtom->migrationGroupSize ? 1 : 0 ); 
  
  if(simParams->rigidBonds) {
    fAtom->rigidBondLength = one->fSet.rigidBondLength;
  }else{
    fAtom->rigidBondLength = 0.0;
  }
  
  //Node::Object()->ioMgr->maxAtom=fAtom.id;        
  
  /*if (isOccupancyValid)
      occupancy[i] = tmpf[0];
  if (isBFactorValid)
      bfactor[i] = tmpf[1];*/
  
  /*
  if (tmpResLookup) tmpResLookup =
          tmpResLookup->append(segment_name, residue_number, i);
  */
  
  Real thisAtomMass = fAtom->mass;
  
  if ( simParams->ignoreMass ) {
  } else if (thisAtomMass <= 0.05) {
      fAtom->status |= LonepairAtom;
  } else if (thisAtomMass < 1.0) {
      fAtom->status |= DrudeAtom;
  } else if (thisAtomMass <= 3.5) {
      fAtom->status = HydrogenAtom;
  } else if (thisAtomName[0]=='O' &&
             (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)) {
      fAtom->status = OxygenAtom;
  }
  
  //Move the langevinParam setting which depends on atom's status
  //to the time when each home patch is filled with their atoms
  //(WorkDistrib::fillAtomListForOnePatch so that "langevinParam"
  //could be shared with the "hydVal".
  //--Chao Mei 
}

//Well, the exclusion check signatures could also done on PE0 and
//sent to other processors through send_Molecule/receive_Molecule 
//two procedures.
void Molecule::build_excl_check_signatures(){
   exclChkSigPool = new ExclusionCheck[exclSigPoolSize];
   for(int i=0; i<exclSigPoolSize; i++){
       ExclusionSignature *sig = &exclSigPool[i];
       ExclusionCheck *sigChk = &exclChkSigPool[i];
       if(sig->fullExclCnt){
           if(!sig->modExclCnt){ //only having fullExclusion
               sigChk->min = sig->fullOffset[0];
               sigChk->max = sig->fullOffset[sig->fullExclCnt-1];
           }else{ //have both full and modified exclusion
               int fullMin, fullMax, modMin, modMax;
               
               fullMin = sig->fullOffset[0];
               fullMax = sig->fullOffset[sig->fullExclCnt-1];
           
               modMin = sig->modOffset[0];
               modMax = sig->modOffset[sig->modExclCnt-1];
               
               if(fullMin < modMin)
                   sigChk->min = fullMin;
               else
                   sigChk->min = modMin;
               if(fullMax < modMax)
                   sigChk->max = modMax;
               else
                   sigChk->max = fullMax;
           }        
       }else{
           if(sig->modExclCnt){
               sigChk->min = sig->modOffset[0];
               sigChk->max = sig->modOffset[sig->modExclCnt-1];
           }else{ //both count are 0
               if(CkMyPe()==0)
                   iout << iWARN << "an empty exclusion signature with index "
                     << i << "!\n" << endi;
               continue;
           }
       }           

       sigChk->flags = new char[sigChk->max-sigChk->min+1];
       memset(sigChk->flags, 0, sizeof(char)*(sigChk->max-sigChk->min+1));
       for(int j=0; j<sig->fullExclCnt; j++){
           int dist = sig->fullOffset[j] - sigChk->min;
           sigChk->flags[dist] = EXCHCK_FULL;
       }
       for(int j=0; j<sig->modExclCnt; j++){
           int dist = sig->modOffset[j] - sigChk->min;
           sigChk->flags[dist] = EXCHCK_MOD;
       }
   }
}

void Molecule::load_atom_set(StringList *setfile, const char *setname,
        int *numAtomsInSet, AtomSetList **atomsSet) const {
  if(setfile == NULL) {
    char errmsg[128];
    sprintf(errmsg,"The text input file for %s atoms is not found!", setname);
    NAMD_die(errmsg);
  }
  FILE *ifp = fopen(setfile->data, "r");
  
  if(ifp==NULL){
      char errmsg[128];
      sprintf(errmsg, "ERROR IN OPENING %s ATOMS FILE: %s\n", setname, setfile->data);
      NAMD_die(errmsg);
  }

  char oneline[128];
  int numLocalAtoms = 0;
  AtomSet one;
  AtomSetList *localAtomsSet = new AtomSetList();  
  while(1) {
    int ret = NAMD_read_line(ifp, oneline, 128);
    if(ret!=0) break;
    if(NAMD_blank_string(oneline)) continue;
    bool hasDash = false;
    for(int i=0; oneline[i] && i<128; i++){
      if(oneline[i]=='-') {
        hasDash = true;
        break;
      }
    }
    if(hasDash) {
      sscanf(oneline,"%d-%d", &(one.aid1), &(one.aid2));
      if(one.aid1>one.aid2 || one.aid1<0 || one.aid2<0) {
        char errmsg[512];
        sprintf(errmsg, "The input for %s atoms is wrong: %s\n", setname, oneline);
        NAMD_die(errmsg);
      }
      numLocalAtoms += (one.aid2-one.aid1+1);
    }else{
      sscanf(oneline, "%d", &(one.aid1));
      if(one.aid1<0) {
        char errmsg[512];
        sprintf(errmsg, "The input for %s atoms is wrong: %s\n", setname, oneline);
        NAMD_die(errmsg);      
      }
      one.aid2 = one.aid1;
      numLocalAtoms++;
    }
    localAtomsSet->add(one);
  }
  //sort the localAtomsSet for binary search to decide 
  //whether an atom is in the set or not
  std::sort(localAtomsSet->begin(), localAtomsSet->end());  

  *numAtomsInSet = numLocalAtoms;
  *atomsSet = localAtomsSet;
}

void Molecule::load_fixed_atoms(StringList *fixedfile){
  load_atom_set(fixedfile, "FIXED", &numFixedAtoms, &fixedAtomsSet);
}

void Molecule::load_constrained_atoms(StringList *constrainedfile){
  load_atom_set(constrainedfile, "CONSTRAINED", &numConstraints, &constrainedAtomsSet);
}

Bool Molecule::is_atom_in_set(AtomSetList *localAtomsSet, int aid, int *listIdx) const {
  int idx = localAtomsSet->size();
  int rIdx = 0;
  int lIdx = localAtomsSet->size()-1;
  
  while(rIdx <= lIdx){
    int mIdx = (rIdx+lIdx)/2;
    const AtomSet one = localAtomsSet->item(mIdx);

    if(aid < one.aid1){
      //aid could be in [rIdx, mIdx);
      idx = mIdx;
      lIdx = mIdx-1;
    }else if(aid > one.aid1){
      //aid could be inside the atom set "one" or in (mIdx, lIdx];
      if(aid<=one.aid2){
        //found, aid in the atom set "one"
        if(listIdx) *listIdx = mIdx;
        return 1;
      }else{
        rIdx = mIdx+1;
      }
    }else{
      //found, aid is exactly same with one.aid1
      if(listIdx) *listIdx = mIdx;
      return 1;
    }
  }

  //not found
  if(listIdx) *listIdx = idx;
  return 0;
}

#endif


/************************************************************************/
/*                  */
/*      FUNCTION print_atoms        */
/*                  */
/*  print_atoms prints out the list of atoms stored in this object. */
/*  It is inteded mainly for debugging purposes.      */
/*                  */
/************************************************************************/

void Molecule::print_atoms(Parameters *params)
{
#ifdef MEM_OPT_VERSION
    DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
#else   
  register int i;
  Real sigma;
  Real epsilon;
  Real sigma14;
  Real epsilon14;

  DebugM(2,"ATOM LIST\n" \
      << "******************************************\n" \
                  << "NUM  NAME TYPE RES  MASS    CHARGE CHARGE   FEP-CHARGE"  \
      << "SIGMA   EPSILON SIGMA14 EPSILON14\n" \
        << endi);

  for (i=0; i<numAtoms; i++)
  {
    params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 
        atoms[i].vdw_type);

    DebugM(2,i+1 << " " << atomNames[i].atomname  \
              << " " << atomNames[i].atomtype << " " \
              << atomNames[i].resname  << " " << atoms[i].mass  \
        << " " << atoms[i].charge << " " << sigma \
        << " " << epsilon << " " << sigma14 \
        << " " << epsilon14 << "\n" \
        << endi);
  }
#endif  
}
/*      END OF FUNCTION print_atoms      */

/************************************************************************/
/*                  */
/*      FUNCTION print_bonds        */
/*                  */
/*  print_bonds prints out the list of bonds stored in this object. */
/*  It is inteded mainly for debugging purposes.      */
/*                  */
/************************************************************************/

void Molecule::print_bonds(Parameters *params)
{
#ifdef MEM_OPT_VERSION
    DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
#else   
  register int i;
  Real k;
  Real x0;

  DebugM(2,"BOND LIST\n" << "********************************\n" \
      << "ATOM1 ATOM2 TYPE1 TYPE2      k        x0" \
      << endi);

  for (i=0; i<numBonds; i++)
  {
    params->get_bond_params(&k, &x0, bonds[i].bond_type);

    DebugM(2,bonds[i].atom1+1 << " " \
       << bonds[i].atom2+1 << " "   \
       << atomNames[bonds[i].atom1].atomtype << " "  \
       << atomNames[bonds[i].atom2].atomtype << " " << k \
       << " " << x0 << endi);
  }
  
#endif  
}
/*      END OF FUNCTION print_bonds      */

/************************************************************************/
/*                  */
/*      FUNCTION print_exclusions      */
/*                  */
/*  print_exlcusions prints out the list of exlcusions stored in    */
/*  this object.  It is inteded mainly for debugging purposes.    */
/*                  */
/************************************************************************/

void Molecule::print_exclusions()
{
#ifdef MEM_OPT_VERSION
    DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
#else
  register int i;

  DebugM(2,"EXPLICIT EXCLUSION LIST\n" \
      << "********************************\n" \
            << "ATOM1 ATOM2 " \
      << endi);

  for (i=0; i<numExclusions; i++)
  {
    DebugM(2,exclusions[i].atom1+1 << "  " \
       << exclusions[i].atom2+1 << endi);
  }
#endif
}
/*      END OF FUNCTION print_exclusions    */

/************************************************************************/
/*                  */
/*      FUNCTION send_Molecule        */
/*                  */
/*  send_Molecule is used by the Master node to distribute the      */
/*   structural information to all the client nodes.  It is NEVER called*/
/*   by the client nodes.              */
/*                  */
/************************************************************************/

void Molecule::send_Molecule(MOStream *msg){
#ifdef MEM_OPT_VERSION
//in the memory optimized version, only the atom signatures are broadcast
//to other Nodes. --Chao Mei

  msg->put(numAtoms);

  msg->put(massPoolSize);
  msg->put(massPoolSize, atomMassPool);

  msg->put(chargePoolSize);
  msg->put(chargePoolSize, atomChargePool); 

  //put atoms' signatures
  msg->put(atomSigPoolSize);
  for(int i=0; i<atomSigPoolSize; i++)
      atomSigPool[i].pack(msg);

  //put atom's exclusion signatures
  msg->put(exclSigPoolSize);
  for(int i=0; i<exclSigPoolSize; i++)
      exclSigPool[i].pack(msg);

  msg->put(numHydrogenGroups);      
  msg->put(maxHydrogenGroupSize);      
  msg->put(numMigrationGroups);      
  msg->put(maxMigrationGroupSize);            
  msg->put(isOccupancyValid);
  msg->put(isBFactorValid);
  
  //put names for atoms
  msg->put(atomNamePoolSize);
  for(int i=0; i<atomNamePoolSize;i++) {
    int len = strlen(atomNamePool[i]);
    msg->put(len);
    msg->put(len*sizeof(char), atomNamePool[i]);
  } 
  
  if(simParams->fixedAtomsOn){
    int numFixedAtomsSet = fixedAtomsSet->size();
    msg->put(numFixedAtoms);
    msg->put(numFixedAtomsSet);
    msg->put(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
  }

  if (simParams->constraintsOn) {
    int numConstrainedAtomsSet = constrainedAtomsSet->size();
    msg->put(numConstraints);
    msg->put(numConstrainedAtomsSet);
    msg->put(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
  }
    
#else
  msg->put(numAtoms);
  msg->put(numAtoms*sizeof(Atom), (char*)atoms);
  
  //  Send the bond information
  msg->put(numRealBonds);
  msg->put(numBonds);
 
  if (numBonds)
  {
    msg->put(numBonds*sizeof(Bond), (char*)bonds);
  }

  //  Send the angle information
  msg->put(numAngles);  
  if (numAngles)
  {
    msg->put(numAngles*sizeof(Angle), (char*)angles);
  }  

  //  Send the dihedral information
  msg->put(numDihedrals);
  if (numDihedrals)
  {
    msg->put(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
  }  

  if (simParams->sdScaling) {
    msg->put(num_alch_unpert_Bonds);
    msg->put(num_alch_unpert_Bonds*sizeof(Bond), (char*)alch_unpert_bonds);

    msg->put(num_alch_unpert_Angles);
    msg->put(num_alch_unpert_Angles*sizeof(Angle), (char*)alch_unpert_angles);

    msg->put(num_alch_unpert_Dihedrals);
    msg->put(num_alch_unpert_Dihedrals*sizeof(Dihedral), (char*)alch_unpert_dihedrals);
  }

  //  Send the improper information
  msg->put(numImpropers);  
  if (numImpropers)
  {
    msg->put(numImpropers*sizeof(Improper), (char*)impropers);
  }

  //  Send the crossterm information
  msg->put(numCrossterms);
  if (numCrossterms)
  {
    msg->put(numCrossterms*sizeof(Crossterm), (char*)crossterms);
  }

  // send the hydrogen bond donor information
  msg->put(numDonors);
  if(numDonors)
  {
    msg->put(numDonors*sizeof(Bond), (char*)donors);
  }

  // send the hydrogen bond acceptor information
  msg->put(numAcceptors);
  if(numAcceptors)
  {
    msg->put(numAcceptors*sizeof(Bond), (char*)acceptors);
  }

  //  Send the exclusion information  
  msg->put(numExclusions);
  if (numExclusions)
  {
    msg->put(numExclusions*sizeof(Exclusion), (char*)exclusions);
  }      
  //  Send the constraint information, if used
  if (simParams->constraintsOn)
  {
     msg->put(numConstraints);
     
     msg->put(numAtoms, consIndexes);
     
     if (numConstraints)
     {
       msg->put(numConstraints*sizeof(ConstraintParams), (char*)consParams);
     }
  }
#endif
  
  /* BEGIN gf */
  // Send the gridforce information, if used
  if (simParams->mgridforceOn)
  {
    DebugM(3, "Sending gridforce info\n" << endi);
    msg->put(numGridforceGrids);
    
    for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
      msg->put(numGridforces[gridnum]);
      msg->put(numAtoms, gridfrcIndexes[gridnum]);
      if (numGridforces[gridnum])
      {
       msg->put(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
      }
      GridforceGrid::pack_grid(gridfrcGrid[gridnum], msg);
    }
  }
  /* END gf */
  
  //  Send the stirring information, if used
  if (simParams->stirOn)
  {
     //CkPrintf ("DEBUG: putting numStirredAtoms..\n");
     msg->put(numStirredAtoms);
     //CkPrintf ("DEBUG: putting numAtoms,stirIndexes.. numAtoms=%d\n",numStirredAtoms);
     msg->put(numAtoms, stirIndexes);
     //CkPrintf ("DEBUG: if numStirredAtoms..\n");
     if (numStirredAtoms)
     {
       //CkPrintf ("DEBUG: big put, with (char*)stirParams\n");
       msg->put(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
     }
  }
  
  
  //  Send the moving drag information, if used
  if (simParams->movDragOn) {
     msg->put(numMovDrag);
     msg->put(numAtoms, movDragIndexes);
     if (numMovDrag)
     {
       msg->put(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
     }
  }
  
  //  Send the rotating drag information, if used
  if (simParams->rotDragOn) {
     msg->put(numRotDrag);
     msg->put(numAtoms, rotDragIndexes);
     if (numRotDrag)
     {
       msg->put(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
     }
  }
  
  //  Send the "constant" torque information, if used
  if (simParams->consTorqueOn) {
     msg->put(numConsTorque);
     msg->put(numAtoms, consTorqueIndexes);
     if (numConsTorque)
     {
       msg->put(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
     }
  }
  
  // Send the constant force information, if used
  if (simParams->consForceOn)
  { msg->put(numConsForce);
    msg->put(numAtoms, consForceIndexes);
    if (numConsForce)
      msg->put(numConsForce*sizeof(Vector), (char*)consForce);
  }
  
  if (simParams->excludeFromPressure) {
    msg->put(numExPressureAtoms);
    msg->put(numAtoms, exPressureAtomFlags);
  }
  
#ifndef MEM_OPT_VERSION
  //  Send the langevin parameters, if active
  if (simParams->langevinOn || simParams->tCoupleOn)
  {
    msg->put(numAtoms, langevinParams);
  }
  
  //  Send fixed atoms, if active
  if (simParams->fixedAtomsOn)
  {
    msg->put(numFixedAtoms);
    msg->put(numAtoms, fixedAtomFlags);
  msg->put(numFixedRigidBonds);
  }
  
  if (simParams->qmForcesOn)
  {
    msg->put(numAtoms, qmAtomGroup);
    msg->put(qmNumQMAtoms);
    msg->put(qmNumQMAtoms, qmAtmChrg);
    msg->put(qmNumQMAtoms, qmAtmIndx);
    msg->put(qmNoPC);
    msg->put(qmNumBonds);
    msg->put(qmMeNumBonds);
    msg->put(qmMeNumBonds, qmMeMMindx);
    msg->put(qmMeNumBonds, qmMeQMGrp);
    msg->put(qmPCFreq);
    msg->put(qmNumGrps);
    msg->put(qmNumGrps, qmGrpID);
    msg->put(qmNumGrps, qmCustPCSizes);
    msg->put(qmTotCustPCs);
    msg->put(qmTotCustPCs, qmCustomPCIdxs);
  }
  
  //fepb
  // send fep atom info
  if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
    msg->put(numFepInitial);
    msg->put(numFepFinal);
    msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
  }
  //fepe

  if (simParams->soluteScalingOn) {
    msg->put(numAtoms*sizeof(char), (char*)ssAtomFlags);
    msg->put(ss_num_vdw_params);
    msg->put(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
    msg->put(numAtoms*sizeof(int), (char*)ss_index);
  }

  #ifdef OPENATOM_VERSION
  // needs to be refactored into its own openatom version
  if (simParams->openatomOn ) {
    msg->put(numFepInitial);
    msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
  }
  #endif //OPENATOM_VERSION
  
  // DRUDE: send data read from PSF
  msg->put(is_lonepairs_psf);
  if (is_lonepairs_psf) {
    msg->put(numLphosts);
    msg->put(numLphosts*sizeof(Lphost), (char*)lphosts);
  }
  msg->put(is_drude_psf);
  if (is_drude_psf) {
    msg->put(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
    msg->put(numTholes);
    msg->put(numTholes*sizeof(Thole), (char*)tholes);
    msg->put(numAnisos);
    msg->put(numAnisos*sizeof(Aniso), (char*)anisos);
  }
  msg->put(numZeroMassAtoms);
  // DRUDE

  //LCPO
  if (simParams->LCPOOn) {
    msg->put(numAtoms, (int*)lcpoParamType);
  }
  
  //Send GromacsPairStuff -- JLai
  if (simParams->goGroPair) {
    msg->put(numLJPair);
    msg->put(numLJPair,indxLJA);
    msg->put(numLJPair,indxLJB);
    msg->put(numLJPair,pairC6);
    msg->put(numLJPair,pairC12);
    msg->put(numLJPair,gromacsPair_type);
    msg->put((numAtoms),pointerToLJBeg);
    msg->put((numAtoms),pointerToLJEnd);
    msg->put(numGaussPair);
    msg->put(numGaussPair,indxGaussA);
    msg->put(numGaussPair,indxGaussB);
    msg->put(numGaussPair,gA);
    msg->put(numGaussPair,gMu1);
    msg->put(numGaussPair,giSigma1);
    msg->put(numGaussPair,gMu2);
    msg->put(numGaussPair,giSigma2);
    msg->put(numGaussPair,gRepulsive);
    msg->put((numAtoms),pointerToGaussBeg);
    msg->put((numAtoms),pointerToGaussEnd);
  }
#endif

  // Broadcast the message to the other nodes
  msg->end();
  delete msg;

#ifdef MEM_OPT_VERSION

  build_excl_check_signatures();

  //set num{Calc}Tuples(Bonds,...,Impropers) to 0
  numBonds = numCalcBonds = 0;
  numAngles = numCalcAngles = 0;
  numDihedrals = numCalcDihedrals = 0;
  numImpropers = numCalcImpropers = 0;
  numCrossterms = numCalcCrossterms = 0;
  numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
  // JLai
  numLJPair = numCalcLJPair = 0;
  // End of JLai

#else

  //  Now build arrays of indexes into these arrays by atom      
  build_lists_by_atom();

#endif
}
 /*      END OF FUNCTION send_Molecule      */

    /************************************************************************/
    /*                  */
    /*      FUNCTION receive_Molecule      */
    /*                  */
    /*  receive_Molecule is used by all the clients to receive the  */
    /*   structural data sent out by the master node.  It is NEVER called   */
    /*   by the Master node.            */
    /*                  */
    /************************************************************************/

void Molecule::receive_Molecule(MIStream *msg){
  //  Get the atom information
  msg->get(numAtoms);

#ifdef MEM_OPT_VERSION
//in the memory optimized version, only the atom signatures are recved
//from the master Node. --Chao Mei

  msg->get(massPoolSize);
  if(atomMassPool) delete [] atomMassPool;
  atomMassPool = new Real[massPoolSize];
  msg->get(massPoolSize, atomMassPool);

  msg->get(chargePoolSize);
  if(atomChargePool) delete [] atomChargePool;
  atomChargePool = new Real[chargePoolSize];
  msg->get(chargePoolSize, atomChargePool);

  //get atoms' signatures
  msg->get(atomSigPoolSize);
  if(atomSigPool) delete [] atomSigPool;
  atomSigPool = new AtomSignature[atomSigPoolSize];
  for(int i=0; i<atomSigPoolSize; i++)
      atomSigPool[i].unpack(msg);

  //get exclusions' signatures
  msg->get(exclSigPoolSize);
  if(exclSigPool) delete [] exclSigPool;
  exclSigPool = new ExclusionSignature[exclSigPoolSize];
  for(int i=0; i<exclSigPoolSize; i++)
      exclSigPool[i].unpack(msg);
 
  msg->get(numHydrogenGroups);      
  msg->get(maxHydrogenGroupSize);      
  msg->get(numMigrationGroups);      
  msg->get(maxMigrationGroupSize);      
  msg->get(isOccupancyValid);
  msg->get(isBFactorValid);

   //get names for atoms
  msg->get(atomNamePoolSize);
  atomNamePool = new char *[atomNamePoolSize];
  for(int i=0; i<atomNamePoolSize;i++) {
    int len;
    msg->get(len);
    atomNamePool[i] = nameArena->getNewArray(len+1);
    msg->get(len, atomNamePool[i]);
  }
  
  if(simParams->fixedAtomsOn){
    int numFixedAtomsSet;
    msg->get(numFixedAtoms);
    msg->get(numFixedAtomsSet);
    fixedAtomsSet = new AtomSetList(numFixedAtomsSet);
    msg->get(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
  } 

  if(simParams->constraintsOn){
    int numConstrainedAtomsSet;
    msg->get(numConstraints);
    msg->get(numConstrainedAtomsSet);
    constrainedAtomsSet = new AtomSetList(numConstrainedAtomsSet);
    msg->get(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
  } 

#else
  delete [] atoms;
  atoms= new Atom[numAtoms];  
  msg->get(numAtoms*sizeof(Atom), (char*)atoms);

  //  Get the bond information
  msg->get(numRealBonds);
  msg->get(numBonds);    
  if (numBonds)
  {
    delete [] bonds;
    bonds=new Bond[numBonds]; 
    msg->get(numBonds*sizeof(Bond), (char*)bonds);
  }  

  //  Get the angle information
  msg->get(numAngles);  
  if (numAngles)
  {
    delete [] angles;
    angles=new Angle[numAngles];  
    msg->get(numAngles*sizeof(Angle), (char*)angles);
  }  

  //  Get the dihedral information
  msg->get(numDihedrals);    
  if (numDihedrals)
  {
    delete [] dihedrals;
    dihedrals=new Dihedral[numDihedrals];  
    msg->get(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
  }  

  if (simParams->sdScaling) {
    msg->get(num_alch_unpert_Bonds);
    alch_unpert_bonds=new Bond[num_alch_unpert_Bonds];
    msg->get(num_alch_unpert_Bonds*sizeof(Bond), (char*)alch_unpert_bonds);

    msg->get(num_alch_unpert_Angles);
    alch_unpert_angles=new Angle[num_alch_unpert_Angles];
    msg->get(num_alch_unpert_Angles*sizeof(Angle), (char*)alch_unpert_angles);

    msg->get(num_alch_unpert_Dihedrals);
    alch_unpert_dihedrals=new Dihedral[num_alch_unpert_Dihedrals];
    msg->get(num_alch_unpert_Dihedrals*sizeof(Dihedral), (char*)alch_unpert_dihedrals);
  }
  
  //  Get the improper information
  msg->get(numImpropers);
  if (numImpropers)
  {
    delete [] impropers;
    impropers=new Improper[numImpropers];  
    msg->get(numImpropers*sizeof(Improper), (char*)impropers);
  }
  
  //  Get the crossterm information
  msg->get(numCrossterms);
  if (numCrossterms)
  {
    delete [] crossterms;
    crossterms=new Crossterm[numCrossterms];  
    msg->get(numCrossterms*sizeof(Crossterm), (char*)crossterms);
  }
  
  //  Get the hydrogen bond donors
  msg->get(numDonors);  
  if (numDonors)
  {
    delete [] donors;
    donors=new Bond[numDonors];  
    msg->get(numDonors*sizeof(Bond), (char*)donors);
  }
  
  //  Get the hydrogen bond acceptors
  msg->get(numAcceptors);  
  if (numAcceptors)
  {
    delete [] acceptors;
    acceptors=new Bond[numAcceptors];  
    msg->get(numAcceptors*sizeof(Bond), (char*)acceptors);
  }
  
  //  Get the exclusion information 
  msg->get(numExclusions);  
  if (numExclusions)
  {
    delete [] exclusions;
    exclusions=new Exclusion[numExclusions];  
    msg->get(numExclusions*sizeof(Exclusion), (char*)exclusions);
  }
        
      //  Get the constraint information, if they are active
      if (simParams->constraintsOn)
      {
         msg->get(numConstraints);

         delete [] consIndexes;
         consIndexes = new int32[numAtoms];
         
         msg->get(numAtoms, consIndexes);
         
         if (numConstraints)
         {
           delete [] consParams;
           consParams = new ConstraintParams[numConstraints];
      
           msg->get(numConstraints*sizeof(ConstraintParams), (char*)consParams);
         }
      }
#endif

      /* BEGIN gf */
      if (simParams->mgridforceOn)
      {
         DebugM(3, "Receiving gridforce info\n");
         
         msg->get(numGridforceGrids);
         
         DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
         
         delete [] numGridforces;
         numGridforces = new int[numGridforceGrids];
         
         delete [] gridfrcIndexes;      // Should I be deleting elements of these first?
         delete [] gridfrcParams;
         delete [] gridfrcGrid;
         gridfrcIndexes = new int32*[numGridforceGrids];
         gridfrcParams = new GridforceParams*[numGridforceGrids];
         gridfrcGrid = new GridforceGrid*[numGridforceGrids];
         
         int grandTotalGrids = 0;
         for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
             msg->get(numGridforces[gridnum]);
             
             gridfrcIndexes[gridnum] = new int32[numAtoms];
             msg->get(numAtoms, gridfrcIndexes[gridnum]);
         
             if (numGridforces[gridnum])
             {
                 gridfrcParams[gridnum] = new GridforceParams[numGridforces[gridnum]];
                 msg->get(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
             }
             
             gridfrcGrid[gridnum] = GridforceGrid::unpack_grid(gridnum, msg);
             
             grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
         }
      }
      /* END gf */
      
      //  Get the stirring information, if stirring is  active
      if (simParams->stirOn)
      {
         msg->get(numStirredAtoms);

         delete [] stirIndexes;
         stirIndexes = new int32[numAtoms];
         
         msg->get(numAtoms, stirIndexes);
         
         if (numStirredAtoms)
         {
           delete [] stirParams;
           stirParams = new StirParams[numStirredAtoms];
      
           msg->get(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
         }
      }
      
      //  Get the moving drag information, if it is active
      if (simParams->movDragOn) {
         msg->get(numMovDrag);
         delete [] movDragIndexes;
         movDragIndexes = new int32[numAtoms];
         msg->get(numAtoms, movDragIndexes);
         if (numMovDrag)
         {
           delete [] movDragParams;
           movDragParams = new MovDragParams[numMovDrag];
           msg->get(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
         }
      }
      
      //  Get the rotating drag information, if it is active
      if (simParams->rotDragOn) {
         msg->get(numRotDrag);
         delete [] rotDragIndexes;
         rotDragIndexes = new int32[numAtoms];
         msg->get(numAtoms, rotDragIndexes);
         if (numRotDrag)
         {
           delete [] rotDragParams;
           rotDragParams = new RotDragParams[numRotDrag];
           msg->get(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
         }
      }
      
      //  Get the "constant" torque information, if it is active
      if (simParams->consTorqueOn) {
         msg->get(numConsTorque);
         delete [] consTorqueIndexes;
         consTorqueIndexes = new int32[numAtoms];
         msg->get(numAtoms, consTorqueIndexes);
         if (numConsTorque)
         {
           delete [] consTorqueParams;
           consTorqueParams = new ConsTorqueParams[numConsTorque];
           msg->get(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
         }
      }
      
      // Get the constant force information, if it's active
      if (simParams->consForceOn)
      { msg->get(numConsForce);
        delete [] consForceIndexes;
        consForceIndexes = new int32[numAtoms];
        msg->get(numAtoms, consForceIndexes);
        if (numConsForce)
        { delete [] consForce;
          consForce = new Vector[numConsForce];
          msg->get(numConsForce*sizeof(Vector), (char*)consForce);
        }
      }

      if (simParams->excludeFromPressure) {
        exPressureAtomFlags = new int32[numAtoms];
        msg->get(numExPressureAtoms);
        msg->get(numAtoms, exPressureAtomFlags);
      }

#ifndef MEM_OPT_VERSION
      //  Get the langevin parameters, if they are active
      if (simParams->langevinOn || simParams->tCoupleOn)
      {
        delete [] langevinParams;
        langevinParams = new Real[numAtoms];

        msg->get(numAtoms, langevinParams);
      }

      //  Get the fixed atoms, if they are active
      if (simParams->fixedAtomsOn)
      {
        delete [] fixedAtomFlags;
        fixedAtomFlags = new int32[numAtoms];

        msg->get(numFixedAtoms);
        msg->get(numAtoms, fixedAtomFlags);
        msg->get(numFixedRigidBonds);
      }

      if (simParams->qmForcesOn)
      {
        if( qmAtomGroup != 0)
            delete [] qmAtomGroup;
        qmAtomGroup = new Real[numAtoms];
        
        msg->get(numAtoms, qmAtomGroup);
        
        msg->get(qmNumQMAtoms);
        
        if( qmAtmChrg != 0)
            delete [] qmAtmChrg;
        qmAtmChrg = new Real[qmNumQMAtoms];
        
        msg->get(qmNumQMAtoms, qmAtmChrg);
        
        if( qmAtmIndx != 0)
            delete [] qmAtmIndx;
        qmAtmIndx = new int[qmNumQMAtoms];
        
        msg->get(qmNumQMAtoms, qmAtmIndx);
        
        msg->get(qmNoPC);
        
        msg->get(qmNumBonds);
        
        msg->get(qmMeNumBonds);
        
        if( qmMeMMindx != 0)
            delete [] qmMeMMindx;
        qmMeMMindx = new int[qmMeNumBonds];
        
        msg->get(qmMeNumBonds, qmMeMMindx);
        
        if( qmMeQMGrp != 0)
            delete [] qmMeQMGrp;
        qmMeQMGrp = new Real[qmMeNumBonds];
        
        msg->get(qmMeNumBonds, qmMeQMGrp);
        
        msg->get(qmPCFreq);
        
        msg->get(qmNumGrps);
        
        if( qmGrpID != 0)
            delete [] qmGrpID;
        qmGrpID = new Real[qmNumGrps];
        msg->get(qmNumGrps, qmGrpID);
        
        if( qmCustPCSizes != 0)
            delete [] qmCustPCSizes;
        qmCustPCSizes = new int[qmNumGrps];
        msg->get(qmNumGrps, qmCustPCSizes);
        
        msg->get(qmTotCustPCs);
        
        if( qmCustomPCIdxs != 0)
            delete [] qmCustomPCIdxs;
        qmCustomPCIdxs = new int[qmTotCustPCs];
        msg->get(qmTotCustPCs, qmCustomPCIdxs);
      }
    
//fepb
      //receive fep atom info
      if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
        delete [] fepAtomFlags;
        fepAtomFlags = new unsigned char[numAtoms];

        msg->get(numFepInitial);
        msg->get(numFepFinal);
        msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
      }
//fepe

//soluteScaling
      if (simParams->soluteScalingOn) {
        delete [] ssAtomFlags;
        delete [] ss_vdw_type;
        delete [] ss_index;
        ssAtomFlags = new unsigned char[numAtoms];
        ss_vdw_type = new int [params->get_num_vdw_params()];
        ss_index = new int [numAtoms];
        msg->get(numAtoms*sizeof(unsigned char), (char*)ssAtomFlags);
        msg->get(ss_num_vdw_params);
        msg->get(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
        msg->get(numAtoms*sizeof(int), (char*)ss_index);
      }
//soluteScaling
#ifdef OPENATOM_VERSION
      // This needs to be refactored into its own version
      if (simParams->openatomOn) {
        delete [] fepAtomFlags;
        fepAtomFlags = new unsigned char[numAtoms];

        msg->get(numFepInitial);
        msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
#endif //OPENATOM_VERSION

      // DRUDE: receive data read from PSF
      msg->get(is_lonepairs_psf);
      if (is_lonepairs_psf) {
        msg->get(numLphosts);
        delete[] lphosts;
        lphosts = new Lphost[numLphosts];
        msg->get(numLphosts*sizeof(Lphost), (char*)lphosts);
      }
      msg->get(is_drude_psf);
      if (is_drude_psf) {
        delete[] drudeConsts;
        drudeConsts = new DrudeConst[numAtoms];
        msg->get(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
        msg->get(numTholes);
        delete[] tholes;
        tholes = new Thole[numTholes];
        msg->get(numTholes*sizeof(Thole), (char*)tholes);
        msg->get(numAnisos);
        delete[] anisos;
        anisos = new Aniso[numAnisos];
        msg->get(numAnisos*sizeof(Aniso), (char*)anisos);
      }
      msg->get(numZeroMassAtoms);
      // DRUDE

  //LCPO
  if (simParams->LCPOOn) {
    delete [] lcpoParamType;
    lcpoParamType = new int[numAtoms];
    msg->get(numAtoms, (int*)lcpoParamType);
  }

  //Receive GromacsPairStuff -- JLai

  if (simParams->goGroPair) {
    msg->get(numLJPair);
    delete [] indxLJA;
    indxLJA = new int[numLJPair];
    msg->get(numLJPair,indxLJA);
    delete [] indxLJB;
    indxLJB = new int[numLJPair];
    msg->get(numLJPair,indxLJB);
    delete [] pairC6;
    pairC6 = new Real[numLJPair];    
    msg->get(numLJPair,pairC6);
    delete [] pairC12;
    pairC12 = new Real[numLJPair];
    msg->get(numLJPair,pairC12);
    delete [] gromacsPair_type;
    gromacsPair_type = new int[numLJPair];
    msg->get(numLJPair,gromacsPair_type);
    delete [] pointerToLJBeg;
    pointerToLJBeg = new int[numAtoms];
    msg->get((numAtoms),pointerToLJBeg);
    delete [] pointerToLJEnd;
    pointerToLJEnd = new int[numAtoms];
    msg->get((numAtoms),pointerToLJEnd);
    // JLai
    delete [] gromacsPair;
    gromacsPair = new GromacsPair[numLJPair];
    for(int i=0; i < numLJPair; i++) {
        gromacsPair[i].atom1 = indxLJA[i];
        gromacsPair[i].atom2 = indxLJB[i];
        gromacsPair[i].pairC6  = pairC6[i];
        gromacsPair[i].pairC12 = pairC12[i];
        gromacsPair[i].gromacsPair_type = gromacsPair_type[i];
    }
    //
    msg->get(numGaussPair);
    delete [] indxGaussA;
    indxGaussA = new int[numGaussPair];
    msg->get(numGaussPair,indxGaussA);
    delete [] indxGaussB;
    indxGaussB = new int[numGaussPair];
    msg->get(numGaussPair,indxGaussB);
    delete [] gA;
    gA = new Real[numGaussPair];
    msg->get(numGaussPair,gA);
    delete [] gMu1;
    gMu1 = new Real[numGaussPair];
    msg->get(numGaussPair,gMu1);
    delete [] giSigma1;
    giSigma1 = new Real[numGaussPair];
    msg->get(numGaussPair,giSigma1);
    delete [] gMu2;
    gMu2 = new Real[numGaussPair];
    msg->get(numGaussPair,gMu2);
    delete [] giSigma2;
    giSigma2 = new Real[numGaussPair];
    msg->get(numGaussPair,giSigma2);
    delete [] gRepulsive;
    gRepulsive = new Real[numGaussPair];
    msg->get(numGaussPair,gRepulsive);
    delete [] pointerToGaussBeg;
    pointerToGaussBeg = new int[numAtoms];
    msg->get((numAtoms),pointerToGaussBeg);
    delete [] pointerToGaussEnd;
    pointerToGaussEnd = new int[numAtoms];
    msg->get((numAtoms),pointerToGaussEnd);
    //
  }
#endif

      //  Now free the message 
      delete msg;

#ifdef MEM_OPT_VERSION

      build_excl_check_signatures();

    //set num{Calc}Tuples(Bonds,...,Impropers) to 0
    numBonds = numCalcBonds = 0;
    numAngles = numCalcAngles = 0;
    numDihedrals = numCalcDihedrals = 0;
    numImpropers = numCalcImpropers = 0;
    numCrossterms = numCalcCrossterms = 0;
    numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
    // JLai
    numLJPair = numCalcLJPair = 0;
    // End of JLai

#else

      //  analyze the data and find the status of each atom
      build_atom_status();
      build_lists_by_atom();      

      
#endif
}
 /*      END OF FUNCTION receive_Molecule    */

/* BEGIN gf */
    /************************************************************************/
    /*                                                                      */
    /*      FUNCTION build_gridforce_params                                 */
    /*                                                                      */
    /*   INPUTS:                                                            */
    /*  gridfrcfile - Value of gridforcefile from config file               */
    /*  gridfrccol - Value of gridforcecol from config file                 */
    /*  gridfrcchrgcol - Value of gridforcechargecol from config file       */
    /*  potfile - Value of gridforcepotfile from config file                  */
    /*  initial_pdb - PDB object that contains initial positions            */
    /*  cwd - Current working directory                                     */
    /*                                                                      */
    // This function builds all the parameters that are necessary to
    // do gridforcing. This involves looking through a PDB object to
    // determine which atoms are to be gridforced, and what the force
    // multiplier is for each atom.  This information is then stored
    // in the arrays gridfrcIndexes and gridfrcParams.
    /************************************************************************/

void Molecule::build_gridforce_params(StringList *gridfrcfile,
                                      StringList *gridfrccol,
                                      StringList *gridfrcchrgcol,
                                      StringList *potfile,
                                      PDB *initial_pdb,
                                      char *cwd)
{
    PDB *kPDB;
    register int i;             //  Loop counters
    register int j;
    register int k;

    DebugM(3,  "Entered build_gridforce_params multi...\n");
//     DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
//     DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
    
    MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
    numGridforceGrids = 0;
    while (mgridParams != NULL) {
        numGridforceGrids++;
        mgridParams = mgridParams->next;
    }
    
    DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
    gridfrcIndexes = new int32*[numGridforceGrids];
    gridfrcParams = new GridforceParams*[numGridforceGrids];
    gridfrcGrid = new GridforceGrid*[numGridforceGrids];
    numGridforces = new int[numGridforceGrids];
    
    int grandTotalGrids = 0;    // including all subgrids
    
    mgridParams = simParams->mgridforcelist.get_first();
    for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
        int current_index=0;    //  Index into values used
        int kcol = 5;           //  Column to look for force constant in
        int qcol = 0;           //  Column for charge (default 0: use electric charge)
        Real kval = 0;          //  Force constant value retreived
        char filename[NAMD_FILENAME_BUFFER_SIZE];       //  PDB filename
        char potfilename[NAMD_FILENAME_BUFFER_SIZE];    //  Potential file name
        
        if (mgridParams == NULL) {
            NAMD_die("Problem with mgridParams!");
        }
        
        // Now load values from mgridforcelist object
        if (mgridParams->gridforceFile == NULL)
        {
            if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required.");
            kPDB = initial_pdb;
        }
        else
        {
            DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
            
            if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
            {
                strcpy(filename, mgridParams->gridforceFile);
            }
            else
            {
                strcpy(filename, cwd);
                strcat(filename, mgridParams->gridforceFile);
            }
        
            kPDB = new PDB(filename);
            if ( kPDB == NULL )
            {
                NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
            }
           
            if (kPDB->num_atoms() != numAtoms)
            {
                NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
            }
        }

        //  Get the column that the force constant is going to be in.  It
        //  can be in any of the 5 floating point fields in the PDB, according
        //  to what the user wants.  The allowable fields are X, Y, Z, O, or
        //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
        //  The default is the 5th field, which is beta (temperature factor)
        if (mgridParams->gridforceCol == NULL)
        {
            kcol = 5;
        }
        else
        {
            if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
            {
                kcol=1;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
            {
                kcol=2;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
            {
                kcol=3;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
            {
                kcol=4;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
            {
                kcol=5;
            }
            else
            {
                NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
            }
        }
    
        //  Get the column that the charge is going to be in.
        if (mgridParams->gridforceQcol == NULL)
        {
            qcol = 0;   // Default: don't read charge from file, use electric charge
        }
        else
        {
            if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
            {
                qcol=1;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
            {
                qcol=2;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
            {
                qcol=3;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
            {
                qcol=4;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
            {
                qcol=5;
            }
            else
            {
                NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
            }
        }
    
        if (kcol == qcol) {
            NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
        }

    
        //  Allocate an array that will store an index into the constraint
        //  parameters for each atom.  If the atom is not constrained, its
        //  value will be set to -1 in this array.
        gridfrcIndexes[gridnum] = new int32[numAtoms];
       
        if (gridfrcIndexes[gridnum] == NULL)
        {
            NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
        }
        
        //  Loop through all the atoms and find out which ones are constrained
        for (i=0; i<numAtoms; i++)
        {
            //  Get the k value based on where we were told to find it
            switch (kcol)
            {
            case 1:
                kval = (kPDB->atom(i))->xcoor();
                break;
            case 2:
                kval = (kPDB->atom(i))->ycoor();
                break;
            case 3:
                kval = (kPDB->atom(i))->zcoor();
                break;
            case 4:
                kval = (kPDB->atom(i))->occupancy();
                break;
            case 5:
                kval = (kPDB->atom(i))->temperaturefactor();
                break;
            }
           
            if (kval > 0.0)
            {
                //  This atom is constrained
                gridfrcIndexes[gridnum][i] = current_index;
                current_index++;
            }
            else
            {
                //  This atom is not constrained
                gridfrcIndexes[gridnum][i] = -1;
            }
        }
    
        if (current_index == 0)
        {
            //  Constraints were turned on, but there weren't really any constrained
            iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
        }
        else
        {
            //  Allocate an array to hold the constraint parameters
            gridfrcParams[gridnum] = new GridforceParams[current_index];
            if (gridfrcParams[gridnum] == NULL)
            {
                NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
            }
        }
    
        numGridforces[gridnum] = current_index;

        //  Loop through all the atoms and assign the parameters for those
        //  that are constrained
        for (i=0; i<numAtoms; i++)
        {
            if (gridfrcIndexes[gridnum][i] != -1)
            {
                //  This atom has grid force, so get the k value again
                switch (kcol)
                {
                case 1:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
                    break;
                case 2:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
                    break;
                case 3:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
                    break;
                case 4:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
                    break;
                case 5:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
                    break;
                }
            
                //  Also get charge column
                switch (qcol)
                {
                case 0:
#ifdef MEM_OPT_VERSION
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
#else
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
#endif
                    break;
                case 1:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
                    break;
                case 2:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
                    break;
                case 3:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
                    break;
                case 4:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
                    break;
                case 5:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
                    break;
                }
            }
        }
       
        //  If we had to create new PDB objects, delete them now
        if (mgridParams->gridforceFile != NULL)
        {
            delete kPDB;
        }
    
        //  Now we fill in our grid information
    
        // Open potential file
        if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
        {
            strcpy(potfilename, mgridParams->gridforceVfile);
        }
        else
        {
            strcpy(potfilename, cwd);
            strcat(potfilename, mgridParams->gridforceVfile);
        }
    
//        iout << iINFO << "Allocating grid " << gridnum
//             << "\n" << endi;
        
        DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n");
        gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams);
        
        grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
        DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi);
        
        // Finally, get next mgridParams pointer
        mgridParams = mgridParams->next;
    }
}
/* END gf */